数字信号处理的MATLAB实现

昆明理工大学信息工程与自动化学院学生实验报告（2011—2012 学年第二学期）课程名称：数字信号处理开课实验室：信自楼111 2012 年 5 月 31 日年级、专业、班生医学号姓名成绩实验项目名称数字信号处理的matlab 实现指导教师教师评语教师签名：年月日一．实验目的熟练掌握matlab的基本操作。

了解数字信号处理的MATLAB实现。

二．实验设备安装有matlab的PC机一台。

三．实验内容.1.求信号x(n)=cos(6.3Пn/3)+cos(9.7Пn/30)+cos(15.3Пn/30),0≤n≤29的幅度频谱.2. 用冲击响应不变法设计一个Butterworth低通数字滤波器,要求参数为:Wp=0.2Пαp=1dBWs=0.3Пαs=15dB3.用双线性变换法设计一个Chebyshev高通IIR滤波器,要求参数为:Wp=0.6Пαp=1dBWs=0.4586Пαs=15dB4.用窗函数法设计一个低通FIR滤波器,要求参数为:Wp=0.2Пαp=0.3dBWs=0.25Пαs=50dB5.用频率抽样法设计一个带通FIR滤波器,要求参数为:W1s=0.2П W1p=0.35П W2p=0.65П W2s=0.8Пαs=60dB αp=1dB6.根据 4 点矩形序列，( n ) = [1 1 1 1] 。

做 DTFT 变换，再做 4 点 DFT 变换。

然后分别补零做 8 点 DFT 及 16 点 DFT。

7.调用filter解差分方程，由系统对u(n)的响应判断稳定性8编制程序求解下列系统的单位冲激响应和阶跃响应。

y[n]+ 0.75y[n -1]+ 0.125y[n -2] = x[n]- x[n -1]四．实验源程序1.n=[0:1:29];x=cos(6.3*pi*n/30)+cos(9.7*pi*n/30)+cos(15.3*pi*n/30);subplot(2,2,1);stem(n,x);title('信号x(n),0<=n<=29');X=fft(x);magX=abs(X(1:1:16));k=0:1:15;subplot(2,2,2);stem(k,magX);title('FFT幅度');2.wp=0.2*pi; %通带边界数字频率（Hz）ws=0.3*pi; %阻带边界数字频率（Hz）Rp=1; %通带波动（dB）As=15; %阻带波动（dB）T=1; %取T=1OmegaP=wp*T; %模拟原型通带边界频率OmegaS=ws*T; %模拟原型阻带边界频率[cs,ds]=afd_butt(OmegaP,OmegaS,Rp,As);%调用模拟Butterworth滤波器的设计[b,a]=imp_invr(cs,ds,T); %调用模拟到数字滤波器转换函数%[b0,B,A]=dir2cas(b,a); %直接型转换为级联型系数的计算[H,w]=freqz(b,a,501); %调用计算频率特性函数mag=abs(H);pha=angle(H);subplot(2,2,1);plot(w/pi,mag);title('幅频特性');subplot(2,2,2);plot(w/pi,-pha/pi);title('相频特性');3.ws=0.4586*pi; %通带边界数字频率（Hz)wp=0.6*pi; %阻带边界数字频率（Hz）Rp=1; %通带波动（dB）As=15; %阻带衰减（dB）[N,wn]=cheb1ord(wp/pi,ws/pi,Rp,As);%Chebyshew滤波器参数计算[b,a]=cheby1(N,Rp,wn,'high'); %数字Chebyshew高通滤波器设计%[b0,B,A]=dir2cas(b,a); %直接型转换为级联型系数计算[H,w]=freqz(b,a,501); %调用计算频率特性函数mag=abs(H);pha=angle(H); %计算幅频和相频特性subplot(2,2,1);plot(w/pi,mag);title('幅频特性');subplot(2,2,2);plot(w/pi,-pha/pi);title('相频特性');4.p=0.2*pi;ws=0.3*pi;tr_width=ws-wp;M=ceil(6.6*pi/tr_width)+1;n=[0:1:M-1];wc=(ws+wp)/2;hd=ideal_lp(wc,M);w_ham=(hamming(M))';h=hd.*w_ham;[H,w]=freqz(h,[1],501);mag=abs(H);pha=angle(H);db=20*log10((mag+eps)/max(mag));delta_w=2*pi/1000;Rp=-(min(db(1:1:wp/delta_w+1)));As=-round(max(db(ws/delta_w+1:1:501)));subplot(2,2,1);stem(n,h);title('单位抽样响应');subplot(2,2,2);plot(w/pi,db);title('幅频特性(db)');Rp=-(min(db(1:1:wp/delta_w+1)))As=-round(max(db(ws/delta_w+1:1:501)))5.M=40;alpha=(M-1)/2;l=0:M-1;wl=(2*pi/M)*l;T1=0.109021;T2=0.59417456;Hrs=[zeros(1,5),T1,T2,ones(1,7),T2,T1,zeros(1,9),T1,T2,ones(1,7),T2,T1,zeros(1, 4)];Hdr=[0,0,1,1,0,0];wdl=[0,0.2,0.35,0.65,0.8,1];k1=0:floor((M-1)/2);k2=floor((M-1)/2)+1:M-1;angH=[-alpha*(2*pi)/M*k1,alpha*(2*pi)/M*(M-k2)];H=Hrs.*exp(j*angH);h=real(ifft(H,M));[db,mag,pha,grd,w]=freqz_m(h,1);[Hr,ww,a,L]=Hr_Type2(h);subplot(2,2,1);stem(l,h);title('单位抽样响应');subplot(2,2,2);plot(w/pi,db);axis([0,1,-100,10]);title('幅频特性')set(gca,'XTickMode','manual','XTick',[0,0.2,0.35,0.65,0.8,1]);6.% x(n)=[1 1 1 1],长度为4的离散序列,求DTFT[x(n)];n=0:3;x=[1 1 1 1];k=0:1000;w=(2*pi/1000)*k;X=x*(exp(-j*2*pi/1000)).^(n'*k);magX=abs(X);angX=angle(X);realX=real(X);imagX=imag(X);hold onfigure(1);subplot(2,1,1);plot(k/500,magX);gridxlabel('frequency in pi units');title('Magnitude Part of DTFT') subplot(2,1,2);plot(k/500,angX/pi*180);gridxlabel('freqency in pi units');title('Anagle Part of DTFT')hold off%N=4,x=[1 1 1 1] 求DFT[x(n)]N=4; k=0:N-1X=dft(x,N);magX=abs(X)phaX=angle(X)*180/pihold onfigure(2);subplot(2,1,1);stem(k,magX);gridxlabel('k');ylabel('|X(k)|');title('Magmitude of the DFT:N=4') subplot(2,1,2);stem(k,phaX);gridxlabel('k');ylabel('Degrees');title('Angle of the DFT:N=4')hold off% N=8 x=[1 1 1 1 0 0 0 0],在后面补4个0再求DFT.x=[1,1,1,1,zeros(1,4)];N=8;k=0:N-1X=dft(x,N);magX=abs(X)phaX=angle(X)*180/pihold onfigure(3);subplot(2,1,1);stem(k,magX);gridxlabel('k');ylabel('|X(k)|');title('Magmitude of the DFT:N=8')subplot(2,1,2);stem(k,phaX);gridxlabel('k');ylabel('Degrees');title('Angle of the DFT:N=8')hold off% N=16 x=[1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0], 在后面补12个零后求DFT. x=[1,1,1,1,zeros(1,12)];N=16;k=0:N-1X=dft(x,N);magX=abs(X)phaX=angle(X)*180/pihold onfigure(4);subplot(2,1,1);stem(k,magX);gridxlabel('k');ylabel('|X(k)|');title('Magmitude of the DFT:N=16')subplot(2,1,2);stem(k,phaX);gridxlabel('k');ylabel('Degrees');title('Angle of the DFT:N=16')hold off7.A=[1,-0.9];B=[0.05,0.05]; %系统差分方程系数向量B和Ax1n=[1 1 1 1 1 1 1 1 zeros(1,60)]; %产生信号x1(n)=R8(n)x2n=ones(1,60); %产生信号x2(n)=u(n)hn=impz(B,A,40); %求系统单位脉冲响应h(n)subplot(3,1,1);stem(hn,'.'); %调用函数tstem绘图title('(a) 系统单位脉冲响应h(n)');box ony1n=filter(B,A,x1n); %求系统对x1(n)的响应y1(n)subplot(3,1,2);stem(y1n,'.');title('(b) 系统对R8(n)的响应y1(n)');box ony2n=filter(B,A,x2n); %求系统对x2(n)的响应y2(n)subplot(3,1,3);y='y2(n)';stem(y2n,'.');title('(c) 系统对u(n)的响应y2(n)');box on8.clear;n=[-20:20];b=[1,-1];a=[1,0.75,0.125];x1=(n==0);y1=filter(b,a,x1)附：前面用到的函数（1）function [b,a]=afd_butt(Wp,Ws,Rp,As);if Wp<=0error('Passband edge must be large than 0')endif Ws<=Wperror('stopband edge must be larger than Passband edge ')endif (Rp<=0)||(As<0)error('PB ripple and/or SB attenuation must be larger than 0 ') endN=ceil((log10((10^(Rp/10)-1)/(10^(As/10)-1)))/(2*log10(Wp/Ws))); fprintf('\n***Butterworth Filter Order =%2.0f\n',N)OmegaC=Wp/((10^(Rp/10)-1)^(1/(2*N)));[b,a]=u_buttap(N,OmegaC);（2）function [b0,B,A]=dir2cas(b,a);b0=b(1);b=b/b0;a0=a(1);a=a/a0;b0=b0/a0;M=length(b);N=length(a);if N>Mb=[b zeros(1,N-M)];else if M>Na=[a zeros(1,M-N)];N=M;elseNM=0;endK=floor(N/2);B=zeros(K,3);A=zeros(K,3);if K*2==N;b=[b 0];a=[a 0];endbroots=cplxpair(roots(b));aroots=cplxpair(roots(a));for i=1:2:2*KBrow=broots(i:1:i+1,:);Brow=real(poly(Brow));B(fix((i+1)/2),:)=Brow;Arow=aroots(i:1:i+1,:);Arow=real(poly(Arow));A(fix((i+1)/2),;)=Arow;end（3）function [db,mag,pha,grd,w] = freqz_m(b,a);% Modified version of freqz subroutine% ------------------------------------% [db,mag,pha,grd,w] = freqz_m(b,a);% db = Relative magnitude in dB computed over 0 to pi radians % mag = absolute magnitude computed over 0 to pi radians% pha = Phase response in radians over 0 to pi radians% grd = Group delay over 0 to pi radians% w = 501 frequency samples between 0 to pi radians% b = numerator polynomial of H(z) (for FIR: b=h)% a = denominator polynomial of H(z) (for FIR: a=[1])%[H,w] = freqz(b,a,1000,'whole');H = (H(1:1:501))'; w = (w(1:1:501))';mag = abs(H);db = 20*log10((mag+eps)/max(mag));pha = angle(H);% pha = unwrap(angle(H));grd = grpdelay(b,a,w);% grd = diff(pha);% grd = [grd(1) grd];% grd = [0 grd(1:1:500); grd; grd(2:1:501) 0];% grd = median(grd)*500/pi;（4）function [Hr,w,b,L] = Hr_Type2(h);% Computes Amplitude response of Type-2 LP FIR filter% ---------------------------------------------------% [Hr,w,b,L] = Hr_Type2(h)% Hr = Amplitude Response% w = frequencies between [0 pi] over which Hr is computed % b = Type-2 LP filter coefficients% L = Order of Hr% h = Type-2 LP impulse response%M = length(h);L = M/2;b = 2*[h(L:-1:1)];n = [1:1:L]; n = n-0.5;w = [0:1:500]'*pi/500;Hr = cos(w*n)*b';（5）function hd=ideal_lp(wc,M)alpha=(M-1)/2;n=[0:1:(M-1)];m=n-alpha+eps;hd=sin(wc*m)./(pi*m);（6）function[b,a]=imp_invr(c,d,T)[R,p,k]=residue(c,d);p=exp(p*T);[b,a]=residuez(R,p,k);b=real(b');a=real(a')（7）function [b,a]=u_buttap(N,OmegaC)[z,p,k]=buttap(N);p=p*OmegaC;k=k*OmegaC^N;B=real(poly(z));b=k*B;b0=k;a=real(poly(p));五．实验结果1.实验内容1的运行结果如下图所示：2.实验内容2的运行结果如下所示：a=1.0000 -3.3635 5.0684 -4.2759 2.1066 -0.5706 0.06613.实验内容3的运行结果如下所示：4.实验内容4的运行结果如下所示：Rp =0.0357As = 505.实验内容5的运行结果如下所示：6.实验内容6的运行结果如下所示：k =0 1 2 37.实验内容7的运行结果如下所示：8.实验内容8的运行结果如下所示：六．总结体会这次实验让我获益匪浅，它让我深刻体会到实验前的理论知识准备的重要性。