离散时间信号通过线性时不变系统

数 字 信 号 处 理

实 验 报 告

实验名称： 离散时间信号通过线性时不变系统 姓 名： 专 业： 年 级： 学 号： 指导教师：

1

一、实验目的 1、通过本实验，进一步加深对DFT算法原理和基本性质的理解，熟悉FFT算法原理和FFT子程序的应用。 2、掌握应用FFT对信号进行频谱分析的方法。 3、通过本实验进一步掌握频谱采样定理。 4、了解应用FFT进行信号频谱分析过程中可能出现的问题，以便在实际中正确应用FFT。 二、实验原理及方法 1、一个连续时间信号xa(t)的频谱可以用它的傅立叶变换表示 Xa(j)2、对信号进行理想采样，得到采样序列 xa(t)ejtdt x(n)xa(nT) 3、以T为采样周期，对x(n)进行Z变换 X(z)x(n)zn 4、当zej时，得到序列傅立叶变换SFT X(e)x(n)ejn j5、为数字角频率 T6、已经知道：  Fs12m X(e)Xa[j()] ( 2-6 ) TTTj7、序列的频谱是原模拟信号的周期延拓，即可以通过分析序列的频谱，得到相应连续信号的频谱。（信号为有限带宽，采样满足Nyquist定理） 8、无线长序列可以用有限长序列来逼近，对于有限长序列可以使用离散傅立叶变换（DFT）。可以很好的反映序列的频域特性，且易于快速算法在计算机上实现。 当序列x(n)的长度为N时，它的离散傅里叶变换为： X(k)DFT[x(n)]j2Nx(n)Wn0N1knN 其中WNe，它的反变换定义为： 1N1kn x(n)IDFT[X(k)]X(k)WN Nk0 2

比较Z变换式 ( 2-3 ) 和DFT式 ( 2-7 )，令zWN则 X(z)k |kZWNnkx(n)WNDFT[x(n)] N0N1因此有X(k)X(z)|kzWN k WN是Z平面单位圆上幅角为2k的点，也即是将单位圆N等分后的第k点。所以X(k)是x(n)的Z变N换在单位圆上等距采样，或者说是序列傅里叶变换的等距采样。 9、DFT是对序列傅里叶变换的等距采样，因此可以用于序列的频谱分析。在运用DFT进行频谱分析的过程中有可能产生三种误差，这里给出三种误差的定性讨论。 三种误差：混叠现象、泄露现象、栅栏效应 1) 混叠现象 ( 2-6 )式说明序列的频谱是原模拟信号的频谱的周期延拓，周期为21。因此当采样频率Fs小于两倍信号TT（这里指是信号）最大频率时，经过采样就会发生频谱混叠，这使得采样后的信号序列频谱不能真实地反映原信号的频谱。所以在利用DFT分析连续信号的频谱时，必须注意这一问题。避免混叠现象的唯一方法是保证采样速率足够高，使频谱交叠现象不致出现。也就是说，在确定采样频率之前，必须对信号的性质有所了解，一般在采样前，信号通过一个防混叠低通滤波器。 2) 泄漏现象 实际中的信号序列往往很长，为了方便我们往往用截短的序列来近似它们，这样可以使用较短的DFT来对信号进行频谱分析，这种截短等价于给原信号序列乘以一个矩形窗函数。 泄漏是不能与混叠完全分离开的，因为泄漏导致频谱的扩散，从而造成混叠，为了减小泄漏的影响，可以选择适当的窗函数，是频谱的扩散减到最小。 3) 栅栏效应 因为DFT是对单位圆上Z变换的均匀采样，所以他不可能将频谱视为一个连续函数。这样就产生了栅栏效应，就一定意义上看，DFT来观看频谱就好像通过一个尖桩的栅栏来观看一个图景一样，只能在离散点上看到真实频谱，这样就可能发生一些频谱的峰点或谷点被“尖桩的栅栏”所挡住，不能被我们观察到。减小栅栏效应的一个方法就是借助在原序列的末端添补一些零值，从而变动DFT的点数。这一方法实际上是人为地改变了对真实谱采样的点数和位置，相当于搬动了每一根“尖桩栅栏”的位置，从而使得频谱的峰点或者谷点暴露出来。当然，这是每根谱线所对应的频率和原来的不同了。 综上所述，DFT可以用于信号的频谱分析，但必须注意可能产生的误差，在应用过程中要尽可能减少和消除这些误差的影响。 M 快速傅里叶变换FFT是为了减少DFT运算次数的一种算法，常用的FFT是以2为基数的，其长度N2，它的效率高，程序简单，使用也十分方便，当要变换的序列长度不等于2的正整数次方时，可以用末尾补零的方法，使其长度延长到2的整数次方。 三、实验结果及分析 注：图均为上面时域波形，下面频域波形 1、 高斯序列 exp[(np)2/qxa(n)0

0nN1 else3

代码1: /************FFT***********/ #include #include #include #include #define N 1000 typedef struct { double real; double img; }complex; void fft(); /*快速傅里叶变换*/ void initW(); void change(); void add(complex ,complex ,complex *); /*复数加法*/ void mul(complex ,complex ,complex *); /*复数乘法*/ void sub(complex ,complex ,complex *); /*复数减法*/ void divi(complex ,complex ,complex *);/*复数除法*/ void output(); /*输出结果*/ complex x[N], *W;/*输出序列的值*/ int size_x=0;/*输入序列的长度，只限2的N次方*/ double p=16,q=2,PI;double z[N]; int main() { int i; system(\"cls\"); PI=atan(1)*4; size_x=512; for(i=0;i<32;i++) { x[i].real=exp(-(i-p)*(i-p)/q);x[i].img=0;} for(i=32;i4

{ l=1<0 ) { j=j<<1; j|=(k & 1); k=k>>1; } if(j>i) { temp=x[i]; x[i]=x[j]; x[j]=temp; } } } void output() /*输出结果*/ { int i; double y[N];

5

for(i=0;ireal=a.real+b.real; c->img=a.img+b.img; } void mul(complex a,complex b,complex *c) { c->real=a.real*b.real - a.img*b.img; c->img=a.real*b.img + a.img*b.real; } void sub(complex a,complex b,complex *c) { c->real=a.real-b.real; c->img=a.img-b.img; } void divi(complex a,complex b,complex *c) { c->real=( a.real*b.real+a.img*b.img )/( b.real*b.real+b.img*b.img); c->img=( a.img*b.real-a.real*b.img)/(b.real*b.real+b.img*b.img); } P=16,N=32,q=2,FFT点数为512 6

P=16,N=32,q=10 FFT点数为512 P=16,N=32,q=30 FFT点数为512 时域：q取值的增大，信号波形变宽，变矮，在最大值处过度变的平缓。 频域：信号的频谱向低频靠近。 方差q=2 时，信号变化相对快，高频分量大。 方差q=30时，信号变化相对慢，低频分量大。 因为随着q取值的增大，高斯信号逐渐变得平缓，过渡带变得平滑并延长，从而低频分量增加，高频分量减少。 7

P=25,N=32,q=10 FFT点数为512 P=30,N=32,q=10 FFT点数为512 P=32,N=32,q=10 FFT点数为512 时域：p取值的增大，信号波形逐渐向右平移。 频域：信号的频谱中高频分量逐渐增加，频谱泄漏逐渐明显，并逐渐出现频谱混叠现象。当p=32时，能力泄漏至旁边的频率，出现较明显的频谱泄漏与频谱混叠现象。 随着p值增大，信号被截断部分增多，截断部分的过渡带过陡，产生高频分量增多，而造成频谱泄漏与混叠。

8

2、 正弦序列 sin(2fn)xb(n)00nN1 else代码 /************FFT***********/ #include #include #include #include #define N 1000 typedef struct { double real; double img; }complex; void fft(); /*快速傅里叶变换*/ void initW(); void change(); void add(complex ,complex ,complex *); /*复数加法*/ void mul(complex ,complex ,complex *); /*复数乘法*/ void sub(complex ,complex ,complex *); /*复数减法*/ void divi(complex ,complex ,complex *);/*复数除法*/ void output(); /*输出结果*/ complex x[N], *W;/*输出序列的值*/ int size_x=0;/*输入序列的长度，只限2的N次方*/ double f ,PI,z[N]; int main() { int i; system(\"cls\"); PI=atan(1)*4; size_x=32; f=0.0625; for(i=0;i<32;i++) { x[i].real=sin(2*PI*f*i);x[i].img=0;} for(i=32;i9

change(); for(i=0;i0 ) { j=j<<1; j|=(k & 1); k=k>>1; } if(j>i) { temp=x[i]; x[i]=x[j]; x[j]=temp; } } } void output() /*输出结果*/

10

{ int i; double y[N]; for(i=0;i<1000;i++) y[i]=sqrt(x[i].real*x[i].real+x[i].img*x[i].img ); initgraph(1200,600); for(i=0;ireal=a.real+b.real; c->img=a.img+b.img; } void mul(complex a,complex b,complex *c) { c->real=a.real*b.real - a.img*b.img; c->img=a.real*b.img + a.img*b.real; } void sub(complex a,complex b,complex *c) { c->real=a.real-b.real; c->img=a.img-b.img; } void divi(complex a,complex b,complex *c) { c->real=( a.real*b.real+a.img*b.img )/( b.real*b.real+b.img*b.img); c->img=( a.img*b.real-a.real*b.img)/(b.real*b.real+b.img*b.img); } f=0.0625,N=32,FFT点数32 当FFT点数为32时，频谱为单线谱，只在谱峰处有值，其他位置都为0。

11

f=0.0625,N=32,FFT点数512 FFT点数为512时除谱峰以外，其他位置也有值。 出现这种现象是由于栅栏效应引起的，导致采样时只采到谱峰与零值点。利用频谱估计频率时，f峰的位置正确。 mN，m为谱峰的位置，估计值与实际值一致，所以谱 f=0.265625,N=32,FFT点数32

12

f=0.265625,N=32,FFT点数64 f=0.265625,N=64,FFT点数32 f=0.265625,N=64,FFT点数64 N=FFT点数=32、64时没有出现单线谱 N=FFT点数=64的时候出现单线谱 因为当点数为32时，FFT对频域采样点没有采到谱峰位置，而有一定的相位差，其他点采到了各个旁瓣的峰值。而当点数为64点时，正好采样采到谱峰和旁瓣的零点。要使频谱正好采到谱峰，满足此处，f以。 fkFsN。 0.2656251764，所以64点FFT可以采到谱峰，而32点FFT不可 13

f=0.245,N=256 （5）f=1.96kHz 采样频率Fs=8kHz N=256 此时频谱分辨率是多少? 通过FFT离散谱观察到的的信号模拟频率与实际频率相差多少？ 1、当f=0.245时，正弦序列的时域波形发生了假调制现象，这是因为时域点数为2（m0,1,2)，而只有当mfk(m,k为整数）时，对正弦函数的采样才能在每个周期内采到最大值点，从而出现等幅波。 m2 而当f0.245时不能写成这种方式，造成在一个周期内采样时与最大值点有一定的相位差，由于相位的累积，从而造成每个周期内采样最大值的周期性变化， 从整体上看即呈现出假调制现象。 在这种情况下，由于FFT点数为2（m0,1,2)，因此无法使频谱出现单载波。 2、●当利用64点的FFT估计频率时， 频谱分辨率为FmFs80.125kHz N64kfN0.2456415.6816 实际信号频率为 Fs8fk162kHzN64 误差为 ff0.04kHz ●当利用256点的FFT估计频率时： 频谱分辨率为FFs80.03125kHz N256kfN0.24525662.7263 实际信号频率为 Fs8fk631.96875kHzN256 误差为 ff0.00875kHz。 ●当利用512点的FFT估计频率时： 频谱分辨率为FFs80.015625kHz N512kfN0.245512125.44126 实际信号频率为 Fs8fk1261.96875kHzN512

14

误差为 ff0.00875kHz。 结论： （1）FFT的点数越多，信号的频谱分辨率越高，利用频谱估计得到的信号频率 与实际的误差越小。 （2）且在增加FFT点数到512时，频谱分辨率增加，但是实际信号频率与模拟频率误差相等，说明已达取该FFT点数范围内的最大精度。 3、 衰减正弦信号 ensin(2fn)xc(n)00nN1 else a=0.01 f=0.21875 N=32 FFT点数32 a=0.01 f=0.4375 N=32 FFT点数32 15

a=0.01 f=0.5625 N=32 FFT点数32 a=0.01 f=0.21875 N=32 FFT点数256 a=0.01 f=0.4375 N=32 FFT点数256 16

a=0.01 f=0.5625 N=32 FFT点数256 满足Nyquist定理时 ff0 Fs2f0，即f0.5 Fsf=0.5625时不满足Nyquist定理。 随着f的增大，频谱的谱峰逐渐向右平移，两谱峰逐渐向中间靠拢。 因为0.43750.50.0625，0.56250.50.0625， f=0.4375和f=0.5625频谱图关于对称 造成观察到的频谱完全相同，但实际上表示的意义却不相同。 f=0.4375时的谱峰位于n=113处， f=0.5625时的谱峰位于n=143处。 由于存在泄漏现象，出现了高频分量，虽然在f=0.4375时满足Nyquist定理，但实际上已发生了频谱混叠。 4、三角波序列和反三角序列波 （1）三角序列波 n1xd(n)8n00n34n7 else 17

FFT点数8 FFT点数256 （2）反三角序列 4nxd(n)n300n34n7 else 18

FFT点数8 FFT点数为256 FFT点数=8，虽然两者的时域波形不同，但是频域波形却相同，因为二者满足循环移位关系，即xe(n)xd((n4))8R8(n)，从而Xd(k)Xe(k)，这种现象是栅栏效应引起的。 FFT点数=256，对两序列后续补零，导致二者不再满足循环移位关系，Xd(k)Xe(k)，所以频谱不相同。 四、 实验总结 （1）利用FFT来估计模拟信号的频率，FFT的点数越多，信号的频谱分辨率越高，利用频谱估计得到的信号频率 与实际的误差越小。要想提高估计精度，应当使FFT 19

点数在允许的范围内尽可能的大。 （2）当信号发生截断效应时，会产生高频分量，这种情况下会出现频谱混叠和频谱泄漏现象，截断效应越明显，泄漏越大。 （3）对信号抽样时，在满足Nyquist定理的情况下，由于频谱泄漏，会发生频谱混叠，因此，采样频率尽量高。 （4）由于栅栏效应使得有些不同的信号的频谱图相同，这种情况可以通过增加FFT点数来判断。 （5）时域波形的变化情况可以看出频域能量的分布！ 指导教师： 年 月 日 教师评语 20