《数字信号处理实验报告.docx》由会员分享，可在线阅读，更多相关《数字信号处理实验报告.docx（9页珍藏版）》请在三一办公上搜索。
1、数字信号处理实验报告实验一、离散时间系统及离散卷积 1、单位脉冲响应 源程序： function pr1 %定义函数pr1 a=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 到140 h=filter(b,a,x); %调用函数给纵座标赋值 figure(1) %绘图figure 1 (冲激响应) stem(n,h); %在图中绘出冲激 title(冲激响应); %定义标题为:冲激响应 xlabel(n); %绘图横座标为n ylabel
2、(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 到3 l=0:length(a)-1
3、; %l从0 到3 K=5000; k=1:K; w=pi*k/K; %角频率w H=(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到1 xlabel(w(pi); %x座标为 w(pi) ylabel(|H|); %y座标为 angle(H) title(幅度，
4、相位响应); %图的标题为:幅度，相位响应 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到50 u1=stepseq(0,-5,50); %调用stepseq函数声用明u1=u(n) u2=stepseq(10,-5,50); %调用stepseq函数声用明u2=u(n-10)
5、%输入x(n)和冲激响应h(n) x=u1-u2; %x(n)=u(n)-u(n-10) h=(0.9).n).*u1; %h(n)=0.9n*u(n) figure(1) subplot(3,1,1); %绘制第一个子图 stem(n,x); %绘制图中的冲激 axis(-5,50,0,2); %限定横纵座标的范围 title(输入序列); %规定标题为:输入序列 xlabel(n); %横轴为n ylabel(x(n); %纵轴为x(n) subplot(3,1,2); %绘制第二个子图 stem(n,h); %绘制图中的冲激 axis(-5,50,0,2); %限定横纵座标的范围 tit
6、le(冲激响应序列); %规定标题为:冲激响应序列 xlabel(n); %横轴为n ylabel(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)的值为1 function x,n=stepseq(n0,n1,n2) n=n1:n2; x=(n-n0)=0;
7、 %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、 离散傅立叶变换 源程序： 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);
8、 %xlabel(n); %X=dft(x,N); %magX=abs(X); %subplot(2,1,2); %stem(n,X); title(DFT|X|); xlabel(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时的波形绘制冲激 标题为x(n) 横座标为n 调用dft函数计算x(n)的傅里叶变换取变换的幅值 绘制第二个子图DFT|X| 横座标为
9、f(pi) F=50,N=32,T=0.000625时的波形: 2、 快速傅立叶变换 源程序： %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、
10、卷积的快速算法 源程序： function pr6 n=0:15; x=1.n; h=(4/5).n; x(16:32)=0; h(16:32)=0; %到此 x(n)=1, n=015; x(n)=0,n=1632 % h(n)=(4/5)n, n=015; h(n)=0,n=1632 subplot(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)横轴
11、从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
12、=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(
13、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); %在第二个子窗以分贝为单位绘制
14、幅度响应的图形 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,
15、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=ab
16、s(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
17、/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
18、(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); 结果：
