数字信号处理--实验五 用DFT(FFT)对信号进行频谱分析.doc

学生实验报告开课学院及实验室： 电子楼317 2013 年 4 月 29 日学院机械与电气工程学院年级、专业、班姓名学号实验课程名称数字信号处理实验成绩实验项目名称实验五用DFT（FFT）对信号进行频谱分析指导老师一、实验目的学习DFT的基本性质及对时域离散信号进行频谱分析的方法，进一步加深对频域概念和数字频率的理解，掌握MATLAB函数中FFT函数的应用。二、实验原理离散傅里叶变换（DFT）对有限长时域离散信号的频谱进行等间隔采样，频域函数被离散化了，便于信号的计算机处理。 设x(n)是一个长度为M的有限长序列，x(n)的N点傅立叶变换：其中，它的反变换定义为： 令，则有： 可以得到，是Z平面单位圆上幅角为的点，就是将单位圆进行N等分以后第K个点。所以，X(K)是Z变换在单位圆上的等距采样，或者说是序列傅立叶变换的等距采样。时域采样在满足Nyquist定理时，就不会发生频谱混叠。DFT是对序列傅立叶变换的等距采样，因此可以用于序列的频谱分析。如果用FFT对模拟信号进行谱分析，首先要把模拟信号转换成数字信号，转换时要求知道模拟信号的最高截至频率，以便选择满足采样定理的采样频率。一般选择采样频率是模拟信号中最高频率的34倍。另外要选择对模拟信号的观测时间，如果采样频率和观测时间确定，则采样点数也确定了。这里观测时间和对模拟信号进行谱分析的分辨率有关，最小的观测时间和分辨率成倒数关系。最小的采样点数用教材相关公式确定。要求选择的采样点数和观测时间大于它的最小值。如果要进行谱分析的模拟信号是周期信号，最好选择观测时间是信号周期的整数倍。如果不知道信号的周期，要尽量选择观测时间长一些，以减少阶段效应的影响。用FFT对模拟信号作谱分析是一种近似的谱分析，首先一般模拟信号（除周期信号以外）的频谱是连续谱，而用FFT作谱分析得到的是数字谱，因此应该取FFT的点数多一些，用它的包络作为模拟信号的近似谱。另外，如果模拟信号不是严格的带限信号，会因为频谱混叠现象引起谱分析的误差，这种情况下可以预先将模拟信号进行滤波，或者尽量采样频率取高一些。一般频率混叠发生在折叠频率附近，分析时要注意因频率混叠引起的误差。最后要注意一般模拟信号是无限长的，分析时要截断，截断的长度和分辨率有关，但也要尽量取长一些，取得太短会截断引起的误差会很大。举一个极端的例子，一个周期性正弦波，如果所取观察时间太短，例如取小于一个周期，它的波形和正弦波相差太大，肯定误差很大，但如果取得长一些，即使不是周期的倍数，这种截断效应也会小一些。如同理论课教材所讨论的，在运用DFT进行频谱分析的时候可能有三种误差，即：（1）混叠现象当采样率不满足Nyquist定理，经过采样就会发生频谱混叠。这导致采样后的信号序列频谱不能真实的反映原信号的频谱。所以，在利用DFT分析连续信号频谱的时候，必须注意这一问题。避免混叠现象的唯一方法是保证采样的速率足够高，使频谱交叠的现象不出现。这告诉我们，在确定信号的采样频率之前，需要对频谱的性质有所了解。在一般的情况下，为了保证高于折叠频率的分量不会出现，在采样之前，先用低通模拟滤波器对信号进行滤波。（2）泄漏现象实际中的信号序列往往很长，甚至是无限长。为了方便，我们往往用截短的序列来近似它们。这样可以使用较短的DFT来对信号进行频谱分析。这种截短等价于给原信号序列乘以一个矩形窗函数。值得一提的是，泄漏是不能和混叠完全分离开的，因为泄漏导致频谱的扩展，从而造成混叠。为了减少泄漏的影响，可以选择适当的窗函数使频谱的扩散减到最小。（3）栅栏效应因为DFT是对单位圆上Z变换的均匀采样，所以它不可能将频谱视为一个连续函数。这样就产生了栅栏效应，从某种角度看， 用DFT来观看频谱就好像通过一个栅栏来观看一幅景象，只能在离散点上看到真实的频谱。这样的话就会有一些频谱的峰点或谷点被“栅栏”挡住，不能被我们观察到。减小栅栏效应的一个方法是在源序列的末端补一些零值，从而变动DFT的点数。这种方法的实质是改变了真是频谱采样的点数和位置，相当于搬动了“栅栏”的位置，从而使得原来被挡住的一些频谱的峰点或谷点显露出来。注意，这时候每根谱线所对应的频和原来的已经不相同了。从上面的分析过程可以看出，DFT可以用于信号的频谱分析，但必须注意可能产生的误差，在应用过程中要尽可能减小和消除这些误差的影响。DFT运算量较大，快速离散傅里叶变换算法FFT是解决方案。FFT并不是DFT不相同的另一种变换，而是为了减少DFT运算次数的一种快速算法。它是对变换式进行一次次的分解，使其成为若干小点数DFT的组合，从而减小运算量。常用的FFT是以2为基数的，其长度为。它的运算效率高，程序比较简单，使用也十分的方便。当需要进行变换的序列的长度不是2的整数次方的时候，为了使用以2为基的FFT，可以用末尾补零的方法，使其长度延长至2的整数次方。IFFT一般也可以通过FFT程序来完成。三、使用仪器、材料1、硬件：计算机 2、软件：Matlab四、实验步骤(一) 离散信号给定参考实验信号如下：用以8为周期进行周期性延拓形成的周期序列。(1) 分别以变换区间N8，16，32，对进行DFT(FFT)，画出相应的幅频特性曲线；(2) 分别以变换区间N4，8，16，对x2(n)分别进行DFT(FFT)，画出相应的幅频特性曲线；(3) 对x3(n)进行频谱分析，并选择变换区间，画出幅频特性曲线。（二）连续信号1. 实验信号：选择，式中的波形以及幅度特性如图7.1所示。式中频率自己选择。2. 分别对三种模拟信号选择采样频率和采样点数。对，选择采样频率，8，16，采样点数用.计算。对，周期，频率自己选择，采样频率，观测时间，采样点数用计算。图5.1R(t)的波形及其幅度特性对，选择采用频率，采样点数为16，32，64。3.分别对它们转换成序列，按顺序用表示。4.分别对它们进行FFT。如果采样点数不满足2的整数幂，可以通过序列尾部加0满足。5.计算幅度特性并进行打印。五、实验过程原始记录（数据、图表、计算等） (一) 离散信号%n=0:1:10;xn=ones(1,4),zeros(1,7); %输入时域序列向量xn=R4(n)Xk8=fft(xn,8); %计算xn的8点DFTXk16=fft(xn,16); %计算xn的16点DFTXk32=fft(xn,32); %计算xn的32点DFTk=0:7;wk=2*k/8; %产生8点DFT对应的采样点频率(关于归一化值)subplot(2,2,1);stem(n,xn,.);title(a) x_1 (n);xlabel(n);ylabel(x_1 (n);subplot(2,2,2);stem(wk,abs(Xk8),.); %绘制8点DFT的幅频特性图title(b) 8点DFT的幅频特性图);xlabel(omega/pi);ylabel(幅度);k=0:15;wk=2*k/16; %产生16点DFT对应的采样点频率(关于归一化值)subplot(2,2,3);stem(wk,abs(Xk16),.); %绘制16点DFT的幅频特性图title(c) 16点DFT的幅频特性图);xlabel(omega/pi);ylabel(幅度);k=0:31;wk=2*k/32; %产生32点DFT对应的采样点频率(关于归一化值)subplot(2,2,4);stem(wk,abs(Xk32),.); %绘制32点DFT的幅频特性图title(d) 32点DFT的幅频特性图);xlabel(omega/pi);ylabel(幅度);运行结果：%x2(n)n=0:50;xn=cos(pi/4*n); %输入时域序列向量cos(pi/4*n)Xk4=fft(xn,4); %计算xn的4点DFTXk8=fft(xn,8); %计算xn的8点DFTXk16=fft(xn,16); %计算xn的16点DFTk=0:3;wk=2*k/4; %产生4点DFT对应的采样点频率(关于归一化值)subplot(2,2,1);stem(n,xn,.); %绘制4点DFT的幅频特性图title(a) x_2 (n);xlabel(n);ylabel(x_2 (n);subplot(2,2,2);stem(wk,abs(Xk4),.); %绘制4点DFT的幅频特性图title(b) 4点DFT的幅频特性图);xlabel(omega/pi);ylabel(幅度);axis(0,2,0,2) k=0:7;wk=2*k/8; %产生8点DFT对应的采样点频率(关于归一化值)subplot(2,2,3);stem(wk,abs(Xk8),.); %绘制8点DFT的幅频特性图title(c) 8点DFT的幅频特性图);xlabel(omega/pi);ylabel(幅度);k=0:15;wk=2*k/16; %产生32点DFT对应的采样点频率(关于归一化值)subplot(2,2,4);stem(wk,abs(Xk16),.); %绘制16点DFT的幅频特性图title(d) 16点DFT的幅频特性图);xlabel(omega/pi);ylabel(幅度);运行结果：%x3(n)n=0:3;x=ones(1,4),zeros(1,4);xn=x(mod(n,8)+1); %输入xn=R4(n)以8为周期进行周期性延拓形成的周期序列Xk8=fft(xn,8); %计算xn的8点DFTXk16=fft(xn,16); %计算xn的16点DFTXk32=fft(xn,32); %计算xn的32点DFTk=0:7;wk=2*k/8; %产生8点DFT对应的采样点频率(关于归一化值)subplot(3,1,1);stem(wk,abs(Xk8),.); %绘制8点DFT的幅频特性图title(a) 8点DFT的幅频特性图);xlabel(omega/pi);ylabel(幅度);k=0:15;wk=2*k/16; %产生16点DFT对应的采样点频率(关于归一化值)subplot(3,1,2);stem(wk,abs(Xk16),.); %绘制16点DFT的幅频特性图title(c) 16点DFT的幅频特性图);xlabel(omega/pi);ylabel(幅度);k=0:31;wk=2*k/32; %产生32点DFT对应的采样点频率(关于归一化值)subplot(3,1,3);stem(wk,abs(Xk32),.); %绘制32点DFT的幅频特性图title(e) 32点DFT的幅频特性图);xlabel(omega/pi);ylabel(幅度);运行结果：（二）连续信号%x1(n)Fs1=4000; %采样频率Fs2=8000;Fs3=16000;Tp=0.002;N1=Fs1*Tp; %计算采样点数N2=Fs2*Tp;N3=Fs3*Tp;n1=0:N1-1;n2=0:N2-1;n3=0:N3-1;n=-0.0005:0.0001:0.002;xn=(n=0); %x1(n)Xk1=fft(xn,8); %计算xn的8点DFTk1=0:length(abs(Xk1)-1;wk1=2*k1/8; %产生8点DFT对应的采样点频率(关于归一化值)Xk2=fft(xn,16); %计算xn的16点DFTk2=0:length(abs(Xk2)-1;wk2=2*k2/16; %产生16点DFT对应的采样点频率(关于归一化值)Xk3=fft(xn,32); %计算xn的32点DFTk3=0:length(abs(Xk3)-1;wk3=2*k3/32; %产生32点DFT对应的采样点频率(关于归一化值)subplot(3,1,1);stem(wk1,abs(Xk1),.);xlabel(omega/pi);ylabel(幅度);title(X_1 (n) 8点傅立叶变换(f_s =4kHz);subplot(3,1,2);stem(wk2,abs(Xk2),.);xlabel(omega/pi);ylabel(幅度);title(X_2 (n) 16点傅立叶变换(f_s =8kHz);subplot(3,1,3);stem(wk3,abs(Xk3),.);xlabel(omega/pi);ylabel(幅度);title(X_1 (n) 32点傅立叶变换(f_s =16kHz);运行结果：%x2(n),取f=2kHzf=2000;Fs=4*f; %采样频率T=1/f;Tp1=0.5*T;Tp2=T;Tp3=2*T;Ts=1/Fs;N1=Fs*Tp1; %计算采样点数N2=Fs*Tp2;N3=Fs*Tp3;n1=0:N1-1;n2=0:N2-1;n3=0:N3-1;xn1=sin(2*pi*f*n1*Ts+pi/8); xn2=sin(2*pi*f*n2*Ts+pi/8);xn3=sin(2*pi*f*n3*Ts+pi/8);Xk1=fft(xn1,2);k1=0:length(abs(Xk1)-1;wk1=2*k1/2; %产生2点DFT对应的采样点频率(关于归一化值)Xk2=fft(xn2,4);k2=0:length(abs(Xk2)-1;wk2=2*k2/4; %产生4点DFT对应的采样点频率(关于归一化值)Xk3=fft(xn3,8);k3=0:length(abs(Xk3)-1;wk3=2*k3/8; %产生8点DFT对应的采样点频率(关于归一化值)subplot(3,1,1);stem(wk1,abs(Xk1),.);xlabel(omega/pi);ylabel(幅度);title(X_2 (n) 2点傅立叶变换(T_p =0.5T);subplot(3,1,2);stem(wk2,abs(Xk2),.);xlabel(omega/pi);ylabel(幅度);title(X_2 (n) 4点傅立叶变换(T_p =T);subplot(3,1,3);stem(wk3,abs(Xk3),.);xlabel(omega/pi);ylabel(幅度);title(X_2 (n) 8点傅立叶变换(T_p =2T);运行结果：%x3(n)Fs=64;T=1/Fs;N=64;n=0:N-1;x3n=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T);Xk16=fft(x3n,16);k1=0:length(abs(Xk16)-1;wk3=2*k1/8; %产生8点DFT对应的采样点频率(关于归一化值)Xk32=fft(x3n,32);k2=0:length(abs(Xk32)-1;wk2=2*k2/16; %产生16点DFT对应的采样点频率(关于归一化值)Xk64=fft(x3n,64);k3=0:length(abs(Xk64)-1;wk3=2*k3/64; %产生64点DFT对应的采样点频率(关于归一化值)subplot(3,1,1);stem(k1,abs(Xk16),.);xlabel(omega/pi);ylabel(幅度);title(X_3 (n) 16点傅立叶变换);subp