数字信号处理仿真实验

数字信号处理上机实验仿真报告学 院电子工程学院班 级学号姓名实验一：信号、系统及系统响应一、实验目的(1)熟悉连续信号经理想采样前后的频谱变化关系，加深对时域采样定理的理解。(2)熟悉时域离散系统的时域特性。(3)利用卷积方法观系统的时域特性。(4)掌握序列傅里叶变换的计算机实现方法，利用序列的傅里叶变换对连续信号、离散信号及系统响应进行频域分析。二、实验原理采样是连续信号数字处理的第一个关键环节对一个连续信号xt进行a理想采样的过程可(1.1式表示。 ˆxtxtpt (1.1) ˆˆ其中xt为xt的理想采样，pt为周期冲激脉冲，即a aˆpttnT (1.2)nˆxt的傅里叶变换XΩ为a aˆˆTXjΩ1XjΩjΩ (1.3)a a sˆTm将(1.2式代(1.1式并进行傅里叶变换， ˆXjΩxttTejtda ˆxttnTejΩtdtn axTejΩ (1.4)an式中的xT列xn，a

即xnxTaˆˆ Xexne (1.5)n比(1.5和(1.4可知XΩXe )aT为了在数字计算机上观察分析各种序列的频域特性，通常对Xe在,行M为N列xn有XeN1xnejk )n0中k,k,,M1kM一个时域离散线性时不变系统的输出关系为ynxnhnxmhnm (1.8)m上述卷积运算也可以转到频域实现Ye Xe He (1.9)三、实验内容及步骤(1)认真复习采样理论、离散信号与系统、线性卷积、序列的傅里叶变换及性质等有关内容，阅读本实验原理与方法。(2)编制实验用主程序及相应子程序。①信号产生子程序，用于产生实验中要用到的下列信号序列：xa(t)=Ae-atsinΩ0t)u(t)进行采样，可得到采样序列xa(n)=xa(nT)=Ae-anTsinΩ0nT)u(n), ≤n<50其中A为幅度因子，a为衰减因子，Ω0是模拟角频率，T为采样间隔。这些参数都要在实验过程中由键盘输入，产生不同的xa(t)和xa(n)。b.单位脉冲序列：xb(n)=δ(n)c.矩形序列：xc(n)=RN(n),N=10②系统单位脉冲响应序列产生子程序。本实验要用到两种FIR系统。a.ha(n)=R10(n);b.hb(n)=δ(n)+2.5δ(n-1)+2.5δ(n-2)+δ(n-3)③有限长序列线性卷积子程序，用于完成两个给定长度的序列的卷积。可以直接调用MATLAB语言的积数conv。cnv用于两有长序的卷积，它定个列从n=0开始。调用格如：y=ov(,h)(3)调并行验序，完成下实内：①分析样列特。a.取采样率fs1kz,即T=1s。b.改变采样频率，fs=300Hz，观察|X(ejω)|的变化，并做记录(打印曲)；进一步降低采样频率，fs=200Hz，观察频谱混叠是否明显存在，说明原因，并记录(打印)这时的|X(eω)|曲线。②时域离散信号、系统和系统响应分析。a.观察信号xb(n)和系统hb(n的时域和频域特性；利用线性卷积求信号xb(n通过系统hb(n)的响应y(n)，比较所求响应y(n和hb(n的时域及频域特性，注意它们之间有无差别，绘图说明，并用所学理论解释所得结果。b.观察系统ha(n)对信号xc(n)的响应特性。③卷积定理的验证。（4）主程序框图①分析采样序列的特性②时域离散信号、系统和系统响应分析始生和号产生X(w)生X)生束开始生生调用傅氏变换子程序生Xb(ew,bew)和Y(ej)产号像束③卷积定理的验证开始生xb)()利用卷积公式生y(n)调用傅氏变换子程序生Xb(ew,bew和Y(ejw)算YXb出结束四．实验程序及对应波形1.子程序function[XN,n,k]=DFT(xn,N)n=0:N-1;k=-200:200;XN=xn*exp(-j*2*pi/N).^(n'*k); 计算DFT[x(n)]%产生矩形序列functionx=juxing(n2);x=[1,ones(1,n2)];function[x,n]=maichong(n0,n1,n2)n=(n1:n2);x=(n==n0);%产生信号Xa(n)functionx=xn(A,a,w,fs)n=0:50-1;x=A*exp((-a)*n/fs).*sin(w*n/fs).*juxing(49);functionx=u(t);x=(t>=0);%产生脉冲信号function[x,n]=maichong(n0,n1,n2)n=(n1:n2);x=(n==n0);2.主程序①分析采样序列的特性程序代码如下（见后附程序Ex1_1）:A=100;a=200;w0=200;k=-200:200;T=0.001;t=0:T:0.06;N=50;k1=0:1:N;W1max=2*pi*500;W1=W1max*k1/N;w1=W1/pi;xat=A*exp(-a*t).*sin(w0*t).*u(t);Xa=xat*exp(-j*t'*W1);figure(1);subplot(2,1,1);plot(t,xat);xlabel('t');ylabel('xa(t)');title('采样信号a(t)');axis([0,0.06,-5,35]);subplot(2,1,2);plot(w1,abs(Xa));xlabel('w');ylabel('X(jw)');title('xa(t)的频谱');%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%fs=1000HzA=100;a=200;w0=200;k=-200:200;;;n=;%产号x);));,'');('');('');('采号f');;;('');ylabel'X(e^jw));title'xa(n的频谱(fs=1000Hz));%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%fs=300Hz;xan=xn(A,a,w0,fs);%产生信号xa(n)X=DFT(xan,50);figure(3)subplot(2,1,1)n=0:49;stem(n,xan,'.');axis([0,50,-20,50]);xlabel('n');ylabel('xa(n)');title('采样信号fs=300Hz');subplot(2,1,2);plot(w,abs(X));xlabel('w/pi');ylabel'X(e^jw));title'xa(n的频谱(fs=300Hz));%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%fs=200Hz;n=;%产号x);));,'');('');('');('采号f');;;('');ylabel'X(e^jw));title'xa(n的频谱(fs=200Hz));仿真结果如:由图可见在折叠频率w=π,即f=fs/2=500Hz处混叠很小当fs=300Hz时，存在较明显的混叠失真；当fs=200时，发生严重的混叠失真。②时域离散信号、系统和系统响应分析a.观察信号xb(n)和系统hb(n)的时域和频域特性；序Ex1_2）%信号xbN=64xb=zeros(1,N)xb(1)=1;xn1=[0:N-1]figure(1);subplot(2,1,1)stem(xn1,xb'.')title'信号xb时域特性')w=linspace(-5*pi,5*pi,1000)Xb=xb*exp(-j*xn1'*w)figure(1)subplot(2,1,2)plot(w/pi,abs(Xb)'.')title'信号xb频域特性')%ϵÍ信号hbhb=[1,2.5,2.5,1];nhb=[0:3]figure(2);subplot(2,1,1)stem(nhb,hb'.')title'系统函ýhb时域特性')w=linspace(-2.8*pi,2.8*pi,1000)Hb=hb*exp(-j*nhb'*w)figure(2)subplot(2,1,2)plot(w/pi,abs(Hb)'.')title'系统函ýhb频域特性')%响应yb=conv(xb,hb);nyb=[0:length(yb)-1]figure(3);subplot(2,1,1)stem(yb'.')title('响应yb时域特性')w=linspace(-2.8*pi,2.8*pi,1000)Yb=yb*exp(-j*nyb'*w)figure(3)subplot(2,1,2)plot(w/pi,abs(Yb),'.')title('响应Yb频域特性')仿真结果如下:b.观察系统ha(n)对信号xc)的响应特性。实验代码如下(见后附程序Ex1_3):%号xc=0c)n]i)t,'')i('信号xc时应')=)c)i)u)l)ϳ信号ha=0a)n]i)t,'')i('系数ha时应')=0)a)i)u)l,'')i('系数ha频应')%应2;y]i)t,'')i('响应y2时应')=0)2)figure(3)subplot(2,1,2)plot(w/pi,abs(Y2),'.')title('响应Yb频域响应')仿真结果如下:③卷积定理的验证。序%时域卷积定理的验证号N=10w=linspace(-5*pi,5*pi,1000)%系统信号haN=10ha=ones(,N)hn1=[0:N-1]w=linspace(-5*pi,5*pi,1000)%响应w=linspace(-5*pi,5*pi,1000)figure(1)subplot(2,1,1)plot(w/pi,abs(Y2))应Y=Xc.*Hasubplot(2,1,2)plot(w/pi,bs(Y))title('响应Y=Xc.a频域特性')仿真结果如下：两种方法得到的谱是一样的，卷积定理。五、思考题(1)在分析理想采样序列特性的实验中，采样频率不同时，相应理想采样序列的傅里叶变换频谱的数字频率度量是否都相同?它们所对应的模拟频率是否相?为什?答：由T可知，若采样频率不同，则其周期T不同，相应的数字频率也不相同；而因为是同一信号，故其模拟频率保持不变。(2)在卷积定理验证的实验中，如果选用不同的频域采样点数M值，例如，选M=10和M=20，分别列的傅里叶变换，得Y(ek)X(ejk)H(ek),k0,1,,M1a所得结果之间有无差异?为什?得Y(ejk)定，由频域采若M小于采样序列的长度N时域混叠现象。实验二：用FFT作谱分析一、实验目的：（1）进一步加深DFT算法原理和基本性质的理（因为FFT是DFT的一种快速算法，所以FFT的运算结果必然满足DFT的基本性质。（2）熟悉FFT算法原理和FFT子程序的应用。（3）学习用FT对连续信号和时域离散信号进行谱分析的方，了解可能出现误差及其原因便在实际中确应用FFT。二、实验步骤：（1）复习DFT的定义、性质和用DFT作谱分析的有关内容。（2）复习FFT算法原理与编程思想并对照DITFFT运算流图和程序框图，读懂本实验提供的FFT子程。（3）编制信号产生子程序，产生以下典型信号供谱分析用：xn=R4(n)1n1,0nxn8n,4n720,4n,0nxnn3,4n730,其它xn=cos(pi/4*n)4xn=sin(pi/8*n)5xt=cos(pi*8*t)+sin(pi*16*t)+cos(20*pi*t)6应当注意，如果给出的是连续信号xa(t),则首先要根据其最高频率采样速率fs以及由频择采样点数N，然后对其进行软件采样（即计算x(n)=xa(nT),0<=n<=N-1,产生对序列x(n)。对信号x6(t)应则有可能产生较大的分析误差实验根据DFT的隐含周期性思考这个问题。（4）编写主程序。下图给出程序框图，供参考。验提供FFT子程序和通用绘图子程序。始读入度N生信号图用T子程序（函数）计算信号的DTF调用绘图数制Xk）线束三．实验程序及对应波形(1)对2中所给出的信号逐行谱分析。实验代码如下（见后附程序Ex2_1）号1x=[,,,1]n=03]fgr(1;upo(,11)s(1)tite号x1时性)w=isac(5pi5p,00)X=x*x(jn'w)fgr(1)sblt(,,)potwp,bsX))T分析信号x1域性)fgr(1)sblt(,,))t号性)号n=07]fgr(2;upo(,11)号性)w=isac(5pi5p,00)fgr(2)sblt(,,)potwp,bsX))T号域性)fgr(2)sblt(,,)smas2))t号性)号3211234]n=07]fgr(3;upo(,11)号性)w=isac(5pi5p,00)fgr(3)sblt(,,)potwp,bsX))T号性)fgr(3)sblt(,,)smb())t号性)号n=[:5]fgr(4;upo(,11)号性)w=isac(5pi5p,00)fgr(4)sblt(,,)potwp,bsX))T号性)fgr(4)sblt(,,)smas4))t号性)号n=[:5]fgr(5;upo(,11)号性)w=isac(5pi5p,00)fgr(5)sblt(,,)potwp,bsX))fgr(5)sblt(,,)smas5))t号性)号t0000:5fgue();ublt(,11);线2.555])T0.1n=025T;号fgue();ublt(,12);

号线%建离变量%采样周期为T号.5-5])wlnpc(*p5p,50;fgue();ublt(,13);

%画散信线%点%计算DTFTpt/ibX)teDFT分析样信频谱;真果下：(2)令用FFT算8点和16点离散傅里叶变(3)令重复(2。序号n=[:5]fgr(1;upo(,11)号性)w=isac(5pi5p,00)fgr(1)sblt(,,)potwp,bsX))T号域性)号n=[:5]fgr(2;upo(,11)T号域性)w=isac(5pi5p,00)fgr(2)sblt(,,)potwp,bsX))fgr(3;upo(,11)时域特性)n8:7fgr(3;upo(,12)se(n_8asXn_8)的8点t换)n_60:5figure(3);subplot(3,1,3)se(11,b(n_6)title('xn1的16点t换)%x(n)=x4(n)+jx5(n)xn2=x4+j*x5nx2=0:length(xn2)-1figure(4);subplot(3,1,1)stem(nx2,xn2)title('xn2=x4+j*x5性)Xn2_8=fftshift(fft(xn2,8))n8:7figure(4);subplot(3,1,2)se(n_8asXn_8)title('xn2的8点t换)Xn2_16=fftshift(fft(xn2,16))n_60:5figure(4);subplot(3,1,3)se(21,b(n_6)title('xn2的16点t换)真果下：四、思考题(1)在8，x2(n)和x3(n)的幅频特性会相?为?N=16?答在=8x2(n)和x3(n)的幅频特性会相同在N8时x2(n)和x3(n)的幅频不相同。(2)如果周期信号预先不知道，如何用FFT进行谱分?答：如果周期信号的周期预先道。实验三：用窗函数法设计FIR数字滤波器一、实验目的(1)掌握用窗函数法设计FIR数字滤波器的原理和方法。(2)熟悉线性相位FIR数字滤波器特性。(3)了解各窗函滤波特性的响。二、实验原理如果所希望的滤波器的理想频数为H(ej)则其对应d的单位脉冲响应为=h(n) 1=d 2

H(ej)ejd d

(2-1)d dw(n)h(n)h(n)窗函数设计法的基本原理是用有限长单位脉冲响应序列h(n)逼近h(n)。由于h(n)往往是无限长序列，且是非因果的，所以用窗函d dw(n)h(n)h(n)d ＝d (2-2)h(n)就作为实际设计的FIR数字滤波器的单位脉冲响应序列率响应函数H(ej)为：H(ej)＝N1h(n)e (2-3)n0式中，N为所选窗函数w(n)的长度。由第七章可知用窗函数法设计的滤波器性能取决于窗函数(n)的类型及窗口长度N的取值设计过程中要根据对阻带最小衰减和过渡带宽度的要求选择合适的窗函数类型和窗口长度N的窗函数可达到的小衰减和过渡带宽度见第七章。d这样选定窗函数类型和长度N后，求出单位脉冲响应h(n)＝h(n)·(n)，并按式（2-3）求出H(ej)。H(ej)是否满足要求，要d在h(n)尾部加零使长度满足于2用FFT计算H(ej)要观察补零点数增可如果H(ej)不满足要求，则要重新选择窗函数类型和长度N足。如果要求线性相位特性，则h(n)还必须满足h(N1(2-4)根据上式中的正负号和长度N的奇偶性又将性相位FIR滤波设计线性低通特征，可选择h(N1一类，而不能选h(N1n)三、实验内容及步骤(1)复习用窗法设计R数字滤波器一节内容，阅读本实验原理，掌握设计步骤。(2)编写程序。①编写能产生矩型窗、升余弦窗、改进升余弦窗和二阶升余弦窗的窗函数子程序。②编写主程序。其中幅度特性要求用dB表示。三、实验代码及仿真结果序fp200;f=00;ffcpfs;fs5000;ffcpfs;=f;;wc(fpfc)fs;hf1,clca(+);n0:nghh-;w=ispce-i,i100);H=0lg1(bsH);fgr(1)t(')titl'n域性线)fgr(2)sbpot21,)powp,);xlabel('\omega/π');ylabel('|H(e^omega)|/dB')tl'n幅度谱;sbpot21,)po(/i,gl(1p);xlabl\mgπ');ylabel('相位)tl'n相位谱;仿真结果如:N=inut宽度N=');k=input(型1矩窗弦窗3窗4窗择;sublot2,,1);wp/;a(N1/;n0:N1);=as;=nm/m;ifk==1Bbret);eleifk==2Bhann();eleifk==3Bamg);eleifk==4Bbaka();ededddhdh.(';s(d';xbln;yae'(';tite在矩窗下的N=64时h(n)数;sublot2,,2);[,m]feqh,[]104);ma=b();d2l1m+pm(g;pom/,b;xal'/);lbl2l[(j));tl'n幅频性;gi;ph=anleH);sbpo(,,);lma;xbl';lb(';tl'n相频性;sublot2,,4);,;xlbe(w';ylael'(ej)';tl'n幅频性;仿真结果如:五、思考题(1)如果给定通带截止频率和阻带截止频率以及阻带最小衰减，如何用窗法设计线性相位滤波器?写出设计步骤。答：①首先确定模拟带通滤波器的技术指标。②确定归一化低通技术要求③设计归一化低通④将低通转化通(2)如果要求函数法设计带通器，且给定上、下边带截止频率为ω1和2，试求理想带通的单位脉冲响应hd(n)。答：①首先确定所需数字滤波器的类型，确定线性相位数字滤波器的类型。②选择合适的窗函数w(n)③确定数字滤波器的频应函数H(e)H(w)ed()d d④计算理想带通的单位脉冲响应hd()h(n)d

1

2H(eendd1实验四：用双线性变换法设计IIR数字滤波器一、实验目的(1)掌握用双线性变换法设计IIR数字滤波器的原方法。(2)熟悉IIR数字滤波器特性。二、实验内容及步骤(1)复习用双线性变换法设计IIR数字滤波器一节内容，阅读本实验原理，掌握设计步骤。(2)编写程序。①设计产生巴特沃斯滤波器原型，采用双线性变换法设计巴特沃斯数字低通滤波器。②编写主程序。其中幅度特性要求用dB表示。三、上机实验内容(1)实验所需设计的R数字低通滤波器的指标为wp=0.2*pi，ws=0.*pi，Rp=1dB，As=15dB，滤波器采样频率为100Hz。(2)采用巴特滤波器为原型进计。(3)将相指求的IIR滤波器和FR滤波器进行比较。四、实验代码及仿真结果(见后附程序Ex4_1)fs=00;T=1fs;Wp=.2pi;Ws=.3pi;Ap=;As=5;fs=s/s;

采样率%设置归一化带和阻截止率%设置通带大和最衰减[N,c]butor(W,W,ApAs's); %调用butter函数确定巴特沃斯滤波器阶数[B,A]=butter(N,Wc,'s');W=linspace(0,pi,400*pi);hf=freqs(B,A,W);hf=20*log10(abs(hf));subplot(2,1,1);plot(W/pi,hf);线gridon;tt(IR巴特沃斯模拟滤波器');xlabel(频率/Hz');ylabel(幅度');[D,C]=bilinear(B,A,fs);Hz=freqz(D,C,W);Hz=20*log10(abs(Hz));subplot(2,1,2);

调用butter函数设计巴特沃斯滤波器指定一段频率值%计算模拟滤波器的幅频响应%绘出巴特沃斯模拟滤波器的幅频特性曲调用双线性变换法%返回频率响应plot(W/pi,Hz); 绘出巴特沃斯数字低通滤波器的幅频特性曲