《数字信号处理上机实验指导书(38页).docx》由会员分享，可在线阅读，更多相关《数字信号处理上机实验指导书(38页).docx（38页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、-数字信号处理上机实验指导书-第 38 页数字信号处理上机实验指导书电子与信息工程学院目录实验一离散信号产生及频谱的绘制4实验二时域采样与频域采样9实验三离散傅立叶变换及谱分析16实验四用FFT对信号作频谱分析21实验五双线性变换法设计IIR数字滤波器32实验六 FIR数字滤波器设计与软件实现38附录：实验用MATLAB语言工具箱函数简介45引言“数字信号处理”是一门理论和实验密切结合的课程，为了深入地掌握课程内容，应当在学习理论的同时，做习题和上机实验。上机实验不仅可以帮助学生深入地理解和消化基本理论，而且能锻炼初学者的独立解决问题的能力。所以，根据本课程的重点要求编写了四个实验。第一章是全
2、书的基础内容，抽样定理、时域离散系统的时域和频域分析以及系统对输入信号的响应是重要的基本内容。由于第一章大部分内容已经在前期信号与系统课程中学习完，所以可通过实验一帮助学生温习以上重要内容，加深学生对“数字信号处理是通过对输入信号的一种运算达到处理目的” 这一重要概念的理解。这样便可以使学生从信号与系统课程顺利的过渡到本课程的学习上来。第二章、三章DFT、FFT是数字信号处理的重要数学工具，它有广泛的使用内容。采用实验二、实验三加深理解DFT的基本概念、基本性质。FFT是它的快速算法，必须学会使用。数字滤波器的基本理论和设计方法是数字信号处理技术的重要内容。学习这一部分时，应重点掌握IIR和F
3、IR两种不同的数字滤波器的基本设计方法。IIR滤波器的单位冲激响应是无限长的，设计方法是先设计模拟滤波器，然后再通过SZ平面转换，求出相应的数字滤波器的系统函数。这里的平面转换有两种方法，即冲激响应不变法和双线性变换法，后者没有频率混叠的缺点，且转换简单，是一种普遍应用的方法。FIR滤波器的单位冲激响应是有限长的，设计滤波器的目的即是求出符合要求的单位冲激响应。窗函数法是一种基本的，也是一种重要的设计方法。学习完第七章后可以进行实验四。由于数字信号处理实验的主要目的是验证数字信号处理的有关理论，进一步理解巩固所学理论知识。所以，实现实验用的算法语言可以有许多种，但为了提高实验效率，要求学生用编
4、程效率比C语言高好几倍的MATLAB语言。下面介绍MATLAB的主要特点。（有关MATLAB的启动、程序运行和有关信号处理工具箱函数等内容将放到最后附录中介绍。）MATLAB是一种交互式的以矩阵为基本数据结构的系统。在生成矩阵对象时，不要求明确的维数说明。所谓交互式，是指MATLAB的草稿纸编程环境。即用户每输入一条命令并按回车键，MATLAB系统便解释执行之，并显示执行结果。根据该结果，用户立即知道刚输入的命令的正确性，或利用中间结果进行其他处理等。与C语言或FORTRON语言做科学数值计算的程序设计相比较，利用MATLAB可节省大量的编程时间。将其用于数字信号处理实验，则可大大提高实验效率
5、，在有限的上机时间内，实验内容可增加几倍。例如，C语言FFT子程序有70多行，而用MATLAB只调用一个fft函数即可实现对序列进行FFT计算。另外，MATLAB的工具箱及图形显示（打印）功能，可满足各层次人员直观、方便的进行分析、计算和设计工作，从而可大大节省时间。例如，序列的卷积、滤波，系统函数H(z)的幅频特性和相频特性等计算，均有现成的工具箱函数。而用其它算法语言完成这些计算的编程比较麻烦，且程序较长。由于上述特点，在美国一些大学里，MATLAB已成为辅助教学的有益工具。MATLAB已成功地用于数字信号处理课程中的问题分析、实验、滤波器设计及计算机模拟。附录中所介绍的信号处理工具箱函数
6、及绘图函数基本可满足本教材所要求的上机实验需要。对序列进行谱分析的MATLAB程序及运行结果见附录。实验一离散信号产生及频谱的绘制一、实验目的: （1）熟悉Matlab环境。（2）掌握 Matlab 中一些基本函数的建立方法（2）通过编程绘制的幅度相位谱加深理解系统的特性二、实验内容1、编写程序，产生以下离散序列：（1）f(n)=(n) (-3n4)n1=-3;n2=4;n0=0;n=n1:n2;x=n=n0;stem(n,x,filled);axis(n1,n2,0,1.1*max(x);xlabel(时间(n);ylabel(幅度x(n);title(单位脉冲序列);（2）f(n)=u(n
7、) (-5n=n0;stem(n,x,filled);axis(n1,n2,0,1.1*max(x);xlabel(时间(n);ylabel(幅度x(n);title(单位阶跃序列);box（3） (0n16)n1=16;a=-0.1;w=1.6*pi;n=0:n1;x=exp(a+j*w)*n);subplot(2,2,1);plot(n,real(x);title(复指数信号的实部);subplot(2,2,3);stem(n,real(x),filled);title(复指数序列的实部);subplot(2,2,2);plot(n,imag(x);title(复指数信号的虚部);subp
8、lot(2,2,4);stem(n,imag(x),filled);title(复指数序列的虚部);box %box on 加边框/ box off不加边框/ box 加或不加切换，注意只对上面一个图有效（4）f(n)=3sin(n/4) (0nM，比原序列尾部多N-M个零点；如果NM，z则=IDFT发生了时域混叠失真，而且的长度N也比x(n)的长度M短，因此。与x(n)不相同。在数字信号处理的应用中，只要涉及时域或者频域采样，都必须服从这两个采样理论的要点。 对比上面叙述的时域采样原理和频域采样原理，得到一个有用的结论，这两个采样理论具有对偶性：“时域采样频谱周期延拓，频域采样时域信号周期延
9、拓”。因此放在一起进行实验。二、实验内容及步骤（1）时域采样理论的验证。给定模拟信号，式中A=444.128，=50，=50rad/s，它的幅频特性曲线如图3.1图2.1 的幅频特性曲线现用DFT(FFT)求该模拟信号的幅频特性，以验证时域采样理论。按照的幅频特性曲线，选取三种采样频率，即=1kHz，300Hz，200Hz。观测时间选。为使用DFT，首先用下面公式产生时域离散信号，对三种采样频率，采样序列按顺序用，表示。因为采样频率不同，得到的，的长度不同， 长度（点数）用公式计算。选FFT的变换点数为M=64，序列长度不够64的尾部加零。X(k)=FFTx(n) ， k=0,1,2,3,-,
10、M-1式中k代表的频率为 。要求： 编写实验程序，计算、和的幅度特性，并绘图显示。观察分析频谱混叠失真。（2）频域采样理论的验证。给定信号如下：编写程序分别对频谱函数在区间上等间隔采样32和16点，得到：再分别对进行32点和16点IFFT，得到：分别画出、的幅度谱，并绘图显示x(n)、的波形，进行对比和分析，验证总结频域采样理论。提示：频域采样用以下方法容易变程序实现。 直接调用MATLAB函数fft计算就得到在的32点频率域采样 抽取的偶数点即可得到在的16点频率域采样，即。 当然也可以按照频域采样理论，先将信号x(n)以16为周期进行周期延拓，取其主值区（16点），再对其进行16点DFT(
11、FFT),得到的就是在的16点频率域采样。四、思考题： 如果序列x(n)的长度为M，希望得到其频谱在上的N点等间隔采样，当NM时， 如何用一次最少点数的DFT得到该频谱采样？五、实验报告及要求a）运行程序打印要求显示的图形，。b) 分析比较实验结果，简述由实验得到的主要结论c) 简要回答思考题d) 附上程序清单和有关曲线。六、实验程序清单1 时域采样理论的验证程序清单% 时域采样理论验证程序exp2a.mTp=64/1000;%观察时间Tp=64微秒%产生M长采样序列x(n)% Fs=1000;T=1/Fs;Fs=1000;T=1/Fs;M=Tp*Fs;n=0:M-1;A=444.128;al
12、ph=pi*50*20.5;omega=pi*50*20.5;xnt=A*exp(-alph*n*T).*sin(omega*n*T);Xk=T*fft(xnt,M); %M点FFTxnt)yn=xa(nT);subplot(3,2,1);tstem(xnt,yn);%调用自编绘图函数tstem绘制序列图box on;title(a) Fs=1000Hz);k=0:M-1;fk=k/Tp;subplot(3,2,2);plot(fk,abs(Xk);title(a) T*FTxa(nT),Fs=1000Hz);xlabel(f(Hz);ylabel(幅度);axis(0,Fs,0,1.2*ma
13、x(abs(Xk)% Fs=300Hz和 Fs=200Hz的程序与上面Fs=1000Hz完全相同。2 频域采样理论的验证程序清单%频域采样理论验证程序exp2b.mM=27;N=32;n=0:M;%产生M长三角波序列x(n)xa=0:floor(M/2); xb= ceil(M/2)-1:-1:0; xn=xa,xb;Xk=fft(xn,1024);%1024点FFTx(n), 用于近似序列x(n)的TFX32k=fft(xn,32);%32点FFTx(n)x32n=ifft(X32k);%32点IFFTX32(k)得到x32(n)X16k=X32k(1:2:N);%隔点抽取X32k得到X16
14、(K)x16n=ifft(X16k,N/2);%16点IFFTX16(k)得到x16(n)subplot(3,2,2);stem(n,xn,.);box ontitle(b) 三角波序列x(n);xlabel(n);ylabel(x(n);axis(0,32,0,20)k=0:1023;wk=2*k/1024;%subplot(3,2,1);plot(wk,abs(Xk);title(a)FTx(n);xlabel(omega/pi);ylabel(|X(ejomega)|);axis(0,1,0,200)k=0:N/2-1;subplot(3,2,3);stem(k,abs(X16k),.)
15、;box ontitle(c) 16点频域采样);xlabel(k);ylabel(|X_1_6(k)|);axis(0,8,0,200)n1=0:N/2-1;subplot(3,2,4);stem(n1,x16n,.);box ontitle(d) 16点IDFTX_1_6(k);xlabel(n);ylabel(x_1_6(n);axis(0,32,0,20)k=0:N-1;subplot(3,2,5);stem(k,abs(X32k),.);box ontitle(e) 32点频域采样);xlabel(k);ylabel(|X_3_2(k)|);axis(0,16,0,200)n1=0:
16、N-1;subplot(3,2,6);stem(n1,x32n,.);box ontitle(f) 32点IDFTX_3_2(k);xlabel(n);ylabel(x_3_2(n);axis(0,32,0,20)%=tstem 程序清单=function tstem(xn,yn)%时域序列绘图函数% xn:信号数据序列，yn:绘图信号的纵坐标名称（字符串）n=0:length(xn)-1;stem(n,xn,.);box onxlabel(n);ylabel(yn);axis(0,n(end),min(xn),1.2*max(xn)七、实验程序运行结果1 时域采样理论的验证程序运行结果exp
17、2a.m如图3.2所示。由图可见，采样序列的频谱的确是以采样频率为周期对模拟信号频谱的周期延拓。当采样频率为1000Hz时频谱混叠很小；当采样频率为300Hz时，在折叠频率150Hz附近频谱混叠很严重；当采样频率为200Hz时，在折叠频率110Hz附近频谱混叠更很严重。图2.22 时域采样理论的验证程序exp2b.m运行结果如图3.3所示。图2.3该图验证了频域采样理论和频域采样定理。对信号x(n)的频谱函数X(ej)在0，2上等间隔采样N=16时， N点IDFT得到的序列正是原序列x(n)以16为周期进行周期延拓后的主值区序列：由于NM，频域采样定理，所以不存在时域混叠失真，因此与x(n)相
18、同。实验三离散傅立叶变换及谱分析一、 实验目的1掌握离散傅里叶变换的计算机实现方法。2检验实序列傅里叶变换的性质。3掌握计算序列的圆周卷积的方法。4熟悉连续信号经理想采样前后的频谱变化关系，加深对时域采样定理的理解。5学习用DFT对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差，以便在实际中正确应用DFT。二、 实验内容1实现离散傅里叶变换。2计算序列圆周卷积。3计算实序列傅里叶变换并检验DFT性质。4实现连续信号傅里叶变换以及由不同采样频率采样得到的离散信号的傅里叶变换。5实现补零序列的傅里叶变换。6实现高密度谱和高分辨率谱，并比较二者的不同。三、 实验报告要求见各程序要求%
19、以下为4个扩展函数% （1）离散傅立叶变换 采用矩阵相乘的方法function Xk=dft(xn,N)n=0:1:N-1;k=0:1:N-1;WN=exp(-j*2*pi/N);nk=n*k;WNnk=WN.(nk);Xk=xn*WNnk;%（2）逆离散傅立叶变换function xn=idft(Xk,N)n=0:1:N-1;k=0:1:N-1;WN=exp(-j*2*pi/N);nk=n*k;WNnk=WN.(-nk);xn=(Xk*WNnk)/N;% (3) 实序列的分解 % 实序列可分解为共扼对称分量% 和共扼反对称分量function xec,xoc=circevod(x)N=len
20、gth(x);n=0:(N-1);xec=0.5*(x+x(mod(-n,N)+1); %根据上面的公式计算，mod函数为取余xoc=0.5*(x-x(mod(-n,N)+1);% (4) 序列的循环移位 function y=cirshftt(x,m,N)if length(x)N error(N mustbe = the length of x) %要求移位周期大于信号长度endx=x zeros(1,N-length(x);n=0:1:N-1;n=mod(n-m,N);y=x(n+1);%例1 本例检验实序列的性质DFTxec(n)=ReX(k)DFTxoc(n)=ImX(k)%设 x(
21、n)=10*(0.8).n 0=nNerror(N mustbe = the length of x1)endif length(x2)Nerror(N mustbe = the length of x2)endx1=x1 zeros(1,N-length(x1); %将x1,x2补0成为N长序列x2=x2 zeros(1,N-length(x2);m=0:1:N-1;x2=x2(mod(-m,N)+1); %该语句的功能是将序列翻褶，延拓，取主值序列H=zeros(N,N);for n=1:1:N %该for循环的功能是得到x2序列的循环移位矩阵 H(n,:)=cirshftt(x2,n-1
22、,N); %和我们手工计算圆周卷积得到的表是一致的endy=x1*H %用矩阵相乘的方法得到结果% 例3 本例验证采样定理%令，绘制其傅立叶变换。用不同频率对其进行采样，分别画出离散时间傅立叶变换。已给出采样频率为时的的程序 %实验报告要求：（1）请写出时的程序 (2)将实验结果画出 (3)实验结果说明什么（采样频率不同结果有何不同）。Dt=0.00005; %步长为0.00005st=-0.005:Dt:0.005;xa=exp(-1000*abs(t); %取时间从-0.005s到0.005s这段模拟信号Wmax=2*pi*2000; %信号最高频率为2*2000K=500; %频域正半轴
23、取500个点进行计算k=0:1:K;W=k*Wmax/K; % 求模拟角频率Xa=xa*exp(-j*t*W)*Dt; %计算连续时间傅立叶变换（利用矩阵运算实现） Xa=real(Xa); %取实部W=-fliplr(W),W(2:501); %将角频率范围扩展为从-到+Xa=fliplr(Xa),Xa(2:501);%A = 1 3 5 7 9 则fliplr(A)=9 7 5 3 1 subplot(2,2,1);plot(t*1000,xa); %画出模拟信号，横坐标为时间（毫秒），纵坐标为幅度xlabel(time(millisecond);ylabel(xa(t);title(an
24、olog signal);subplot(2,2,2);plot(W/(2*pi*1000),Xa*1000);%画出连续时间傅立叶变换 xlabel(frequency(kHZ); %横坐标为频率（kHz）ylabel(xa(jw); %纵坐标为幅度title(FT);%下面为采样频率5kHz时的程序Ts=0.0002; %采样间隔为n=-25:1:25;x=exp(-1000*abs(n*Ts); %离散时间信号K=500;k=0:1:K;w=pi*k/K; %w为数字频率X=x*exp(-j*n*w); %计算离散时间傅立叶变换（序列的傅立叶变换）X=real(X);w=-fliplr(
25、w),w(2:K+1);X=fliplr(X),X(2:K+1);subplot(2,2,3);stem(n*Ts*1000,x); %画出采样信号（离散时间信号）xlabel(time(millisecond);gtext(Ts=0.2ms); %该语句可以将引号中的内容放置在figure中的任何地方，只需 %将十字的中心放在想放置内容的地方，然后按鼠标即可。ylabel(x1(n);title(discrete signal);subplot(2,2,4);plot(w/pi,X); %画出离散时间傅立叶变换xlabel(frequency(radian); %横坐标为弧度ylabel(x
26、1(jw);title(DTFT);%例4 本例说明补零序列的离散傅立叶变换%序列，已给出序列的傅立叶变换程序和将原序列补零到10长序列的DFT%实验报告要求： （1）编写将序列补零到20长后计算DFT的程序(2)给出实验结果(3)实验结果说明什么（即序列补零后进行DFT,频谱有何变化）n=0:4;x=ones(1,5); %产生矩形序列k=0:999;w=(pi/500)*k;X=x*(exp(-j*pi/500).(n*k); %计算离散时间傅立叶变换Xe=abs(X); %取模subplot(3,2,1);stem(n,x);ylabel(x(n); %画出矩形序列subplot(3,2
27、,2);plot(w/pi,Xe);ylabel(|X(ejw)|); %画出离散时间傅立叶变换N=10;x=ones(1,5),zeros(1,N-5); %将原序列补零为10长序列n=0:1:N-1;X=dft(x,N); %进行DFTmagX=abs(X);k=(0:length(magX)-1)*N/length(magX);subplot(3,2,3);stem(n,x);ylabel(x(n); %画出补零序列subplot(3,2,4);stem(k,magX); %画出DFT结果axis(0,10,0,5);ylabel(|X(k)|);%例5 本题说明高密度谱和高分辨率谱之间
28、的区别，高密度谱是信号补零后得到的，虽然谱线相当密，但是因为信号有效长度不变，所以其分辨率也不变，因此还是很难看出信号的频谱成分。高分辨率谱是将信号有效长度加长，因此分辨率提高，可以看出信号的成分。%有一个序列为 （该序列周期计算可得40）%（1）下面给出有10个有效采样点序列的DFT程序%（2）请写出将第一问中的10长序列补零到40长，计算其DFT%（3）采样n=0：39，计算有40个有效采样点的序列的DFT%实验报告要求： （1）请编写将有10个有效采样点的序列补零到40长后计算DFT的程序 (2) 请编写计算有40个有效采样点的序列的DFT的程序 (3) 将实验结果画出并分析实验结果说明
29、什么M=10;n=0:M-1;x=2*cos(0.35*pi*n)+cos(0.5*pi*n);subplot(2,1,1);stem(n,x);title(没有足够采样点的信号);Y=dft(x,M);k1=0:M-1;w1=2*pi/M*k1;subplot(2,1,2);stem(w1/pi,abs(Y);title(信号的频谱);M=10;n=0:M-1;b=40;s=0:b-1;x=2*cos(0.35*pi*n)+cos(0.5*pi*n), zeros(1,30);subplot(2,1,1);stem(s,x);title(有足够采样点的信号);Y=dft(x,b);k1=0:
30、b-1;w1=2*pi/b*k1;subplot(2,1,2);stem(w1/pi,abs(Y);title(信号的频谱);实验四用FFT对信号作频谱分析一、实验目的1、理解用FFT对周期序列进行频谱分析时所面临的问题并掌握其解决方法。2、掌握用时域窗函数加权处理的技术。3、理解用FFT对非周期信号进行频谱分析所面临的问题并掌握其解决方法。二、实验原理与计算方法、对周期序列进行频谱分析应注意的问题k k(a)时域周期整数倍截断 (b)时域非周期整数倍截断图 4-1 周期函数的幅频曲线k图 4-2 矩形窗的频谱对时间序列作FFT时，实际上要作周期延拓（如果取长序列的一段进行计算还要先作截断）。
31、周期序列是无限长时间序列，如果截断区间刚好就是该序列周期的整数倍，那么在进行周期延拓后，将还原出原来的周期序列，由此可以较精确地计算出的该周期序列的频谱。反之，如果截断区间并不是该序列周期的整数倍，那么在进行周期延拓后，就不可能还原出原来的周期序列，由此计算出的频谱与该周期序列的频谱存在误差，而且误差的大小与截断区间的选取直接相关，如图4-1所示，其中幅度频谱的量值表示为，以dB(分贝)为单位。这种误差是由于周期序列与矩形截断序列相乘在频域产生二者的频谱卷积形成的。矩形窗的频谱是抽样函数序列，如图4-2所示。除了k = 0处主瓣内集中了大部分的能量外，两旁的较小峰值处的旁瓣也分散了一部分能量，它与周期序列频谱卷积的结果使原来集中的频谱展宽，称为频率泄漏。如果对已知周期的信号作频谱分析，在进行时域截断时，完全可以选取其周期的整数倍裁取，从而可以避免这种频率泄漏的发生。不过，通常需要进行频谱分析的信号是周期未知的信号，或随机信号，无法判断它的周期值，为了尽量避免频率泄漏对结果的影响，在作时间截断时，就应选取其频谱的旁瓣较小的截断函数，以减轻泄漏问题。2、时域窗函数的应用作为截断函数，矩形窗在作时间截断时，对所截取区间内的信号不加以任何影响