通信10-3李国庆 324.doc

数字信号处理实验指导书赵红怡姓名：李国庆学号：101010603242013年5月实验一 离散傅里叶变换（DFT）对确定信号进行谱分析一 实验目的1 加深对DFT算法原理和基本性质的理解。2 熟悉DFT算法和原理的编程方法。3 学习用DFT对信号进行谱分析的方法，了解可能出现的误差及其原因，以便在实际中正确利用。二 实验原理 一个连续信号的频谱可以用其傅里叶变换表示，即 若对进行理想采样可得采样序列 对进行DTFT，可得其频谱为：其中数字频率与模拟频率的关系为： 的DFT为 若是限带信号，且在满足采样定理的条件下，是的周期延拓，是在单位圆上的等间隔采样值，即。 为在计算机上分析计算方便，常用来近似，这样对于长度为N的有限长序列（无限长序列也可用有限长序列来逼近），便可通过DFT求其离散频谱。三 实验步骤1 复习DFT变换及其性质。2 编写DFT变换程序，分析DFT结果。 3按实验内容要求，上机实验，并写出实验报告。四 实验内容1用DFT对下列序列进行谱分析。 （1） （2）%a08.mN=100;n=0:1:N-1;xn=3*sin(0.04*pi*n);subplot(3,2,1);stem(n,xn);title(连续信号序列);Xk=dft(xn,N); %离散傅里叶变换magXk=abs(Xk) %求幅度谱phaXk=angle(Xk); %求相位谱k=0:N-1;subplot(3,2,3);stem(k,magXk);title(幅度谱);subplot(3,2,5);stem(k,phaXk);title(相位谱);xn1=1,1,1,1,0,0,0,0,0,0,0,0;N1=length(xn1);n1=0:1:N1-1;subplot(3,2,2);stem(n1,xn1); title(离散信号序列); %信号序列Xk1=dft(xn1,N1); %离散傅里叶变换magXk1=abs(Xk1) phaXk1=angle(Xk1);k1=0:N1-1;subplot(3,2,4);stem(k1,magXk1);title(幅度谱);subplot(3,2,6);stem(k1,phaXk1);title(相位谱);2为了说明高密度频谱和高分辨率频谱之间的区别，考察序列 （1）当0n10时，确定并画出x(n)的离散傅里叶变换。 （2）当0n100时，确定并画出x(n)的离散傅里叶变换。 %a09.mn=0:1:99;x=cos(0.48*pi*n)+cos(0.52*pi*n);n1=0:1:9;y1=x(1:1:10);subplot(3,3,1);stem(n1,y1);title(信号x(n),0=n=9);axis(0,10,-2.5,2.5);text(10.2,-2.5,n)Y1=fft(y1,10);magY1=abs(Y1(1:1:6);k1=0:1:5;w1=2*pi/10*k1;subplot(3,3,4);stem(w1/pi,magY1);title(FFT幅度);axis(0,1,0,10);text(1.05,0,w)xlabel(以pi为单位的频率)subplot(3,3,7);stem(w1/pi,angle(Y1(1:1:6),title(FFT相位谱);n2=0:1:99;y2=x(1:1:10),zeros(1,90);subplot(3,3,2);stem(n2,y2);title(信号x(n),0=n=9+90zeros);axis(0,100,-2.5,2.5);text(100.5,-2.5,n)Y2=fft(y2,100);magY2=abs(Y2(1:1:51);k2=0:1:50;w2=2*pi/100*k2;subplot(3,3,5);stem(w2/pi,magY2);title(FFT幅度);axis(0,1,0,10);text(1.05,0,w)xlabel(以pi为单位的频率)subplot(3,3,8);stem(w2/pi,angle(Y2(1:1:51),title(FFT相位谱);axis(0,1,-5,5);subplot(3,3,3);stem(n,x);title(信号x(n),0=n=99);axis(0,100,-2.5,2.5);text(100.5,-2.5,n)Y3=fft(x,100);magY3=abs(Y3(1:1:51);k3=0:1:50;w3=2*pi/100*k3;subplot(3,3,6);stem(w3/pi,magY3);title(FFT幅度);axis(0,1,0,60);text(1.05,0,w)xlabel(以pi为单位的频率)subplot(3,3,9);stem(w3/pi,angle(Y3(1:1:51),title(FFT相位谱);axis(0,1,-5,5);五 实验子程序1. DFT验证频谱泄漏特性程序(参见程序s11.m)2x(n)的10点离散傅里叶变换(参见程序s12.m)六 思考题1 若 在N=8和N=16两种情况下，和的幅频特性会相同吗？答：不相同。七 实验报告要求1 简述实验目的和实验原理。2 结合所给的程序，绘出实验序列的时域、频域图形。3 分析结果并与理论结果比较，总结实验所得结果。3 简要回答思考题。实验二 用FFT进行信号的谱分析一 实验目的1 加深对FFT算法原理和基本性质的理解。2 熟悉FFT算法和原理的编程方法。3 学习用FFT对信号进行谱分析的方法，了解可能出现的误差及其原因，以便在实际中正确利用。二 实验原理N点序列的N点DFT可以表示成利用系数的周期性和对称性 可采用FFT算法来计算的DFT。当时称为基2的FFT算法DIT-FFT。三 实验步骤1 复习FFT算法及其性质。2 编写FFT变换程序，分析FFT结果。 3按实验内容要求，上机实验，并写出实验报告。四 实验内容信号如下，对下面所给出的信号逐个进行谱分析。 %当N=8时的谱分析N=8;n=0:N-1;N1=16%对x(n)=R4(n)的谱分析xn=2,4,6,8;Xk=fft(xn,N);N2=length(xn);n2=0:N2-1;k=(0:length(abs(Xk)-1)*N/length(angle(Xk);subplot(3,4,1);stem(n2,xn);title(信号序列);subplot(3,4,5);stem(k,abs(Xk);title(x1(n)的幅度谱);subplot(3,4,9);stem(k,angle(Xk);title(x1(n)的相位谱);%对x(n)=cos(pi/4*n)的谱分析xn1=cos(pi/4*n);Xk1=fft(xn1,N);k1=(0:length(abs(Xk1)-1)*N/length(angle(Xk1);subplot(3,4,2);stem(n,xn1);title(信号序列);subplot(3,4,6);stem(k1,abs(Xk1);title(x2(n)的幅度谱);subplot(3,4,10);stem(k1,angle(Xk1);title(x2(n)的相位谱);%对x(n)=sin(8/pi*n)的谱分析xn2=sin(pi/8*n);Xk2=fft(xn2,N);k2=(0:length(abs(Xk2)-1)*N/length(angle(Xk2);subplot(3,4,3);stem(n,xn2);title(信号序列);subplot(3,4,7);stem(k2,abs(Xk2);title(x3(n)的幅度谱);subplot(3,4,11);stem(k2,angle(Xk2);title(x3(n)的相位谱);%对cos(8*pi*t)+cos(16*pi*t)+cos(20*pi*t)的谱分析fs=64;n1=0:N1-1;t=n1/fs;xn3=cos(8*pi*t)+cos(16*pi*t)+cos(20*pi*t);Xk3=fft(xn3,N1);k3=(0:length(abs(Xk3)-1)*N/length(abs(Xk3);subplot(3,4,4);stem(n1,xn3);title(信号序列);subplot(3,4,8);stem(k3,abs(Xk3);title(x4(n)的幅度谱);subplot(3,4,12);stem(k3,angle(Xk3);title(x4(n)的相位谱); %当N=16时的谱分析N=8;n=0:N-1;N1=16%对x(n)=R4(n)的谱分析xn=2,4,6,8;Xk=fft(xn,N1);N2=length(xn);n2=0:N2-1;k=(0:length(abs(Xk)-1)*N1/length(angle(Xk);subplot(3,4,1);stem(n2,xn);title(信号序列);subplot(3,4,5);stem(k,abs(Xk);title(x1(n)的幅度谱);subplot(3,4,9);stem(k,angle(Xk);title(x1(n)的相位谱);%对x(n)=cos(pi/4*n)的谱分析xn1=cos(pi/4*n);Xk1=fft(xn1,N1);k1=(0:length(abs(Xk1)-1)*N1/length(angle(Xk1);subplot(3,4,2);stem(n,xn1);title(信号序列);subplot(3,4,6);stem(k1,abs(Xk1);title(x2(n)的幅度谱);subplot(3,4,10);stem(k1,angle(Xk1);title(x2(n)的相位谱);%对x(n)=sin(8/pi*n)的谱分析xn2=sin(pi/8*n);Xk2=fft(xn2,N1);k2=(0:length(abs(Xk2)-1)*N1/length(angle(Xk2);subplot(3,4,3);stem(n,xn2);title(信号序列);subplot(3,4,7);stem(k2,abs(Xk2);title(x3(n)的幅度谱);subplot(3,4,11);stem(k2,angle(Xk2);title(x3(n)的相位谱);%对cos(8*pi*t)+cos(16*pi*t)+cos(20*pi*t)的谱分析fs=64;n1=0:N1-1;t=n1/fs;xn3=cos(8*pi*t)+cos(16*pi*t)+cos(20*pi*t);Xk3=fft(xn3,N1);k3=(0:length(abs(Xk3)-1)*N1/length(abs(Xk3);subplot(3,4,4);stem(n1,xn3);title(信号序列);subplot(3,4,8);stem(k3,abs(Xk3);title(x4(n)的幅度谱);subplot(3,4,12);stem(k3,angle(Xk3);title(x4(n)的相位谱);五 实验参考程序(参见程序s2.m)六 实验报告要求1 简述实验目的和实验原理。2 结合所给的程序，绘出实验序列的时域、频域图形。3 分析结果并与理论结果比较，总结实验所得结果。实验三 IIR数字滤波器的设计一 实验目的1加深对巴特沃思与切比雪夫数字滤波器基本原理的理解。2掌握双线性变换法设计巴特沃思与切比雪夫IIR数字滤波器的方法。3熟悉双线性变换法确定巴特沃思与切比雪夫IIR数字滤波器系统函数的方法。二 实验原理设计的基本思想是将模拟滤波器的频率响应压缩为数字滤波器的频率响应，模拟和数字的频率变换关系为，s平面到z平面的变换关系为，和的关系为：。设计的方法是根据技术指标，对先进行预畸，也就是计算：设计模拟滤波器使之满足技术指标，然后根据双线性变换公式将模拟滤波器转换成所希望的数字滤波器，即。三 实验步骤1复习双线性变换法设计IIR数字滤波器基本原理。2编写双线性变换法设计IIR数字滤波器程序。3按实验内容要求，上机实验，并写出实验报告。四 实验内容1 设计Butterworth数字低通滤波器 要求：wp=0.2*pi;wr=0.4*pi;Ap=3;Ar=20;T=1;Omegap=(2/T)*tan(wp/2);Omegar=(2/T)*tan(wr/2);ep=sqrt(10(Ap/10)-1);Ripple=sqrt(1/(1+ep*ep);Attn=1/(10(Ar/20);%双线性变换法设计Butterworth数字低通滤波器cs,ds=afd_butt(Omegap,Omegar,Ap,Ar)b,a=bilinear(cs,ds,T)C,B,A=dir2cas(b,a)db,mag,pha,w=freqz_m(b,a);subplot(2,3,1);plot(w/pi,mag);title(数字滤波器幅度响应|H(w)|); axis(0,1,0,1.1);subplot(2,3,2);plot(w/pi,db);title(数字滤波器衰减响应(dB);axis(0,1,-40,5);subplot(2,3,3);plot(w/pi,pha/pi);title(数字滤波器相位响应);b,a=afd_butt(wp,wr,Ap,Ar); C,B,A=sdir2cas(b,a);db,mag,pha,w=freqs_m(b,a,pi);subplot(2,3,4);plot(w/pi,mag);title(模拟滤波器幅度响应|H(Omega)|);subplot(2,3,5);plot(w/pi,db);title(模拟滤波器衰减响应(dB);axis(0,1,-40,5);subplot(2,3,6);plot(w/pi,pha/pi);title(模拟滤波器相位响应);2 设计Chebshev 1型数字低通滤波器(要求同3)wp=0.2*pi;wr=0.3*pi;Ap=1;Ar=15;T=1;Omegap=(2/T)*tan(wp/2);Omegar=(2/T)*tan(wr/2);ep=sqrt(10(Ap/10)-1);Ripple=sqrt(1/(1+ep*ep);Attn=1/(10(Ar/20);%双线性变换法设计Chebshev 1型数字低通滤波器cs,ds=afd_chb1(Omegap,Omegar,Ap,Ar);b,a=bilinear(cs,ds,T);C,B,A=dir2cas(b,a)db,mag,pha,w=freqz_m(b,a);subplot(2,3,1);plot(w/pi,mag);title(数字滤波器幅度响应|H(w)|);axis(0,1,0,1.1);subplot(2,3,2);plot(w/pi,db);title(数字滤波器衰减响应(dB);axis(0,1,-40,5);subplot(2,3,3);plot(w/pi,pha/pi);title(数字滤波器相位响应);%双线性变换法设计Chebshev 1型模拟低通滤波器b,a=afd_chb1(wp,wr,Ap,Ar)C,B,A=sdir2cas(b,a)db,mag,pha,w=freqs_m(b,a,pi);subplot(2,3,4);plot(w/pi,mag);title(模拟滤波器幅度响应|H(Omega)|);axis(0,1,0,1.1);subplot(2,3,5);plot(w/pi,db);title(数字滤波器衰减响应(dB);axis(0,1,-40,5);subplot(2,3,6);plot(w/pi,pha/pi);title(数字滤波器相位响应); 3设计Chebshev 2型数字低通滤波器 要求：wp=0.2*pi;wr=0.3*pi;Ap=1;Ar=15;T=1;Omegap=(2/T)*tan(wp/2);Omegar=(2/T)*tan(wr/2);ep=sqrt(10(Ap/10)-1);Ripple=sqrt(1/(1+ep*ep);Attn=1/(10(Ar/20);%双线性变换法设计Chebshev 2型数字低通滤波器cs,ds=afd_chb2(Omegap,Omegar,Ap,Ar);b,a=bilinear(cs,ds,T);C,B,A=dir2cas(b,a)db,mag,pha,w=freqz_m(b,a);subplot(2,3,1);plot(w/pi,mag);title(数字滤波器幅度响应|H(w)|);axis(0,1,0,1.1);subplot(2,3,2);plot(w/pi,db);title(数字滤波器衰减响应(dB);axis(0,1,-40,5);subplot(2,3,3);plot(w/pi,pha/pi);title(数字滤波器相位响应);%双线性变换法设计Chebshev 2型模拟低通滤波器b,a=afd_chb2(wp,wr,Ap,Ar)C,B,A=sdir2cas(b,a)db,mag,pha,w=freqs_m(b,a,pi);subplot(2,3,4);plot(w/pi,mag);title(模拟滤波器幅度响应|H(Omega)|);axis(0,1,0,1.1);subplot(2,3,5);plot(w/pi,db);title(数字滤波器衰减响应(dB);axis(0,1,-40,5);subplot(2,3,6);plot(w/pi,pha/pi);title(数字滤波器相位响应);实验结果五 参考程序1. 双线性变换法设计Butterworth数字低通滤波器(参见程序s31.m)2. 双线性变换法设计Chebshev 1型数字低通滤波器(参见程序s32.m)3. 双线性变换法设计设计Chebshev 2型数字低通滤波器(参见程序s33.m)六思考题用双线性变换法设计数字滤波器过程中，T的取值对设计结果有无影响？为什么？答：对结果无影响。因为在将S平面变换到Z平面时是s=c1-e-11+e-1到z=c+sc-s 的变换过程，与T无关。七实验报告要求 1简述实验目的和实验原理。2结合所给的程序，绘出实验图形。3分析结果并与理论结果比较，总结实验所得结果。4简要回答思考题。实验四 FIR数字滤波器的设计一 实验目的1 加深对FIR数字滤波器设计基本原理的理解。2 掌握窗函数法、频率取样法设计FIR数字滤波器的方法。3 熟悉各种窗函数对滤波器特性的影响。4 熟悉过渡带取样点对滤波器的影响，了解过渡点设置的方法。二 实验原理 若表示理想的低通滤波器，这个滤波器对应的单位脉冲响应为，关于对称。 为了得到一个因果的线性相位FIR滤波器，必须在两边同时截取得到长度为M的，即 根据窗函数的不同，可以设计具有不同过度带宽度和阻带最小衰减的FIR滤波器。FIR滤波器的系统函数为，对于长度为N的序列的频谱，其N个频域取样点可以不失真地代表它，也能够表达，即。第一种频率采样：第一种频率采样：三 实验步骤1复习窗函数法、频率取样法设计FIR数字滤波器的基本原理。2编写窗函数法、频率取样法设计FIR数字滤波器程序。3按实验内容要求，上机实验，并写出实验报告。四 实验内容1根据下列技术指标，窗函数法设计一个FIR低通滤波器： 求单位脉冲响应和频率响应，观察和带宽以及阻带最小衰减。wp=0.2*pi;wr=0.4*pi;tr_width=wr-wp %低通M=ceil(6.6*pi/tr_width)+2n=0:1:M-1;wc=(wr+wp)/2;hd=ideal_lp(wc,M);w_ham=(hamming(M);h=hd.*w_ham;db,mag,pha,w=freqz_m(h,1);delta_w=2*pi/1000;Ap=-(min(db(1:1:wp/delta_w+1)Ar=-round(max(db(wr/delta_w+1:1:501)subplot(2,2,1);stem(n,hd);title(理想脉冲响应)axis(0 M-1 -0.1 0.3);xlabel(n);ylabel(hd(n);subplot(2,2,2);stem(n,w_ham);title(海明窗)axis(0 M-1 0 1.1);xlabel(n);ylabel(w(n);subplot(2,2,3);stem(n,h);title(实际脉冲响应)axis(0 M-1 -0.1 0.3);xlabel(n);ylabel(h(n);subplot(2,2,4);plot(w/pi,db);title(幅度响应(单位:dB)axis(0 1 -100 10);xlabel(频率:pi);ylabel(分贝数);实验结果3 用矩形窗函数boxcar()设计低通、带通FIR数字滤波器，判断类型。（1）低通实验程序wp=0.2*pi;wr=0.4*pi;tr_width=wr-wp %低通M=ceil(1.8*pi/tr_width)+2n=0:1:M-1;wc=(wr+wp)/2;hd=ideal_lp(wc,M);w_ham=(boxcar(M); %矩形窗h=hd.*w_ham;db,mag,pha,w=freqz_m(h,1);delta_w=2*pi/1000;Ap=-(min(db(1:1:wp/delta_w+1)Ar=-round(max(db(wr/delta_w+1:1:501)subplot(2,2,1);stem(n,hd);title(理想脉冲响应)axis(0 M-1 -0.1 0.3);xlabel(n);ylabel(hd(n);subplot(2,2,2);stem(n,w_ham);title(海明窗)axis(0 M-1 0 1.1);xlabel(n);ylabel(w(n);subplot(2,2,3);stem(n,h);title(实际脉冲响应)axis(0 M-1 -0.1 0.3);xlabel(n);ylabel(h(n);subplot(2,2,4);plot(w/pi,db);title(幅度响应(单位:dB)axis(0 1 -100 10);xlabel(频率:pi);ylabel(分贝数);实验结果（2）带通实验程序wp1=0.2*pi;wr1=0.4*pi;wp2=0.6*pi;wr2=0.4*pi;tr_width=wr1-wp1 %带通M=ceil(1.8*pi/tr_width)+2n=0:1:M-1;wc1=(wr1+wp1)/2;wc2=(wr2+wp2)/2;hd=ideal_bdp(wc1,wc2,M); %带通w_ham=(boxcar(M); %矩形窗h=hd.*w_ham;db,mag,pha,w=freqz_m(h,1);delta_w=2*pi/1000;Ap=-(min(db(1:1:wp1/delta_w+1)Ar=-round(max(db(wr2/delta_w+1:1:501)subplot(2,2,1);stem(n,hd);title(理想脉冲响应)axis(0 M-1 -0.1 0.3);xlabel(n);ylabel(hd(n);subplot(2,2,2);stem(n,w_ham);title(矩形窗)axis(0 M-1 0 1.1);xlabel(n);ylabel(w(n);subplot(2,2,3);stem(n,h);title(实际脉冲响应)axis(0 M-1 -0.1 0.3);xlabel(n);ylabel(h(n);subplot(2,2,4);plot(w/pi,db);title(幅度响应(单位:dB)axis(0 1 -100 10);xlabel(频率:pi);ylabel(分贝数);ideal_bdp(wc1,wc2,M)函数程序为：function hd=ideal_bdp(wc1,wc2,M);alpha=(M-1)/2;n=0:1:(M-1);m=n-alpha+eps;hd=sin(wc1*pi*m)-sin(wc2*pi*m)./(pi*m);实验结果4 用三角窗函数triang()和汉宁窗函数hanning()设计低通FIR数字滤波器。（选作）5 根据下列技术指标，频率采样法1设计一个FIR低通滤波器：1. 频率采样法1程序%a25.mM=21;alpha=(M-1)/2;l=0:M-1;wl=(2*pi/M)*l;Hrs=1,1,1,1,zeros(1,14),1,1,1;Hdr=1,1,0,0;wdl=0,0.35,0.35,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,w=freqz_m(h,1);Hr,ww,a,L=hr_type1(h);subplot(2,2,1);plot(wl(1:11)/pi,Hrs(1:11),.,wdl,Hdr);axis(0 1 -0.1 1.1);title(频率样本H（k）：M=21)subplot(2,2,2);stem(l,h);axis(-1 M -0.1 0.4);title(脉冲响应h(n);xlabel();subplot(2,2,3);plot(ww/pi,Hr,wl(1:11)/p