实验一-FFT与谱分析(共34页)

1、精选优质文档-倾情为你奉上数字信号处理实验报告姓名：朱道萍班级：通信1401学号：9专心-专注-专业实验一 FFT与谱分析1.实验目的(1) 加深对DFT算法原理和物理意义的理解。(2)熟悉FFT算法，增强对FFT结果的分析能力。(3)掌握用FFT对连续时间信号和离散时间信号进行谱分析的方法，了解误差及其产生的原因。2.实验原理(1)复习 DFT 的定义、性质和用 DFT 作谱分析的有关内容。(2)复习 FFT 算法原理与编程思想。(3)查询Matlab中fft函数的用法，学会通过改变fft的参数计算不同DFT。(4)编制信号产生子程序，产生一下信号：(5) 编写实验程序： 参考图1-1所示的

2、程序流程图，编写实验程序。3.实验内容 对 2 中所给出的11个信号逐一进行谱分析。各信号的FFT变换区间N、连续时间信号x6(t)的采样频率fs等参数的参考取值如下： a) x1n,x2n, x3n, x4n, x5n , x7n, x8n: N=8, 16 b) x6(t): fs=64Hz，N=16, 32, 64 c) x9n, x10n, x11n: N=16 x1n:N=8n=0:N-1;x1=1 1 1 1 0 0 0 0;m=0:2*N-1;x2=1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0;f1=fft(x1,N);f2=fft(x2,2*N);subplo

3、t(2,2,1)stem(n,x1)subplot(2,2,2)stem(n,abs(f1)subplot(2,2,3)stem(m,x2)subplot(2,2,4)stem(m,abs(f2)x2n:N=8 n=0:N-1; x1=1 2 3 4 4 3 2 1;m=0:2*N-1; x2=1 2 3 4 4 3 2 1 0 0 0 0 0 0 0 0;f1=fft(x1,N);f2=fft(x2,2*N);subplot(2,2,1)stem(n,x1) subplot(2,2,2)stem(n,abs(f1)subplot(2,2,3)stem(m,x2)subplot(2,2,4)s

4、tem(m,abs(f2)x3n:N=8n=0:N-1; x1=4 3 2 1 1 2 3 4;m=0:2*N-1; x2=4 3 2 1 1 2 3 4 0 0 0 0 0 0 0 0;f1=fft(x1,N);f2=fft(x2,2*N)subplot(2,2,1)stem(n,x1)subplot(2,2,2)stem(n,abs(f1)subplot(2,2,3)stem(m,x2) subplot(2,2,4)stem(m,abs(f2)x4n:N=8;n=0:N-1; x1=cos(0.25*pi*n);m=0:2*N-1; x2=cos(0.25*pi*m); f1=fft(x1

5、,N);f2=fft(x2,2*N); subplot(2,2,1)stem(n,x1) subplot(2,2,2)stem(n,abs(f1) subplot(2,2,3)stem(m,x2) subplot(2,2,4)stem(m,abs(f2)x5n:N=8;n=0:N-1; x1=sin(0.125*pi*n);m=0:2*N-1; x2=sin(0.125*pi*m);f1=fft(x1,N);f2=fft(x2,2*N);subplot(2,2,1)stem(n,x1) subplot(2,2,2)stem(n,abs(f1)subplot(2,2,3)stem(m,x2)su

6、bplot(2,2,4)stem(m,abs(f2)x6n:N=16;fs=64;n=0:N-1; x1=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs);m=0:2*N-1; x2=cos(8*pi*m/fs)+cos(16*pi*m/fs)+cos(20*pi*m/fs);l=0:4*N-1; x3=cos(8*pi*l/fs)+cos(16*pi*l/fs)+cos(20*pi*l/fs);f1=fft(x1,N);f2=fft(x2,2*N);f3=fft(x3,4*N);subplot(3,2,1)stem(n,x1)subplot(3,2

7、,2)stem(n,abs(f1)subplot(3,2,3)stem(m,x2)subplot(3,2,4)stem(m,abs(f2)subplot(3,2,5)stem(l,x3) subplot(3,2,6)stem(l,abs(f3)x7n,x8n:N=8; n=0:N-1; x51=sin(0.125*pi*n);x41=cos(0.25*pi*n);x71=x41+x51;m=0:2*N-1; x52=sin(0.125*pi*m);x42=cos(0.25*pi*m);x72=x42+x52;f41=fft(x41,N);f51=fft(x51,N);f71=fft(x71,N

8、);f42=fft(x42,2*N);f52=fft(x52,2*N);f72=fft(x72,2*N);figure(1)subplot(2,2,1)stem(n,real(f71)subplot(2,2,2)stem(n,f41)subplot(2,2,3)stem(n,imag(f71)subplot(2,2,4)stem(n,imag(f51) figure(2)subplot(2,2,1)stem(m,real(f72)subplot(2,2,2)stem(m,f42)subplot(2,2,3)stem(m,imag(f72)subplot(2,2,4)stem(m,imag(f5

9、2) x9n,x10n:N=8;n=0:N-1; x51=sin(0.125*pi*n);x41=cos(0.25*pi*n);x81=x41+j*x51;m=0:2*N-1; x52=sin(0.125*pi*m);x42=cos(0.25*pi*m);x82=x42+j*x52;f41=fft(x41,N);f51=fft(x51,N);f81=fft(x81,N);f83=f81;for i=2:N f83(i)=conj(f81(N+2-i);end f42=fft(x42,2*N);f52=fft(x52,2*N);f82=fft(x82,2*N);f84=f82;for i=2:2

10、*N f84(i)=conj(f82(2*N+2-i);end figure(1)subplot(2,2,1)stem(n,(f81+f83)/2)subplot(2,2,2)stem(n,f41)subplot(2,2,3)stem(n,abs(f81-f83)*1i/2)subplot(2,2,4)stem(n,abs(f51) figure(2)subplot(2,2,1)stem(m,(f82+f84)/2)subplot(2,2,2)stem(m,f42)subplot(2,2,3)stem(m,abs(f82-f84)*1i/2)subplot(2,2,4)stem(m,abs(f

11、52)x11n:N=16;n=0:N-1; x1=1 2 3 4 0 0 0 0 0 0 0 0 0 0 0 0;x2=1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4;x3=1 0 0 0 2 0 0 0 3 0 0 0 4 0 0 0;f1=fft(x1,N);f2=fft(x2,N);f3=fft(x3,N); subplot(3,2,1)stem(n,x1)subplot(3,2,2)stem(n,abs(f1)subplot(3,2,3)stem(n,x2)subplot(3,2,4)stem(n,abs(f2)subplot(3,2,5)stem(n,x3)subpl

12、ot(3,2,6)stem(n,abs(f3)4. 思考题(1)在N=8时，x2n和x3n的幅度谱会相同吗？为什么？N=16时，x2n和x3n的幅度谱会 相同吗？为什么？答：不相同，由DFT的定义可知,xn的作用是对相同的K值的幅度值，所以虽然不影响频谱分，布但是影响频率幅值.由于x2n和x3n的数值分布范围不同,所以导致在相同K值时幅值不同，因此x2n和x3n的幅频特性不相同。N=16时也不相同。(2) 比较x4n的8点和16点的DFT波形，这说明什么？再比较x5n的8点和16点的DFT 波形，这说明了什么？分析x4n和x5n的这些频谱异同产生的原因。答：x4n的8点和16点波形相同，x5n

13、的8点和16点波形相同，说明仅当函数为偶函数 时，8点和16点DFT波形相同。(3) 比较x6(t)的64Hz采样时16，32，64点DFT波形，分析它们的异同产生的原因。答：波形相同，因为是对同一个数字频谱的采样。(4) 比较x6(t)的16Hz采样和64Hz采样时64点DFT波形，这说明什么？答：64Hz是16Hz扩张四倍后的波形，说明频率改变频谱的压缩或扩张。(5)找出X7k与X4k=DFTx4n、X5k=DFTx5n的关系。答：X7k的实部和X4k相等，X7k的虚部和X5k的虚部相等。(6)找出X8k与X4k=DFTx4n、X5k=DFTx5n的关系。答：X4k=(X8k+ X8*n-

14、k)/2, X5k=j/2(X8k- X8*n-k)。(7) 比较x9n、x10n和x11n的DFT结果，分析信号末尾补零、周期性延拓和时域插值在频 谱上的变化。答：信号末尾补零频谱不发生变化，周期性延拓频谱变化是周期性的，时域插值频谱压缩。(8)如果周期信号的周期预先不知道，如何用FFT进行谱分析？答：周期信号的周期预先不知道时,可先截取M点进行DFT,再将截取长度扩大1倍截取,比较 结果,如果二者的差别满足分析误差要求,则可以近似表示该信号的频谱,如果不满足误差 要求就继续将截取长度加倍,重复比较,直到结果满足要求。5.分析与结论 FFT变换即快速傅里叶变换的性质同DFT即离散傅里叶变换相

15、同。离散傅里叶变换有两个物理意义，一是，是对该序列的傅里叶变换w的抽样或者说对Z变换单位圆内的抽样。二是，将该序列进行周期延拓后的傅里叶级数变换的主值序列。用FFT作谱分析时，若给出的是连续信号，需要根据其最高频率确定采样速率，根据频率分辨率选取采样点数N，对其进行软件采样，产生对应序列。若给出的是周期序列，则应该截取周期的整数倍进行谱分析，这样可以尽量减少实验造成的误差。实验二 IIR数字滤波器设计1.实验目的(1)熟悉用双线性变换法设计IIR数字滤波器的原理与方法。(2)掌握数字滤波器的计算机仿真方法。(3)通过观察对实际心电图信号的滤波作用，获得数字滤波的感性知识。2.实验原理(1)II

16、R数字滤波器的系统函数可以写成 设计经典的IIR数字滤波器，就是要确定该系统函数中的2组系数： Num: b0,b1,.,bM Den: a1,a2,.,aN 也就是说，这两组系数确定了，IIR数字滤波器的理论设计就完成了。接下去可以进行 结构变换和系数量化。(2) IIR 数字滤波器的设计可以借助模拟滤波器的设计方法。先设计一个模拟滤波器，然后转 换成数字滤波器。(3)模拟滤波器的原型有很多种。常用的低通模拟滤波器原型有四种：巴特沃斯（Butterworth） 滤波器、切比雪夫 I型（Chebyshev ）滤波器、切比雪夫II型（Chebyshev ）滤波器、 椭圆（Elliptic）滤波器

17、。这些滤波器的增益（Gain）曲线如图 2-1 所示。(4) 滤波器设计指标可分为模拟指标和数字指标，绝对指标和相对指标。一般来说，具体指 标有通带的边界频率、阻带的边界频率，通带最大衰减，阻带最小衰减四种。低通滤波器 指标具体含义如图2-2所示。 (5)根据指标设计模拟滤波器，得到传输函数Hc(s)。用双线性变换法， 其中，为了防止频响曲线发生畸变，设计指标的模拟频率与数字频率间的关系需做预畸变：3.实验内容(1) 复习有关巴特沃斯、切比雪夫 I 型、切比雪夫 II 型、椭圆四种模拟滤波器设计和用双线 性变换法设计 IIR 数字滤波器的内容。(2) 人体心电图信号在测量过程中往往受到工业高频

18、（主要是50Hz市电及其高次谐波）干扰， 所以必须经过低通滤波处理后，才能作为判断心脏功能的有用信息。 xn为一实际心电图信号采样序列样本: xn=-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,-38,-60,-84,-90,-66,-32,-4,-2,-4,8,12, 12,10,6,6,6,4,0,0,0,0,0, -2,-4,0,0,0 -2 -2,0,0,-2,-2,-2,-2,0， 画出该信号的时域波形和频谱。从频谱图中可以看到该信号存在高频干扰，噪声主要 集中在大于 0.36rad 的频段。(3) 用双线性变换法设

19、计低通 IIR 数字滤波器，滤除信号xn中的干扰成分。根据待处理信号 的频谱，确定滤波器的设计指标。例如： 通带0,0.2内，最大衰减小于 1dB, 阻带0.3,内，最小衰减大于 15dB。(4) 利用Matlab函数设计IIR滤波器 a. 先采用巴特沃兹原型设计。根据设计指标，调用Matlab信号处理工具箱函数buttord() 和butter()，得到以巴特沃兹模拟低通为原型的数字滤波器H(z)。以0.02为采样间隔， 画出数字滤波器在频率区间0,)上的频率响应特性曲线。 b. 编写滤波器仿真程序，调用filter()函数计算H(z)对心电图信号采样序列xn的响应yn。 画出滤波后心电图信

20、号的时域波形和频谱。 c. 调用Matlab的cheb1ord()，cheby1()，cheb2ord()，cheby2()，ellipord()，ellip()等函数，分别得到切比雪夫I型、切比雪夫II型、椭圆滤波器为模拟原型的数字滤波器。重复步骤(4)a，(4)b。(5)利用Matlab的滤波器设计分析工具FDATool（Filter Design & Analysis Tool）设计IIR滤波器。 a. 在Matlab的“Commad Window”里输入fdatool，打开滤波器设计分析工具。 b. 在设计窗口中输入设计指标，如图2-3所示： 其中“Match exactly”这项，如

21、果选择“stopband”，则得到的滤波器阻带指标正好匹配，通带留有裕量；如果选择“passband”，则得到的滤波器通带指标正好匹配，阻带留有裕量。c. 点击 Design Filter ，就可以得到IIR滤波器。设计窗口中会显示设计得到的滤波器的频响 曲线。d. 在“Current Filter Information”一栏里，可以看到滤波器的结构信息，IIR滤波器可以有很 多种结构，选择不同的结构，得到的滤波器系数形式也各不相同。滤波器的结构可以通过 菜单“Edit”转换。如同一个6阶直接II型转置型IIR滤波器，图2-4是单一节点结构（Single Section），设计结果为2组系数

22、，分子分母各7个系数；而图2-5是二阶节点级联结构（Second-Order Sections，SOS），设计结果为3个二阶节点，每个节点分子分母各3个系数。 e. 点击菜单“File”-“Export”，可以把设计得到的滤波器系数导出到Matlab的工作空间（Workspace）或者存为文件。f. 使用设计得到的滤波器对心电图信号做滤波。可参考步骤(4)b。g. 在“Design Method”的下拉菜单中选用不同的模拟滤波器原型，分别设计相应的IIR数字滤波器。分别画出这些数字滤波器在频率区间0, )上的频率响应特性曲线。滤波前：Wp=0.2*pi; Ws=0.3*pi; Rp=1; Rs

23、=15; x=-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,-38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,4,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0;n=0:length(x)-1; N,Wn=buttord(Wp/pi,Ws/pi,Rp,Rs); b,a=butter(N,Wn); y=filter(b,a,x); figure(1);freqz(b,a,50); figure(2)subplot(2,2,1) stem

24、(n,x,.);axis(0 56 -100 50);hold on;xlabel(it n);ylabel(it xn);title(Cardiogram sequence);subplot(2,2,2) stem(n,y,.);axis(0 56 -100 50);xlabel(it n);ylabel(it yn);title(Cardiogram sequence filtered);m=0:pi/1024:pi/1024*1023;subplot(2,2,3) plot(m,abs(fft(x,1024);axis(0 511 -50 500);xlabel(itomega);yla

25、bel(it X(ejomega);title(Spectrum of cardiogram);subplot(2,2,4) plot(m,abs(fft(y,1024);axis(0 511 -50 500);xlabel(it omega);ylabel(it Y(ejomega);title(Spectrum of cardiogram filtered);巴特沃兹滤波器 Wp=0.2*pi;Ws=0.3*pi;Rp=1;Rs=15;x=-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,-38,-60,-84,-90,-66

26、,-32,-4,-2,-4,8,12,12,10,6,6,6,4,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0;n=0:55;N,Wn=cheb1ord(Wp/pi,Ws/pi,Rp,Rs);b,a=cheby1(N,Rp,Wn);y=filter(b,a,x); figure(1);freqz(b,a,50);figure(2)subplot(2,2,1)stem(n,x,.);axis(0 56 -100 50);hold on;xlabel(n);ylabel(xn);title(Cardiogram sequence);subplot(2,

27、2,2)stem(n,y,.);axis(0 56 -100 50);xlabel(n);ylabel(yn);title(Cardiogram sequence filtered); m=0:1023; subplot(2,2,3)plot(m,abs(fft(x,1024);axis(0 511 -50 500);xlabel(w);ylabel(X(ejw);title(Spectrum of cardiogram);subplot(2,2,4)plot(m,abs(fft(y,1024);axis(0 511 -50 500);xlabel(w);ylabel(Y(ejw);title

28、(Spectrum of cardiogram filtered)切比雪夫滤波器切比雪夫滤波器Wp=0.2*pi;Ws=0.3*pi;Rp=1;Rs=15; x=-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,-38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,4,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0;n=0:55;N,Wn=cheb2ord(Wp/pi,Ws/pi,Rp,Rs);b,a=cheby2(N,Rs,Wn);y=fi

29、lter(b,a,x); figure(1);freqz(b,a,50);figure(2)subplot(2,2,1)stem(n,x,.);axis(0 56 -100 50);hold on;xlabel(n);ylabel(xn);title(Cardiogram sequence); subplot(2,2,2)stem(n,y,.);axis(0 56 -100 50);xlabel(n);ylabel(yn);title(Cardiogram sequence filtered); m=0:1023;subplot(2,2,3)plot(m,abs(fft(x,1024);axi

30、s(0 511 -50 500);xlabel(w);ylabel(X(ejw);title(Spectrum of cardiogram);subplot(2,2,4)plot(m,abs(fft(y,1024);axis(0 511 -50 500);xlabel(w);ylabel(Y(ejw);title(Spectrum of cardiogram filtered);椭圆滤波器滤波后的时域波形频谱和频率响应特性曲线Wp=0.2*pi;Ws=0.3*pi;Rp=1;Rs=15;x=-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0

31、,-16,-38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,4,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0;n=0:55;N,Wn=ellipord(Wp/pi,Ws/pi,Rp,Rs);b,a=ellip(N,Rp,Rs,Wn);y=filter(b,a,x);figure(1);freqz(b,a,50); figure(2)subplot(2,2,1)stem(n,x,.);axis(0 56 -100 50);hold on;xlabel(n);ylabel(xn);title(Cardi

32、ogram sequence);subplot(2,2,2)stem(n,y,.);axis(0 56 -100 50);xlabel(n);ylabel(yn);title(Cardiogram sequence filtered);m=0:1023;subplot(2,2,3)plot(m,abs(fft(x,1024);axis(0 511 -50 500);xlabel(w);ylabel(X(ejw);title(Spectrum of cardiogram);subplot(2,2,4)plot(m,abs(fft(y,1024);axis(0 511 -50 500);xlabe

33、l(w);ylabel(Y(ejw);title(Spectrum of cardiogram filtered);FDATool设计IIR滤波器4. 思考题(1) 用双线性变换法设计数字滤波器过程中，变换公式 中采样间隔 T 的取值，对 设计结果有无影响？为什么？答：无影响。依靠双线性变换是建立起来s平面和z平面的单值映射关系，因此可以有效避 免频谱混叠现象，无论T取何值都是单值映射关系，对设计结果不会有影响(2) 设计过程中，双线性变换法得到的数字滤波器的频响曲线与模拟滤波器的频响曲线是否 存在畸变？是什么原因造成的？答：存在畸变，因为连续时间频率变量和离散时间频率变量的变换必定是非线性的

34、，所以会 产生畸变。5.分析与结论 双线性变换法是将S平面压缩变换到S1平面，再由变换关系Z=exp(S1*T)将S1变换到Z平面上，这样S平面和Z平面之间就可以建立起一一对应的单值关系，从而消除多值变换性。上述滤波前后的心电图信号波形知，数字滤波器的滤波过程就是将输入序列通过一定的运算变换成输出序列，滤波的作用就是通过一定的计算消除或削弱一些不需要的频率分量的干扰。本次实验中用了双线性变换法、冲击响应不变法及直接设计法设计了IIR滤波器。通过本次实验，我更直观的看到了滤波器的特性，对理论有了更深刻的认识。实验三 FIR数字滤波器设计1.实验目的(1)掌握用窗函数法设计FIR 数字滤波器的原理

35、和方法。(2)熟悉线性相位FIR 数字滤波器特性。(3)了解窗函数类型及窗口长度对滤波特性的影响。2.实验原理顾名思义，FIR数字滤波器的单位冲激响应（Impulse Reponse）hn是有限长的，所以设计任务就是确定单位冲激响应hn的N个值（滤波器的阶数为M=N-1）。FIR数字滤波器处理信号的基本算法就是hn的N个值和输入信号做卷积运算。FIR数字滤波器的设计方法有很多种，比较常用的有窗函数法。窗函数法设计FIR滤波器步骤：(1)根据设计指标写出理想滤波器的频响函数Hd(ejw) 其中M为滤波器阶数。(2) 求出理想滤波器的脉冲响应(3) 根据阻带衰减确定窗函数形状： 阻带衰减 a. 如

36、果使用Blackman族的窗函数，查表3-1，选择最大逼近误差绝对值大雨阻带衰减A，同时近似过渡带宽最小的窗函数类型。 b. 如果使用Kaiser窗，计算参数(4) 根据过渡带宽计算窗函数长度 过渡带宽 =|s-p| 假设滤波器阶数为M，即窗函数长度N=M+1。 a. 如果使用Blackman族的窗函数，查表3-1，通过近似过渡带宽求出M的值。 b. 如果使用Kaiser窗，M=(A-8)/2.285，根据线性相位理论，设计LPF或者BPF，M取 整数即可。设计HPF或带BSF，M必须取偶数。(5) 给理想滤波器的脉冲响应加窗，得到实际滤波器的系数。Hn=hdnwn(6) 计算实际滤波器的频响

37、H(ejw)，验证是否达到设计指标。如果达到指标，涉及结束；如果没有达到，返回步骤（3），更改窗函数类型或者承诺长度，重新设计。3.实验内容(1)复习用窗函数法设计 FIR 数字滤波器的方法，阅读本实验原理，掌握设计步骤。(2)分别 Blackman Family 中的最合适的窗设计FIR 滤波器，要求： a.高通滤波器 b.滤波器指标 p=0.7, s=0.54, p=3dB, s=40dB 根据指标设计滤波器，绘制相应的幅频特性曲线，观察是否达到设计指标。如果没有达到 指标，则更改窗口类型或者窗口长度，重新设计。(3) 设计指标不变，换用 Kaiser 窗设计。绘制相应的幅频特性曲线，观察

38、是否达到设计指标。 如果没有达到指标，则更改窗口类型或者窗口长度，重新设计。(4) 将步骤(2)中窗口类型改为 Blackman Family 中的其他 4 种，窗口长度选用 N=41 (即M=40)，设计滤波器。Hanning Windowclose allclearwp=0.7*pi;ws=0.54*pi;rp=3;rs=40;trans_width=abs(wp-ws); wc=(wp+ws)/2/pi; M=ceil(6.2*pi/trans_width)/2)*2; r=max(-20*log(1-power(10,(-rp/20),rs);h=fir1(M,wc, high,hanning(M+1) H=fft(h,512); plot(0:511/256,20*log10(abs(H) axis(0,1,-125,5);title(Hanning Window)ylabe