数字信号处理实验报告(全)

实验一、离散时间系统及离散卷积1、单位脉冲响应源程序：function pr1() %定义函数pr1a=[1,-1,0.9]; %定义差分方程y(n)-y(n-1)+0.9y(n-2)=x(n) b=1;x=impseq(0,-40,140); %调用impseq函数n=-40:140; %定义n从-40 到140h=filter(b,a,x); %调用函数给纵座标赋值figure(1) %绘图figure 1 (冲激响应) stem(n,h); %在图中绘出冲激title('冲激响应'); %定义标题为:'冲激响应'xlabel('n'); %绘图横座标为nylabel('h(n)'); %绘图纵座标为h(n)figure(2) %绘图figure 2[z,p,g]=tf2zp(b,a); %绘出零极点图zplane(z,p)function [x,n]=impseq(n0,n1,n2)%声明impseq函数n=[n1:n2];x=[(n-n0)==0];结果：Figure 1:Figure 2:2、离散系统的幅频、相频的分析源程序：function pr2()b=[0.0181,0.0543,0.0543,0.0181];a=[1.000,-1.76,1.1829,-0.2781];m=0:length(b)-1; %m从0 到3l=0:length(a)-1; %l从0 到3K=5000;k=1:K;w=pi*k/K; %角频率wH=(b*exp(-j*m'*w))./(a*exp(-j*l'*w));%对系统函数的定义magH=abs(H); %magH为幅度angH=angle(H); %angH为相位figure(1)subplot(2,1,1); %在同一窗口的上半部分绘图plot(w/pi,magH); %绘制w(pi)-magH的图形grid;axis([0,1,0,1]); %限制横纵座标从0到1xlabel('w(pi)'); %x座标为 w(pi)ylabel('|H|'); %y座标为 angle(H)title('幅度，相位响应'); %图的标题为:'幅度，相位响应' subplot(2,1,2); %在同一窗口的下半部分绘图plot(w/pi,angH); %绘制w(pi)-angH的图形grid; %为座标添加名称xlabel('w(pi)'); %x座标为 w(pi)ylabel('angle(H)'); %y座标为 angle(H)结果：3、卷积计算源程序：function pr3()n=-5:50; %声明n从-5到50u1=stepseq(0,-5,50); %调用stepseq函数声用明u1=u(n)u2=stepseq(10,-5,50); %调用stepseq函数声用明u2=u(n-10) %输入x(n)和冲激响应h(n)x=u1-u2; %x(n)=u(n)-u(n-10)h=((0.9).^n).*u1; %h(n)=0.9^n*u(n)figure(1)subplot(3,1,1); %绘制第一个子图stem(n,x); %绘制图中的冲激axis([-5,50,0,2]); %限定横纵座标的范围title('输入序列'); %规定标题为:'输入序列'xlabel('n'); %横轴为nylabel('x(n)'); %纵轴为x(n)subplot(3,1,2); %绘制第二个子图stem(n,h); %绘制图中的冲激axis([-5,50,0,2]); %限定横纵座标的范围title('冲激响应序列'); %规定标题为:'冲激响应序列'xlabel('n'); %横轴为nylabel('h(n)'); %纵轴为h(n)%输出响应[y,ny]=conv_m(x,n,h,n); %调用conv_m函数subplot(3,1,3); %绘制第三个子图stem(ny,y);axis([-5,50,0,8]);title('输出响应'); %规定标题为:'输出响应'xlabel('n');ylabel('y(n)'); %纵轴为y(n)%stepseq.m子程序%实现当n>=n0时x(n)的值为1function [x,n]=stepseq(n0,n1,n2)n=n1:n2;x=[(n-n0)>=0];%con_m的子程序%实现卷积的计算function [y,ny]=conv_m(x,nx,h,nh)nyb=nx(1)+nh(1);nye=nx(length(x))+nh(length(h));ny=[nyb:nye];y=conv(x,h);结果：实验二、离散傅立叶变换与快速傅立叶变换1、离散傅立叶变换（DFT）源程序：function pr4()F=50;N=64;T=0.000625;n=1:N;x=cos(2*pi*F*n*T); %x(n)=cos(pi*n/16)subplot(2,1,1); %绘制第一个子图x(n)stem(n,x); %绘制冲激title('x(n)'); %标题为x(n)xlabel('n'); %横座标为nX=dft(x,N); %调用dft函数计算x(n)的傅里叶变换magX=abs(X); %取变换的幅值subplot(2,1,2); %绘制第二个子图DFT|X|stem(n,X);title('DFT|X|');xlabel('f(pi)'); %横座标为f(pi)%dft的子程序%实现离散傅里叶变换function [Xk]=dft(xn,N)n=0:N-1;k=0:N-1;WN=exp(-j*2*pi/N);nk=n'*k;WNnk=WN.^nk;Xk=xn*WNnk;结果：F=50,N=64,T=0.000625时的波形F=50,N=32,T=0.000625时的波形:2、快速傅立叶变换（FFT）源程序：%function pr5()F=50;N=64;T=0.000625;n=1:N;x=cos(2*pi*F*n*T); %x(n)=cos(pi*n/16) subplot(2,1,1);plot(n,x);title('x(n)');xlabel('n'); %在第一个子窗中绘图x(n)X=fft(x);magX=abs(X);subplot(2,1,2);plot(n,X);title('DTFT|X|');xlabel('f(pi)'); %在第二个子图中绘图x(n)的快速傅%里叶变换结果：3、卷积的快速算法源程序：function pr6()n=0:15;x=1.^n;h=(4/5).^n;x(16:32)=0;h(16:32)=0;%到此 x(n)=1, n=0~15; x(n)=0,n=16~32% h(n)=(4/5)^n, n=0~15; h(n)=0,n=16~32subplot(3,1,1);stem(x);title('x(n)');axis([1,32,0,1.5]); %在第一个子窗绘图x(n)横轴从1到32,纵轴从0到1.5 subplot(3,1,2);stem(h);title('h(n)');axis([1,32,0,1.5]); %在第二个子窗绘图h(n)横轴从1到32,纵轴从0到1.5 X=fft(x); %X(n)为x(n)的快速傅里叶变换H=fft(h); %H(n)为h(n)的快速傅里叶变换Y=X.*H; %Y(n)=X(n)*H(n)%Y=conv(x,h);y=ifft(Y); %y(n)为Y(n)的傅里叶反变换subplot(3,1,3) %在第三个子窗绘图y(n)横轴从1到32,纵轴从0到6 stem(abs(y));title('y(n=x(n)*h(n))');axis([1,32,0,6]);结果：实验三、IIR数字滤波器设计源程序：function pr7()wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;T=1;Fs=1/T;OmegaP=(2/T)*tan(wp/2); %OmegaP(w)=2*tan(0.1*pi) OmegaS=(2/T)*tan(ws/2); %OmegaS(w)=2*tan(0.15*pi)ep=sqrt(10^(Rp/10)-1);Ripple=sqrt(1/(1+ep.^2));Attn=1/10^(As/20);N=ceil((log10((10^(Rp/10)-1)/(10^(As/10)-1)))/(2*log10(OmegaP/OmegaS) ));OmegaC=OmegaP/((10.^(Rp/10)-1).^(1/(2*N)));[cs,ds]=u_buttap(N,OmegaC);[b,a]=bilinear(cs,ds,Fs);[mag,db,pha,w]=freqz_m(b,a);subplot(3,1,1); %在第一个子窗绘制幅度响应的图形plot(w/pi,mag);title('幅度响应');xlabel('w(pi)');ylabel('H');axis([0,1,0,1.1]);set(gca,'XTickmode','manual','XTick',[0,0.2,0.35,1.1]);set(gca,'YTickmode','manual','YTick',[0,Attn,Ripple,1]);grid;subplot(3,1,2); %在第二个子窗以分贝为单位绘制幅度响应的图形plot(w/pi,db);title('幅度响应(dB)');xlabel('w(pi)');ylabel('H');axis([0,1,-40,5]);set(gca,'XTickmode','manual','XTick',[0,0.2,0.35,1.1]);set(gca,'YTickmode','manual','YTick',[-50,-15,-1,0]);grid;subplot(3,1,3); %在第三个子窗绘制相位响应的图形plot(w/pi,pha);title('相位响应');xlabel('w(pi)');ylabel('pi unit');%axis([0,1,0,1.1]);set(gca,'XTickmode','manual','XTick',[0,0.2,0.35,1.1]);set(gca,'YTickmode','manual','YTick',[-1,0,1]);grid;function [b,a]=u_buttap(N,OmegaC)[z,p,k]=buttap(N);p=p*OmegaC;k=k*OmegaC.^N;B=real(poly(z));b0=k;b=k*B;a=real(poly(p));function [mag,db,pha,w]=freqz_m(b,a)[H,w]=freqz(b,a,1000,'whole');H=(H(1:501))';w=(w(1:501))';mag=abs(H);db=20*log10((mag+eps)/max(mag));pha=angle(H);结果：实验四、FIR数字滤波器的设计源程序：function pr8()wp=0.2*pi;ws=0.3*pi;tr_width=ws-wp;M=ceil(6.6*pi/tr_width)+1;n=0:M-1;wc=(ws+wp)/2;alpha=(M-1)/2;m=n-alpha+eps;hd=sin(wc*m)./(pi*m);w_ham=(hamming(M))';h=hd.*w_ham;[mag,db,pha,w]=freqz_m(h,[1]);delta_w=2*pi/1000;Rp=-(min(db(1:wp/delta_w+1)));As=-round(max(db(ws/delta_w+1:501)));subplot(2,2,1);stem(n,hd);title('理想冲激响应');axis([0,M-1,-0.1,0.3]);ylabel('hd(n)');subplot(2,2,2);stem(n,h);title('实际冲激响应');axis([0,M-1,-0.1,0.3]);ylabel('h(n)');subplot(2,2,3);plot(w/pi,pha);title('滤波器相位响应');axis([0,1,-pi,pi]);ylabel('pha');set(gca,'XTickmode','manual','XTick',[0,0.2,0.3,1.1]); set(gca,'YTickmode','manual','YTick',[-pi,0,pi]); grid;subplot(2,2,4);plot(w/pi,db);title('滤波器幅度响应');axis([0,1,-100,10]);ylabel('H(db)');set(gca,'XTickmode','manual','XTick',[0,0.2,0.3,1.1]); set(gca,'YTickmode','manual','YTick',[-50,-15,0]);function [mag,db,pha,w]=freqz_m(b,a)[H,w]=freqz(b,a,1000,'whole');H=(H(1:501))';w=(w(1:501))';mag=abs(H);db=20*log10((mag+eps)/max(mag));pha=angle(H);结果：。