基于插值的非均匀信号的傅里叶变换算法研究
目录
一、 A/D卡设计 ....................................... 1
1.1 基于PCI 总线的A/D卡 . ......................... 1
1.2 基于USB 总线的A/D卡 . ......................... 2
二、 非均匀离散傅立叶变换 ............................. 4
三、 不同的插值算法 .................................. 6
1. 拉格朗日多项式插值 . ............................. 6
2. 三次样条插值 . ................................... 7
3. 牛顿插值........................................ 8
四、 主要算法及程序 ................................. 10
1. 拉格朗日算法 . .................................. 10
2. 三次样条插值 . .................................. 10
3. Newton算法 .................................... 12
五、 算法结果及比较分析 ............................. 14
六、 心得体会 ....................................... 19
七、 参考文献 ....................................... 20
一、A/D卡设计
1.1 基于PCI 总线的A/D卡
1、PCI 总线的含义
PCI 是由Intel 公司1991年推出的一种局部总线。从结构上看，PCI 是在CPU 和原来的系统总线之间插入的一级总线，具体由一个桥接电路实现对这一层的管理，并实现上下之间的接口以协调数据的传送。管理器提供了信号缓冲，使之能支持10种外设，并能在高时钟频率下保持高性能，它为显卡、声卡、网卡、MODEM 等设备提供了连接接口，它的工作频率为33MHz/66MHz。
PCI 是Peripheral Component Interconnect(外设部件互连标准）的缩写，它是目前个人电脑中使用最为广泛的接口，几乎所有的主板产品上都带有这种插槽。PCI 插槽也是主板带有最多数量的插槽类型，在目前流行的台式机主板上，ATX 结构的主板一般带有5～6个PCI 插槽，而小一点的MATX 主板也都带有2～3个PCI 插槽，可见其应用的广泛性。
PCI 总线是一种不依附于某个具体处理器的局部总线。管理器提供了信号缓冲，使之能支持10种外设，并能在高时钟频率下保持高性能。PCI 总线也支持总线主控技术，允许智能设备在需要时取得总线控制权，以加速数据传送。
图1.1 典型的PCI 系统总线构成
2、PCI 总线的基本含义
不同于ISA 总线，PCI 总线的地址总线与数据总线是分时复用的。这样做的好处是，
一方面可以节省接插件的管脚数，另一方面便于实现突发数据传输。在做数据传输时，由一个PCI 设备做发起者（主控，Initiator 或Master ），而另一个PCI 设备做目标（从设备，Target 或Slave ）。总线上的所有时序的产生与控制，都由Master 来发起。PCI 总线在同一时刻只能供一对设备完成传输，这就要求有一个仲裁机构（Arbiter ），来决定在谁有权力拿到总线的主控权。
当PCI 总线进行操作时，发起者（Master ）先置REQ#，当得到仲裁器（Arbiter ）的许可时（GNT#），会将FRAME#置低，并在AD 总线上放置Slave 地址，同时C/BE#放置命令信号，说明接下来的传输类型。所有PCI 总线上设备都需对此地址译码，被选中的设备要置DEVSEL#以声明自己被选中。然后当IRDY#与TRDY#都置低时，
可以传输数据。当Master 数据传输结束前，将FRAME#置高以标明只剩最后一组数据要传输，并在传完数据后放开IRDY#以释放总线控制权。
1.2 基于USB 总线的A/D卡
1.2.1 USB总线介绍
USB 总线为通用串行总线，USB 接口位于PS/2接口和串并口之间，允许外设在开机状态下热插拔，最多可串接下来127个外设，传输速率可达480Mb/S，P 它可以向低压设备提供5伏电源，同时可以减少PC 机I/O接口数量。USB 是基于通用连接技术，实现外设的简单快速连接，达到方便用户、降低成本、扩展PC 连接外设范围的目的。数据采集就是把来自各种传感器的信号数据实时地、准确地测量或汇集起来, 用计算机进行实时处理或记录存储, 实时完成测试和控制功能。数据采集系统结构通过微机的标准接口连接各种功能模块、仪器仪表和传感器, 组成测量系统。
1.2.2 USB接口电路设计
图1.2 USB接口电路结构图
R3是上拉电阻器，它可使USB 口的D+端上拉到DS2490S 的VB 端，表示USB 主机系
统是高速设备，同时这个上拉电阻器告诉主机有USB 设备插入。该上拉电阻器的设置对适配器的影响很大，它的负载值和1-Wire 网络的总长决定1-Wire 总线电压上升到5 V的速度。经过实验测试选择R3的阻值为27 Ω±lO%。R1、R2为USB 数据线保护电阻器。L 、L2具有禁止高频干扰并且减弱EMI 辐射的功能。LF33CV 为3.3 V 电压稳压器，与周围元件C1、C2组成强上拉部分，给EEPROM 或温度传感器等器件提供额外的电源。
1.2.3 USB接口的数据采集系统的设计实现
数据采集系统使用采集卡进行数据采集, 然后经过A/D转换器供计算机加工处理。基于USB 接口的数据采集与频谱分析系统本系统结构由硬件部分和软件部分组成, 硬件部分主要有计算机、I/ O 接口设备. 计算机作为硬件平台的核心可采用台式机, 系统采用的I/ O 设备为A/ D 数据采集卡, 该采集卡是一种基于USB 总线数据采集产品, 可与带USB 接口的各种台式计算机、笔记本电脑、工控机连接构成高性能的数据采集测量系统. 整个系统主要由4部分组成：USB 接口芯片及外围电路、控制电路、数据缓冲电路和A ／D 转换电路。USB 接口芯片选择了Cypress 公司的EZ-USB 2131Q，该芯片内嵌8051控制器，因此整个系统以EZ-USB 控制器为核心，由EZ-USB 经控制电路实现对A ／D 转换电路和数据缓冲电路的控制，模拟信号转换后的数据送入数据缓冲器，当数据缓冲器存满之后，通知EZ-USB 控制器，由主机取出数据。整个系统框图如图1.3所示。
图1.3 系统框图
1.2.4 A／D 转换电路
声卡是计算机对语音信号进行加工的重要部件，它具有对信号滤波、放大、采样保持、A/D和D/A转换等功能。系统中A ／D 转换芯片采用了MAXIM 公司的MAX122，该芯片是12 b 的高速的A ／D 转换器。在完全转换模式下，他的转换时间可以达到2．6μs，采样率为333 kS ／s 。MAX122有5种工作模式，在数据采集系统中，采用了模式2即连续转换模式。在这种模式下，每次转换需要13～14个时钟脉冲节拍，转换可以不间断地进行，但是需要提供开始转换使能信号，并且要保证使能信号和时钟信号同步，读信号和片选始终处于有效状态。数据输出使能信号一直有效，在转换结束时产生新的数据。
二、非均匀离散傅立叶变换
在过去的几十年里，随着DFT 的处理算法的发展，出现了对非均匀离散傅立叶变换的需求（NDFT ）。对NDFT 来说，数据在时域或者频域进行非均匀采样，，相应的采样点在N 个任意但不相同的位置。因此NDFT 可以看做是DFT 的一般形式，而DFT 是NDFT 的一种特殊情况。
NDFT 问题主要分为五大类：
第 1 类：由非等间隔分布的频域采样值估算均匀网格分布的空间域值；
第 2 类：由等间隔分布的频率值估算非等间隔指定点的空间域值；
第 3 类：由非等间隔分布的空间域值估算均匀网格分布的频域值；
第 4 类：由均匀分布的空间域值估算指定频率点的频域值；
第 5 类：由非均匀分布的空间域值估算非均匀分布的频域值。
弹光调制干涉信号的相位是非均匀分布的，要实现傅里叶变换，即要处理 NDFT 的第三类问题。因此，本章将以第 3 类 NDFT 问题为例分析 NDFT 算法。
设有连续的周期信号 x(t) ，其傅里叶变换 X ( f ) 为：
（2.1） T 为信号的周期；
将 x(t) 等相位采样，则有数据{xn} ， n = 0,1, 2,N -1，其傅里叶变换为 X n ( f ) 。均匀采样信号的离散傅里叶变换，则是将上式积分转换为求和相加的形式。
（2.2） Ts 为采样间隔，均匀采样时为常数。在计算频谱时是否引入采样间隔常数 Ts 都不影响频谱的检测。因此，均匀信号离散傅里叶变换表达式简化为：
（2.3） 同样，对于连续周期信号 x(t) 进行非等相位采样，其离散傅里叶变换表达式为
（2.4）
在 NDFT 中，每个采样段积分区间的宽度不相等。而采样信号各个频谱的大小和
采样间隔成比例关系。所以，在非均匀采样信号的傅里叶变换式中，需引入积分区间宽度
。
利用非均匀采样点拟合整个曲线，再对曲线进行均匀采样，经快速傅里叶变换重建。
三、不同的插值算法
1. 拉格朗日多项式插值
已知有 n+1 个数据点，满足
（3.1）
由这些数据点构造 n 次多项式函数，需满足条件：
（3.2）
则有
（3.3）
l k (x ) 称为 n 次插值的基函数。
于是，满足条件(2.3)的插值多项式 L n (x ) ，可表示为
（3.4） 由 l k (x ) 的定义可知，
（3.5）
公式(2.5)描述的插值多项式 Ln (x) ，称为拉格朗日(Lagrange)插值多项式。 拉格朗日插值多项式结构整齐，容易进行理论分析。然而在计算中，当插值点发生变化时基函数则需要重新计算，计算繁琐。这时可以用重心拉格朗日插值法或牛顿插值法来代替。在插值点比较多的时候，拉格朗日插值多项式的次数可能比较高，将引起龙格现象[，解决的办法是分段用低阶的插值多项式。在本项目研究中，插值的点数比较多，所以采用的是三次或四次分段插值多项式。
重心拉格朗日插值法是拉格朗日插值法的一种改进。设
（3.6）
可以将式(2.6)改写为：
（3.7）
定义重心权：
（3.8）
式(2.8)可简化为
（3.9）
式(2.9)可改写为：
（3.10）
式(2.10)描述了重心拉格朗日插值多项式。
2. 三次样条插值
已知有 n 个节点
(1) S(xi )=y I (i=0,1,2,…n) , 且有。
(2) S(x) 在每个小区间[xi , x i +1] 上是三次多项式；
(3) S(x) 在每个小区间[xi , x i +1] 上有连续的一阶和二阶导数；
则称 S(x) 为三次样条插值函数。由于 S(x) 在每个区间都是三次样条插值函数，共有 n 个区间，因此应确定 4n 个参数。
根据三次样条插值函数 S(x) 连续性，以及它的一阶、二阶导数的连续性、节点函数值可确定 4n-2个条件，因此还需确定两个条件。这两个条件一般选择区间的边界条件。边界条件一般有三种：
(1) 已知函数在两端点的导数值， S (x0 ) = f 0 ， S (xn ) = f n
(2) 已知函数在两端点的二阶导数值， S (x0 ) = f 0 ， S (xn ) = f n ，在特殊情况下，S '' '' '' '' ' ' (x0 ) = S (xn ) = 0 ''
(3) 当函数 S(x) 为周期函数时，有S(x0 ) = S(xn ) ， S (x0 ) = S (xn ) ， S (x0 ) = S (xn ) 三次样条插值算法的基本思想是在三次样条中，寻找三次多项式用来逼近每对数据点间的曲线。逼近两点间曲线的三次多项式可以有多条，为使得到的三次多项式唯一，对三次多项式做了条件限定。三次样条插值算法具有良好的收敛性、稳定性，以及二阶光滑性。
3. 牛顿插值
假设有n+1个不同的节点及函数在节点上的值（x 0，y 0），……（x n ，y n ），插值多项式有如下形式： ' ' '' ''
P （=α0+α（+α（(x -x 1) +⋯⋯+α（⋯⋯（x -x n ）n x ）1x -x 0）2x -x 0）n x -x 0）（x -x 1）
（3.11）
其中系数αi （i=0,1,2……n ）为特定系数，可由插值样条P （=y i （i=0,1,2……n ）n x i ）
确定。
根据均差的定义，把x 看成[a,b]上的一点, 可得
f(x)= f（x 0）+f[x 0，x 1]（x -x 0）
f[x, x 0]= f[x 0，x 1]+f[x,x 0，x 1] （x -x 1）
……
f[x, x 0，…x n -1]= f[x, x 0，…x n ]+ f[x, x 0，…x n ](x-xn ) (3.12) 综合以上式子，把后一式代入前一式，可得到：
f(x)= f（x 0）+f[x 0，x 1]（x -x 0）+ f[x 0，x 1，x 2]（x -x 0）（x -x 1）+ …+ f[x, x 0，…x n ]（x -x 0）…(x-xn )+ f[x, x 0，…x n ]ωn = Nn （x ）+R （其中 n x ）N n （x ）= f（x 0）+f[x 0，x 1]（x -x 0）+ f[x 0，x 1，x 2]（x -x 0）（x -x 1）+ …+ f[x, x 0，…x n ]（x -x 0）…(x-xn ) （3.13）
= f(x)- Nn （x ）= f[x, x 0，…x n ]ωn （2.14） R （n x ）
ωn =（x -x 0）…(x-xn )
Newton 插值的系数αi （i=0,1,2……n ）可以用差商表示。一般有
αk =f [x 0，x 1⋯⋯x k ] （k=0,1,2，……，n ） （2.15）
把（2.15）代入（1）得到满足插值条件N （（i=0,1,2，……n ）的n 次=f （x i ）n x i ）Newton 插值多项式
N n （x ）=f（x 0）+f[x 0，x 1]（x -x 1）+f[x 0，x 1，x 2]（x -x 1）（x -x 2）+……+f[x 0，x 1⋯⋯x n ]（x -x 1）（x -x 2）…（x -x n -1）.
其中插值余项为：
1
R （n x ）=f （x ）-N （x ）=f n +（ξ）
（n +1）！ωn +（1x ）
ξ介于x 0，x 1⋯⋯x k 之间。
四、主要算法及程序
1. 拉格朗日算法
function f = Language(x,y,x0) syms t ;
if (length(x) == length(y)) n = length(x); else
disp('x 和y 的维数不相等！' ); return ;
end %检错
f = 0.0; for (i = 1:n) l = y(i); for (j = 1:i-1)
l = l*(t-x(j))/(x(i)-x(j)); end ;
for (j = i+1:n)
l = l*(t-x(j))/(x(i)-x(j)); %计算拉格朗日基函数 end ;
f = f + l; %计算拉格朗日插值函数 simplify(f); %化简
if (i==n)
if (nargin == 3)
f = subs(f,'t' ,x0); %计算插值点的函数值 else
f = collect(f); %将插值多项式展开
f = vpa(f,6); %将插值多项式的系数化成6位精度的小数 end end end
2. 三次样条插值
function S=sancichazhi(x,y,t) n=length(x); z=zeros(n,4); z(:,1)=y; %生成前四列的差商表
for i=2:n for j=2:i if j==5 break end
z(i,j)=(z(i,j-1)-z(i-1,j-1))/(x(i)-x(i-j+1)); end end
%生成系数矩阵
b=linspace(2,2,n); for k=1:n if k==1
c(k)=-2; d(k,1)=-12*(x(2)-x(1))*z(4,4); elseif k==n
a(k)=-2; d(k,1)=12*(x(n)-x(n-1))*z(k,4); c(k)=0; else
d(k,1)=6*z(k+1,3);c(k)=(x(k+1)-x(k))/(x(k+1)-x(k-1)); a(k)=1-c(k); end end
H=[a;b;c;d'];
u(1)=b(1); r(1)=c(1); y1(1)=d(1); for k=2:n r(k)=c(k); l(k)=a(k)/u(k-1); u(k)=b(k)-l(k)*r(k-1); y1(k)=d(k)-l(k)*y1(k-1); end
m(n)=y1(n)/u(n); for k=n-1:-1:1
m(k)=(y1(k)-r(k)*m(k+1))/u(k); end
a(1)=[];c(n)=[];
a1=diag(a,-1);b1=diag(b);c1=diag(c,1); M=a1+b1+c1;
%整理三次多项式 for i=2:n
h(i)=x(i)-x(i-1);b1=[x(i),x(i),x(i)];b2=[x(i-1),x(i-1),x(i-1)]; a1=-m(i-1)/(6*h(i))*poly(b1); a2=m(i)/(6*h(i))*poly(b2);
a33=(m(i-1)*(h(i))^2/6-y(i-1))/h(i)*poly(x(i)); a3=[0,0,a33]; a44=-(m(i)*(h(i))^2/6-y(i))/h(i)*poly(x(i-1)); a4=[0,0,a44]; S(i-1,:)=a1+a2+a4+a3; end
%对于给定的向量，求出三次样条插值函数在某些点的值 l=length(t); for i=1:l for j=2:n if
x(1)+(x(n)-x(1))/(n-1)*(j-2)
disp(' 结果为三次样条插值分段函数的降幂系数' ) plot(t,R,'o' ,x,y, 'kd' ) title(' 三次样条插值' ) grid on end
3. Newton 算法
function f = Newton(x,y,x0) syms t ;
if (length(x) == length(y)) n = length(x); c(1:n) = 0.0; else
disp('x 和y 的维数不相等！' ); return ; end
f = y(1); y1 = 0; l = 1;
for (i=1:n-1) for (j=i+1:n)
y1(j) = (y(j)-y(i))/(x(j)-x(i)); end
c(i) = y1(i+1); l = l*(t-x(i));
f = f + c(i)*l; simplify(f); y = y1;
if (i==n-1)
if (nargin == 3) f = subs(f,'t' ,x0); else
f = collect(f); %将插值多项式展开 f = vpa(f, 6); end end end
调用程序如下 N=1024;
x=[0.5,1.8,2.4,3.6,4.7,5.1,6.5,7.3,8.2,9.8]
y=[0.30886552,0.904582781,0.997978435,0.771243677,0.188851681,-0.061169136,-0.807798281,-0.991820585,-0.90593608,-0.128429604] y0=zeros(1,101); x0=[0:0.1:10]; for i=1:101
y0(i)=Language(x,y,x0(i)) end figure; plot(x,y); figure; plot(x0,y0);
figure;
Y = fft(y0,N); %做FFT 变换 Ayy = (abs(Y)); %取模
Ayy=Ayy/(N/2); %换算成实际的幅度 Ayy(1)=Ayy(1)/2;
F=([1:N]-1)*10/N; %换算成实际的频率值
plot(F(1:N/2),Ayy(1:N/2)); %显示换算后的FFT 模值结果 title(' 幅度-频率曲线图' );
五、算法结果及比较分析
1
0.5
-0.5
-1
0246810
a
1
0.5
-0.5
-1
0246810
b
幅度-频率曲线图
0.120.10.080.060.040.020
c
图5.1 拉格朗日插值结果
上图中，我们对sin （π/5t）进行非均匀采样10个点，如图a 。然后对非均匀信号进行插值，插100个点，如图b 。再对插值后的均匀信号进行fft 变换得到其频谱，如图c 。
1
0.5
-0.5
-1
a
1
0.5
-0.5
-1
B
幅度-频率曲线图
0.120.10.080.060.040.020
012
345
图5.2 三次样条插值结果
1
0.5
-0.5
-1
a
1
0.5
-0.5
-1
0246810
b
幅度-频率曲线图
0.120.10.080.060.040.020
c
图5.3 牛顿插值结果
再对刚刚信号进行牛顿插值，其结果如上图5.3所示。
拉格朗日插值算法的公式结构整齐紧凑，在理论分析是十分方便，然而在实际操作中，插值点每增加或减少时，所对应的基本多项式就需要全部计算，于是整个计算公式都会改变，十分繁琐。这时可用牛顿法。此外，当插值点比较多时，多项式次数会比较高，数值会不稳定。复杂度是O （n^2）。
牛顿插值公式是n 次插值多项式的又一种构造，它克服拉格朗日差值的缺点，每增加一个插值节点，只要在原牛顿插值公式增添一项可形成高一次的插值公式。但不能全面反映被插值函数的特性不满足实际问题要求的插值函数在各点与被插值函数相同，或者与被插函数的导数值也相等，这时牛顿插值便不符合要求。
三次样条插值鉴于高次插值不收敛又不稳定的特点, 低次插值既具有收敛性又具有稳定性, 因此低次值更具有实用价值, 但是低次插值的光滑性较差。
六、心得体会
通过这三个星期的课程设计，我学到了很多的东西, 不仅巩固了我以前所学过的知识, 还让我学到很多在书本上所没有学到过的知识。
同时进一步加深了对插值算法的了解和熟练了对Matlab 的使用, 让我对信号处理这门课程有了更加浓厚的兴趣。 因为以前都是基于课本上所学的理论知识，然而通过这次课程设计之后才能真正理解其意义。
在这次课程设计的过程中，我做的部分是插值算法编写，在课程设计开始部分，我遇到不少的问题, 比如刚开始, 无法调用函数，经过后来的改正才可以。还有刚开始由于对非均匀傅立叶变换原理并不是很了解, 不再单单只是懂理论, 理论与实际相结合是很重要的, 只有理论知识是远远不够的, 只有把所学的理论知识与实践相结合起来, 从理论中得出结论。总的来说, 通过这次的课程设计我对语音信号有了全面的认识, 对Matlab 的知识又有了深刻的理解, 让我感受到只有在充分理解课本 知识的前提下, 才能更好的应用这个工具。
这次课程设计使我了解了MATLAB 的使用方法, 学会分析滤波器的优劣和性能, 提高了分析和动手实践能力，同时我相信, 进一步加强对 MATLAB 的学习与研究对我今后的学习将会起到很大的帮助
七、参考文献
[1]张敏娟，王召巴，王志彬. 弹光调制傅里叶变换光谱复原高速数据处理技术.2013
[2] 丛玉良, 王宏志. 数字信号处理原理及其MATLAB 实现(第2版). 北京：电子工业出版社，2009
[3] 张兰勇，孙健，孙晓云等.LabVIEW 程序设计基础与提高. 北京：机械工业出版社，2012
[4] 雷学堂，徐火希. 可直接感受的基于MATLAB 的语音滤波. 合肥学院学报， 2006
[5] 徐明远, 刘增力.MATLAB 仿真在信号处理中的应用. 西安：西安电子科技大学出版社,2007
20