1、 学生实验报告 开课学院及实验室： 电子楼 317 2013 年 4 月 29 日 学院 机械与电气 工程学院 年级、专 业、班 姓名 学号 实验课程名称 数字信号处理实验 成绩 实验项目名称 实验五 用 DFT（ FFT）对信号进行频谱分析 指导老师 一、实验目的 学习 DFT 的基本性质及对时域离散信号进行频谱分析的方法，进一步加深对频域概念和数字频率的理解，掌握 MATLAB 函数中 FFT 函数的应用。 二、实验 原理 离散傅里叶变换（ DFT）对有限长时域离散信号的频谱进行等间隔采样 ，频域函数被离散化了，便于信号的计算机处理。 设 x(n)是一个长度为 M 的有限长序列， x(n)
2、的 N 点傅立叶变换： 其中 NjN eW 2 ，它的反变换定义为： 10 )(1)( Nk nkNWkXNnx 令 kNWz ，则有： 10 )()(NnnkNkNWz WnxzX 可以得到， kNWzzXkX )()( ， kNWz 是 Z 平面单位圆上幅角为 kN 2 的点，就是将单位圆进行 N 等分以后第 K 个点。所以， X(K)是 Z 变换在单位圆上的等距采样，或者说是序列傅立叶变换的等距采样。时域采样在满足 Nyquist 定理时，就不会发生频谱混叠。 DFT 是对序列傅立叶变换的等距采样，因此可以用于序列的频谱分析。 如果 用 FFT 对模拟信号进行谱分析，首先要把模拟信号转换
3、成数字信号，转换时要求知道模拟信号的最高截至频率，以便选择满足采样定理的采样频率。一般选择采样频率是模拟信号中最高频率的 34 倍。另外要选择对模拟信号的观测时间，如果采样频率和观测时间确定，则采样点数也确定了。这里观测时间和对模拟信号进行谱分析的分辨率有关，最小的观测时间和分辨率成倒数关系。最小的采样点数用教材相关公式确定。要求选择的采样点数和观测时间大于它的最小值。 如果要进行谱分析的模拟信号是周期信号，最好选择观测时间是信号周期的整数倍。如果不知道信号的 周期，要尽量选择观测时间长一些，以减少阶段效应的影响。 用 FFT 对模拟信号作谱分析是一种近似的谱分析，首先一般模拟信号（除周期信号
4、以外）的频谱是连续谱，而用 FFT 作谱分析得到的是数字谱，因此应该取 FFT 的点数多一些，用它的包络作为模拟信号的近似谱。另外，如果模拟信号不是严格的带限信号，会因为频谱混叠现象引起谱分析的误差，这种情况下可以预先将模拟信号进行滤波，或者尽量采样频率取高一些。 一般频率混叠发生在折叠频率附近，分析时要注意因频率混叠引起的误差。最后要注意一般模拟信号是无限长的，分析时要截断，截断的长度 和分辨率有关，但也要尽量取长一些，取得太短会截断引起的误差会很大。举一个极端的例子，一个周期性正弦波，如果所取观察时间太短，例如取小于一个周期，它的波形和正弦波相差太大，肯定误差很大，但如果取得长一些，即使不
5、是周期的倍数，这种截断效应也会小一些。 如同理论课教材所讨论的，在运用 DFT 进行频谱分析的时候可能有三种误差，即： （ 1）混叠现象 当采样率不满足 Nyquist 定理，经过采样就会发生频谱混叠。这导致采样后的信号序列频谱不能真实的反映原信号的频谱。所以，在利用 DFT 分析连续信号频谱的时候，必须注意这一问题 。避免混叠现象的唯一方法是保证采样的速率足够高，使频谱交叠的现象不出现。这告诉我们，在确定信号的采样频率之前，需要对频谱的性质有所了解。在一般的情况下，为了保证高于折叠频率的分量不会出现，在采样之前，先用低通模拟滤波器对信号进行滤波。 （ 2）泄漏现象 实际中的信号序列往往很长，
6、甚至是无限长。为了方便，我们往往用截短的序列来近似它们。这样可以使用较短的 DFT 来对信号进行频谱分析。这种截短等价于给原信号序列乘以一个矩形窗函数。值得一提的是，泄漏是不能和混叠完全分离开的，因为泄漏导致频谱的扩展，从而造成混叠。为了 减少泄漏的影响，可以选择适当的窗函数使频谱的扩散减到最小。 （ 3）栅栏效应 因为 DFT 是对单位圆上 Z 变换的均匀采样，所以它不可能将频谱视为一个连续函数。这样就产生了栅栏效应，从某种角度看， 用 DFT 来观看频谱就好像通过一个栅栏来观看一幅景象，只能在离散点上看到真实的频谱。这样的话就会有一些频谱的峰点或谷点被“栅栏”挡住，不能被我们观察到。减小栅
7、栏效应的一个方法是在源序列的末端补一些零值，从而变动 DFT 的点数。这种方法的实质是改变了真是频谱采样的点数和位置，相当于搬动了“栅栏”的位置，从而使得原来被挡住的一些 频谱的峰点或谷点显露出来。注意，这时候每根谱线所对应的频和原来的已经不相同了。 从上面的分析过程可以看出， DFT 可以用于信号的频谱分析，但必须注意可能产生的误差，在应用过程中要尽可能减小和消除这些误差的影响。 DFT 运算量较大，快速离散傅里叶变换算法 FFT 是解决方案 。 FFT 并不是 DFT 不相同 的另一种变换，而是为了减少 DFT 运算次数的一种快速算法。它是对变换式进行一次次的分解，使其成为若干小点数 DF
8、T 的组合，从而减小运算量。常用的 FFT 是以 2 为基数的，其长度为 MN 2 。它的运算效率高，程序比较简单，使用也十分的方 便。当需要进行变换的序列的长度不是 2 的整数次21 j0( ) ( ) ( ) e 0 1N knNN nX k D F T x n x n k N 方的时候，为了使用以 2 为基的 FFT，可以用末尾补零的方法，使其长度延长至 2 的整数次方。 IFFT一般也可以通过 FFT 程序来完成。 三、使用仪器、材料 1、硬件：计算机 2、软件： Matlab 四、实验步骤 (一 ) 离散信号 给定参考实验 信号如下： 3()xn：用 14( ) ( )x n R n
9、 以 8 为周期进行周期性延拓形成的周期序列。 (1) 分别以变换区间 N 8， 16， 32，对 14( ) ( )x n R n 进行 DFT(FFT)，画出相应的幅频特性曲线； (2) 分别以变换区间 N 4， 8， 16， 对 x2(n)分别进行 DFT(FFT)，画出相应的幅频特性曲线； (3) 对 x3(n)进行频谱分析，并选择变换区间，画出幅频特性曲线。 （二）连续信号 1. 实验信号： 1( ) ( )x t R t 选择 1.5ms ，式中 ()Rt 的波形以及幅度特性如图 7.1 所示。 2 ( ) s in ( 2 / 8 )x t ft 式中频率 f 自己选择。 3 (
10、 ) c o s 8 c o s 1 6 c o s 2 0x t t t t 2. 分别对三种模拟信号选择采样频率和采样点数。 对 1()xt ()Rt ，选择采样频率 4sf kHz ， 8kHz ， 16kHz ，采样点数用 . sf 计算。 对 2 ( ) s in ( 2 / 8 )x t ft，周期 1/Tf ，频率 f 自己选择，采样频率 4sff ，观测时间0.5pTT ， T ， 2T ，采样点数用 psTf计算。 图 5.1 R(t)的波形及其幅度特性 对 3 ( ) c o s 8 c o s 1 6 c o s 2 0x t t t t ，选择采用频率 64sf Hz
11、，采样点数为 16， 32， 64。 3. 分别对它们转换成序列，按顺序用 1 2 3( ), ( ), ( )x n x n x n表示。 4. 分别对它们进行 FFT。如果采样点数不满足 2 的整数幂，可以通过序列尾部加 0 满足。 5. 计算幅度特性并进行打印。 五、实验过程原始记录（数据、图表、计算等） (一 ) 离散信号 % 14( ) ( )x n R n n=0:1:10; xn=ones(1,4),zeros(1,7); %输入时域序列向量 xn=R4(n) Xk8=fft(xn,8); %计算 xn 的 8 点 DFT Xk16=fft(xn,16); %计算 xn 的 16
12、 点 DFT Xk32=fft(xn,32); %计算 xn 的 32 点 DFT k=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 点 D
13、FT 对应的采样点频率 (关于归一化值 ) 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(幅度
14、); 运行结果： %x2(n) n=0:50; xn=cos(pi/4*n); %输入时域序列向量 cos(pi/4*n) Xk4=fft(xn,4); %计算 xn 的 4 点 DFT Xk8=fft(xn,8); %计算 xn 的 8 点 DFT Xk16=fft(xn,16); %计算 xn 的 16 点 DFT k=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);
15、 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/1
16、6; %产生 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 点 DFT Xk16=fft(xn,16); %计算
17、 xn 的 16 点 DFT Xk32=fft(xn,32); %计算 xn 的 32 点 DFT k=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),.); %绘制
18、 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;
19、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 点 DFT k1=0:length(abs(Xk1)-1;wk1=2*k1/8; %产生 8 点 DFT 对应的采样点频率 (关于归一化值 ) Xk2=fft(xn,16); %计算 xn 的 16 点 DFT k2=0:length(abs(Xk2)-1;wk2=2*k2/16; %产生 16
20、 点 DFT 对应的采样点频率 (关于归一化值 ) Xk3=fft(xn,32); %计算 xn 的 32 点 DFT k3=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(幅度 )
21、; 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=2kHz f=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:
22、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*k
23、3/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); y
24、label(幅度 ); 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 对应的采样点频率
25、(关于归一化值 ) 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 点傅立叶变换 ); subplot(3,1,2); stem(k2,abs(Xk32),.); xlabel(omega/pi); ylabel(幅度 ); title(X_3 (n) 32 点傅立叶变换 ); subplot(3,1
26、,3); stem(k3,abs(Xk64),.); xlabel(omega/pi); ylabel(幅度 ); title(X_3 (n) 64 点傅立叶变换 ); 运行结果： 六、实验结果及分析 1.分析 DFT 的变换区间对频域分析的作用，并说明 DFT 的物理意义。 答 : 以 (一 ) 离散信号中的 14( ) ( )x n R n 为例。 变换区 间对频域分析的作用 :随着变换区间的增大 ,36 变换区间的频域明显比 8 区间变换的频域容易识别一个周期内部的情况。即变换区间越长 ， 信号的基频越小。 物理意义 :假设 x(n)是 N 点 DFT 是 x(n)的 z 变换在单位圆上
27、的 N 点等间隔采样 ， 则 x(k)为 x(n)的傅里叶变换 x(e(jw)在区间 0,2上的 N 点等间隔采样 。 2.对周期信号的变换区间应该如何选取。如果周期信号的周期预先不知道，如何用 DFT 分析它的频谱。 答 ： 1） 变换区间的选择：对于非周期信号，假设频谱分辨率 F，而频谱分辨率直接和 FFT 的变换区间有关，因为 FFT 能够实现的频率分辨率是 2 /N 因此有最小的 N2 /F，根据此式可以选择 FFT 的变换区间；对于周期信号，周期信号的频谱是离散谱，只有用整数倍周期的长度作 FFT，得到的离散谱才能代表周期信号的频谱。 2） 周期信号的周期预先不知道时 ,可先截取 M 点进行 DFT，再将截取长度扩大 1 倍截取，比较结果，如果二者的差别满足分析误差要求，则可以近似表示该信号的频谱，如果不满足误差要求就继续将截取长度加倍，重复比较，直到结果满足要求。
