数字信号处理实验指导书96745

PAGEPAGE30数字信号处理实验指导书编者：韩萍何炜琨2004年8月目录TOC\o"1-3"\h\z实验一用FFT进行信号谱分析 1实验二频谱细化、频谱泄漏及栅栏效应 4实验三IIR数字滤波器的设计及滤波实现 7实验四FIR数字滤波器的设计及滤波实现 10实验五圆周卷积在信号处理中的应用 14附录 19实验一用FFT进行信号谱分析1、实验目的：进一步加深DFT算法原理和基本性质的理解（因为FFT只是DFT的一种快速算法，所以FFT的运算结果必然使DFT的基本性质）熟悉FFT算法原理和FFT子程序的应用。学习用FFT对连续信号和时域离散信号进行谱分析的算法，了解可能出现的分析误差及其原因，以便在实际中正确应用FFT。2、实验原理DFT是离散付立叶级数（DFS）在时域和频域内的两个主值序列构成的变换对：用DFT来分析有限长序列的频谱特性。FFT是DFT得快速算法，常用的两种快速算法有按时间抽取的FFT算法和按频率抽取的FFT算法，其中按时间抽取的FFT算法如下式：按频率抽取的FFT算法如下：，注意：DFT隐含周期性，、分别可以看作是一个周期为N的有限长序列、的一个周期。FFT并不是一种新的变换。3、实验步骤（1）复习DFT的定义、性质和用DFT作谱分析的有关内容。（2）复习FFT算法原理并读懂FFTMatlab语言的编程思想。（3）编制信号产生子程序，产生一下典型信号供谱分析使用。(4)编写主程序如图所示流程图，供参考。(5)按实验上机内容要求，上机实验并写出实验报告。实验一主程序框图4、上机实验内容（1）、对2中所给出的信号分别进行谱分析，试验参数如下：、、、、进行DFT时的点数N=8对于需进行采样，采样频率（2）令，用FFT计算8点的离散付立叶变换：，并根据DFT的对称性,由求出和，并与（1）中所得的结果进行比较。（3）令，重复（2）。5、思考题N=8时，和的幅频特性是否相同？为什么？N=16呢？如果周期信号的周期预先不知道，如何用FFT进行谱分析？6、实验报告要求简述实验原理及目的。结合实验中得到的给定典型序列幅频特性曲线，与理论结果相比较，并分析说明误差产生的原因以及用FFT进行谱分析时有关参数的选择方法。总结实验所得的主要结论。简要回答思考题。实验二频谱细化、频谱泄漏及栅栏效应1、实验目的在理论学习的基础上，通过本实验加深学生对DFT和FFT变换中频谱细化、频谱泄漏、栅栏效应误差的理解，使学生对用FFT进行信号频谱分析过程中存在的问题有一个明确的认识。2、实验原理频谱细化：在对长度为N的序列进行谱分析时，可以进行N点的DFT(FFT),当序列的点数比较少时,谱线也很稀疏,不容易看清谱的分布特点,若想将序列的频谱看得更加细致,一种方法就是通过在序列末尾补零,再进行DFT(FFT)就可以使得到的谱更加细致。栅栏效应：对采样信号的频谱，为提高计算效率，通常采用FFT算法进行计算，设数据点数为：

N=T/dt=T.fs

则计算得到的离散频率点为:

Xs(fi)，fi=i.fs/N,i=0，1，2，，N/2

这就相当于透过栅栏观赏风景，只能看到频谱的一部分，而其它频率点看不见，因此很可能使一部分有用的频率成份被漏掉，此种现象被称为栅栏效应。频谱泄漏:工程实际中能观测到的信号不可能是无限长的，只能从某时刻开始观测有限时间长度Ｔ的一段，这就相当于用一个窗函数对信号进行截断:

xs(t)=x(t)u(t)

截断后信号的频谱为：

Xs(f)=X(f)*U(f)

若x(t)为一单一频率f0的正弦函数，u(t)为矩形窗函数，按卷积分计算得到信号频谱为：

Xs(f)=[U(f+f0)+U(f-f0)]/2

信号频谱由原来的线谱，变为连续谱，能量泄漏到整个频带，这种现象称之为能量泄漏。3、实验步骤复习频谱细化方法，栅栏效应、频谱泄漏产生的原因。对实验一中的、、分别进行谱的细化。用505Hz正弦波信号的频谱分析来说明栅栏效应所造成的频谱计算误差。用DFT频谱展宽技术对300Hz正弦波信号分别用矩形窗截断和汉宁窗截断，观察其频谱泄漏情况。4、上机实验内容在实验一的基础上,频谱细化时选取的参数:N=16,32,的取样频率同实验一。编制相应的程序。（2）用505Hz正弦波信号的频谱分析来说明栅栏效应所造成的频谱计算误差,参数如下：

设定采样频率:fs=5120Hz，FFT计算点数为512，其离散频率点为：

fi=i.fs/N=i.5120/512=10，i=0，1，2，….，N/2

位于505Hz位置的真实谱峰被挡住看不见，看见的只是它们在相邻频率500Hz或510Hz处能量泄漏的值。

若设fs=2560Hz，则频率间隔df=5Hz，重复上述分析步骤，这时在505位置有谱线，我们就能得到它们的精确值。

从时域看，这个条件相当于对信号进行整周期采样，实际中常用此方法来提高周期信号的频谱分析精度。

编制相应的程序，并绘制频谱图。（3）对300Hz正弦波信号分别用矩形窗截断和汉宁窗截断,计算截断前后该正弦波信号的的频谱,观察其频谱泄漏情况。采样频率fs=1200Hz，窗的长度N=512。编制相应的程序，并绘制加窗前后的频谱图。5、思考题信号整周期采样时，其FFT谱中是否还存在能量泄漏误差?信号整周期采样时，其FFT谱中是否还存在栅栏效应误差?对信号的频谱进行细化时,对序列进行末尾补零点后,序列的长度必须是原序列长度的整数倍吗?6、实验报告要求简述实验原理及目的。(2)总结实验所得结果,给出克服栅栏效应和频谱泄漏的方法。(3)简要回答思考题。实验三IIR数字滤波器的设计及滤波实现1、实验目的熟悉用双线性变换法设计IIR数字滤波器的原理与方法。掌握数字滤波器的设计方法。通过观察对实际心电图信号的滤波过程，获得数字滤波的感性认识。2、实验原理基于模拟滤波器数学模型--微分方程的积分，利用对积分的数值逼近可以得到相应的数字滤波器。假设模拟滤波器的系统函数对模拟滤波器取样后逼近后，得到的数字滤波器系统函数可见，数字滤波器可通过下述变换关系有模拟滤波器得到上述s平面到z平面之间的变换是双线性变换关系，故称之为“双线性”变换。它对高阶系统同样成立。3、实验内容用双线性变换法设计一个巴特沃斯低通IIR数字滤波器。设计参数为:在低通频带内频率低于时,最大衰减小于1dB；在阻带内[，]频率区间上，最小衰减大于15dB。以为采样间隔，打印出数字滤波器在频率区间[0，]上的幅频响应特性曲线。用所设计的滤波器对实际心电图信号采样序列进行仿真滤波处理，并分别打印出滤波前后的心电图信号波形图，观察总结滤波作用与效果。4、实验步骤复习有关巴特沃斯模拟滤波器的设计和用双线性变换法设计IIR数字滤波器的内容，即变换关系。根据教材中p194页例题中给定的设计指标，调用MATLAB信号处理工具箱函数buttord和butter设计出相应的数字滤波器H(Z).编写滤波器仿真程序，完成对实际采集的心电图信号序列（具体数据见下面）的总响应序列，可直接调用MATLABfilter函数实现仿真。[附：人体心电图采样信号在测量过程中往往受到工业高频干扰，所以，必须经过低通滤波处理后，才能作为判断心脏功能的有用信息。下面的序列就是一个实际心电图信号采样序列样本，其中存在高频干扰，实验时，将其作为输入信号，滤除其中的干扰成分。在计算机上运行仿真程序，并调用绘图子程序，完成实验内容（2）和（3）。5、思考题用双线性变换法设计数字滤波器过程中,变换公式中T的取值对设计结果有无影响？为什么？6、实验报告要求简述实验目的和实验原理。由所打印的特性曲线及设计过程简述双线性变换法的特点。对比滤波前后的心电图信号波形,说明数字滤波器的滤波过程与滤波作用。简要回答思考问题。实验四FIR数字滤波器的设计及滤波实现1、实验目的掌握用窗函数法设计FIR数字滤波器的原理和方法。熟悉线性相位FIR数字滤波器特性。了解各种窗函数对滤波特性的影响。2、实验原理如果所希望的滤波器的理想频率响应函数为，则其对应的单位脉冲响应为窗函数设计方法的基本原理是用有限长单位脉冲响应序列逼近。由于往往是无限长序列，且是因果的，所以用窗函数将截断，并进行加权处理，得到就作为实际设计的FIR数字滤波器的单位脉冲响应序列，其频率响应函数为式中，N为所选窗函数的长度。这里需要注意的是：设计出的滤波器性能取决于窗函数的类型和窗口的长度N的取值，设计过程中，要根据对阻带最小衰减和过渡带宽度的要求选择合适的窗函数类型和窗口长度（各类窗函数的频率特性见教材的第三章）。选定窗函数和长度N后，求出单位脉冲响应和。是否满足要求，要进行演算。一般在尾部加零使长度满足2的整次幂，以便用FFT计算，如果观察细节，补领点数增多即可。如果不满足设计要求，则要重新选择窗函数类型和长度N,再次验算，直至满足要求。如果要求线性相位特性，则还必须满足根据上式中的正、负和长度N的奇、偶性又将线性相位FIR滤波器分为四类，要根据所设计的滤波特性正确选择其中一类。例如要设计线性相位低通特性，可选择一类，而不能选一类。3、实验步骤复习用窗函数法设计FIR数字滤波器一节内容，阅读本实验原理，掌握设计步骤。编写程序编写能产生矩形窗、升余弦窗的窗函数子程序。编写主程序。主程序框图如下（仅供参考），其中幅度特性要求用dB表示。设画图时，用打印幅度特性。第点对应的频率。为使曲线包络更接近的幅频特性，DFT变换的区间要选得大一些。例如窗口长度N=33时，可通过在末尾补零的方法，使长度变为64，再进行64点的DFT，则可得到更精确的幅度衰减特性曲线。4、上机实验内容用升余弦窗设计一线性相位低通FIR数字滤波器，截止频率，窗口长度N=15,33。要求在两种窗口长度情况下，分别求出，打印出相应的幅频特性和相频特性曲线，观察3dB带宽和20dB带宽。总结窗口长度N对滤波特性的影响。设计低通FIR数字滤波器时，一般以理想低通滤波特性为逼近函数，即其中，（2）N=33,，用两种窗函数设计线性相位低通滤波器。绘制相应的幅频特性曲线，观察3dB和20dB带宽以及阻带最小衰减，比较两种窗函数的滤波特性的影响。5、思考题如果给定通带截止频率和阻带截止频率以及阻带最小衰减，如何用窗函数法设计线性相位低通滤波器？写出步骤.如果要求用窗函数法设计带通滤波器，且给定上、下边带截止频率和，试求理想带通的单位脉冲响应6、实验报告要求简述实验目的和原理按照实验步骤及要求，比较各种情况下的滤波性能，说明窗口长度N和窗函数类型对滤波特性的影响。总结用窗函数设计法设计FIR滤波器的主要特点。简要回答思考题。 附录实验一用FFT进行信号谱分析1.实验一程序清单（MATLAB语言）clearallclcb=input('请选择信号：\n1:x1(n),2:x2(n),3:x3(n),4:x4(n),5:x5(n),6:x6(n),7:x7=x4+x5,8:x8=x4+jx5:\n','s');b=str2num(b);i=0;closeall;ifb==6N=input('请选择FFT变换区间长度N:16or32or64:\n','s');N=str2num(N);fs=64;n=0:N-1;x=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs);elseN=input('请选择FFT变换区间长度N:8or16:\n','s')N=str2num(N);endifb==1x=[11110000];elseifb==2x=[12344321];elseifb==3x=[43211234];elseifb==4n=0:N-1;x=cos(0.25*pi*n);elseifb==5n=0:N-1;x=sin((pi*n)/8);elseifb==7n=0:N-1;x=cos(n*pi/4)+sin(n*pi/8);elseifb==8n=0:N-1;x=cos(n*pi/4)+j*sin(n*pi/8);endendendendendendendf=fft(x,N);i=i+1;figure(i)printf(x,abs(f),abs(N),abs(b));ifN==16ifb==7k=conj(f);x4=(f+k)/2;figure(i+2)subplot(2,2,1);stem(abs(x4),'.')xlabel('k');ylabel('|Re[X7(k)]|');title('恢复后的|X4(k)|');x5=(f-k)/2;subplot(2,2,3)stem(abs(x5),'.');xlabel('k')ylabel('|jIm[X7(k)]|');title('恢复后的|X5(k)|');endifb==8k(1)=conj(f(1));form=2:Nk(m)=conj(f(N-m+2));endfe=(f+k)/2;fo=(f-k)/(2);xr=ifft(fe,N);b=4;figure(i+1)printf(xr,abs(fe),abs(N),abs(b));xi=ifft(fo,N)/j;b=5;figure(i+2)printf(xi,abs(fo),abs(N),abs(b));endendb=input('请选择信号：\n1:x1(n),2:x2(n),3:x3(n),4:x4(n),5:x5(n),6:x6(n),7:x7=x4+x5,8:x8=x4+jx5,0(退出):\n','s');b=str2num(b);end子函数：通用画图程序functionprintf(j,k,N,b)subplot(2,2,1);if(b~=8)ifb==1|b==2|b==3M=3.2;n=0:7;stem(n,j,'.');n=0:7;m=zeros(8);elseM=N*0.4;n=0:N-1stem(n,j,'.');n=0:N-1;m=zeros(N);endholdon;plot(n,m);t=max(j);xlabel('n');string=['x',num2str(b),'(n)的波形'];ifb==4ylabel('x4(n)');elseifb==5ylabel('x5(n)');endtext(M,(t*0.8),string);endsubplot(2,2,3);n=0:N-1;stem(n,k,'.');t=max(k);xlabel('K');string=['x',num2str(b),'(n)的N=',num2str(N),'点FFT'];text((N*0.4),(t*0.8),string);ifb==4ylabel('|X8e(k)|');elseifb==5ylabel('|X8o(k)|');end2.程序运行结果图1-1x1(n)及其8点和16点DFT图1-2x2(n)及其8点和16点DFT图1-3x3(n)及其8点和16点DFT图1-4x4(n)的8点和16点波形及其DFT图1-5x5(n)8点和16点波形及其DFT图1-6x6(n)16点、32点和64点采样序列波形及其DFT图1-7x7(n)＝x4(n)+x5(n)的8点和16点波形及其DFT图1-8|X4(k)|=|Re[X7(k)]|及|X5k)|=|jIm[X7(k)]|图1-9x7(n)的8点和16点DFT图1-10|X8e(k)|及IDFT[X8e(k)]图1-11|X8O(k)|及IDFT[X8O(k)]实验三IIR数字滤波器的设计及滤波实现1.实验三程序清单clearallclcx=[-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]k=1;closeallfigure(1)subplot(2,2,1)n=0:55;stem(n,x,'.');axis([056-10050]);holdonn=0:60;m=zeros(61);plot(n,m)xlabel('n')ylabel('x(n)');title('心电图信号采样序列x(n)');B=[0.090362*0.090360.09036];A=[1.2686-0.7051];A1=[1.0106-0.3583];A2=[0.9044-0.2155];while(k<=3)y=filter(B,A,x);x=yifk==2A=A1;endifk==3A=A2;endk=k+1;endsubplot(2,2,3);n=0:55;stem(n,y,'.');axis([056-155]);holdonn=0:60;m=zeros(61);plot(n,m);xlabel('n');ylabel('y(n)');title('三级滤波后的心电图信号');A=[0.090360.18720.09036]B1=[1,-1.2686,0.7051];B2=[1,-1.0106,0.3583];B3=[1,-0.9044,0.2155];[H1,w]=freqz(A,B1,100);[H2,w]=freqz(A,B2,100);[H3,w]=freqz(A,B3,100);H4=H1.*(H2);H=H4.*(H3);mag=abs(H);db=20*log10((mag+eps)/max(mag));subplot(2,2,2)plot(w/pi,db);axis([0,0.5,-50,10]);title('滤波器的幅频响应曲线');xlabel('w/n');ylabel('20lg|H(w)|')2.程序运行结果图3-1心电图信号波形图实验四FIR数字滤波器的设计及滤波实现1.实验四程序清单clearallclcN=input('InputThelengthofthewindowsN:[10,64]\n');w=input('INputTheWc:[0.5,pi]\n');b=1;closeall;i=0;n=[0:(N-1)];hd=ideal(w,N);k=input('请选择窗口类型：1（boxcar）2(hamming)3(hanning)4(blackman):\n','s');k=str2num(k);ifk==1B=boxcar(N);string=['Boxcar','N=',num2str(N)];elseifk==2B=hamming(N);string=['Hamming','N=',num2str(N)];elseifk==3B=hanning(N);string=['Hanning','N=',num2str(N)];elseifk==4B=blackman(N);string=['Blackmanr','N=',num2str(N)];endendh=hd.*(B)';[H,m]=freqz(h,[1],1024,'