数字信号课程设计.doc

实验一：用DFT作谱分析一、实验目的1、进一步加深对DFT的基本性质的理解。2、熟悉FFT算法原理和FFT子程序的应用。3、学习用DFT对连续信号和离散信号进行谱分析的方法，了解可能出现的分析误差及其原因，以便在实际中正确应用DFT。二、上级程序：23x1=1 1 1 1;x2=1 2 3 4 4 3 2 1;x3=4,3,2,1,1,2,3,4;figuresubplot(1,3,1);stem(x1);title(序列x1)xlabel(n);ylabel(x1);subplot(1,3,2);stem(x2);title(序列x2)xlabel(n);ylabel(x2);subplot(1,3,3);stem(x3);title(序列x3)xlabel(n);ylabel(x3);N1=8;F1x1=fft(x1,N1);F1x2=fft(x2,N1);F1x3=fft(x3,N1);figuresubplot(1,3,1);stem(abs(F1x1);title(序列x1的8点谱分析)xlabel(n);ylabel(F1x1);subplot(1,3,2);stem(abs(F1x2);title(序列x2的8点谱分析)xlabel(n);ylabel(F1x2);subplot(1,3,3);stem(abs(F1x3);title(序列x3的8点谱分析)xlabel(n);ylabel(F1x3);N2=16;F2x1=fft(x1,N2);F2x2=fft(x2,N2);F2x3=fft(x3,N2);figuresubplot(1,3,1);stem(abs(F2x1);title(序列x1的16点谱分析);xlabel(n);ylabel(F2x1);subplot(1,3,2);stem(abs(F2x2);title(序列x2的16点谱分析);xlabel(n);ylabel(F2x2);subplot(1,3,3);stem(abs(F2x3);title(序列x3的16点谱分析);xlabel(n);ylabel(F2x3);N3=16;n=0:N3-1;x4=cos(n*pi/4);x5=sin(n*pi/8);figure;subplot(1,2,1);stem(n,x4);title(序列x4)xlabel(n);ylabel(x4);subplot(1,2,2);stem(n,x5);title(序列x5)xlabel(n);ylabel(x5);N4=8;F3x4=fft(x4,N4);F3x5=fft(x5,N4);figuresubplot(1,2,1);stem(abs(F3x4);title(序列x4的8点谱分析);xlabel(n);ylabel(F3x4);subplot(1,2,2);stem(abs(F3x5);title(序列x5的8点谱分析);xlabel(n);ylabel(F3x5);N5=16;F4x4=fft(x4,N5);F4x5=fft(x5,N5);figuresubplot(1,2,1);stem(abs(F4x4);title(序列x4的16点谱分析);xlabel(n);ylabel(F4x4);subplot(1,2,2);stem(abs(F4x5);title(序列x5的16点谱分析);xlabel(n);ylabel(F4x5);N6=64;F5x4=fft(x4,N6);F5x5=fft(x5,N6);figurew=2/N6*0:N6-1;subplot(1,2,1);plot(w,abs(F5x4);title(序列x4的64点谱分析);xlabel(w);ylabel(F5x4);subplot(1,2,2);plot(w,abs(F5x5);title(序列x5的64点谱分析);xlabel(w);ylabel(F5x5);N7=256;F6x4=fft(x4,N7);F6x5=fft(x5,N7);figurew=2/N7*0:N7-1;subplot(1,2,1);plot(w,abs(F6x4);title(序列x4的256点谱分析);xlabel(w);ylabel(F6x4);subplot(1,2,2);plot(w,abs(F6x5);title(序列x5的256点谱分析);xlabel(w);ylabel(F6x5);N8=64;t=0:1/64:1/64*(N8-1);x6=cos(8*pi*t)+cos(16*pi*t)+cos(20*pi*t);figureplot(t,x6);title(序列x6)xlabel(t);ylabel(x6);N9=8;F9x6=fft(x6,N9);figurestem(abs(F9x6);title(序列x6的8点谱分析);xlabel(n);ylabel(F9x6);N10=16;F10x6=fft(x6,N10);figurestem(abs(F10x6);title(序列x6的16点谱分析);xlabel(n);ylabel(F10x6);N11=64;F11x6=fft(x6,N11);w=2/N11*0:N11-1;figureplot(w,abs(F11x6);title(序列x6的64点谱分析);xlabel(w);ylabel(F11x6);N12=256;F12x6=fft(x6,N12);w=2/N12*0:N12-1;figureplot(w,abs(F12x6);title(序列x6的256点谱分析);xlabel(w);ylabel(F12x6);三、思考题1、在N=8时，x2(n)和x3(n)的幅频特性相同吗？为什么？N=16呢？答：在N=8时相同，因为x2(n)和x3(n)是圆周移位的关系，两者包含的信息完全相同，故幅频特性相同；在N=16时不相同，因为x2(n)和x3(n)并非圆周移位的关系，两者所包含的信息不等，故幅频特性相同。2、x4(n)和x5(n)在N=8、16、64和256是的频谱是否相同？为什么？答：对x4(n)来说，频谱相同，因为x4(n)的周期是8，在N=8、16、64和256时，取样点数都是周期的整数倍，包含信号的全部信息，故频谱相同；对x5(n)来说，因为x5(n)的周期是16，在N=8时，对信号进行了截断，包含的信息不完整，而在N=16、64和256时抽样点包含信号的全部信息，故频谱相同；3、x6(n)在N=8、16、64和256时能否将三条谱线分开，为什么？答：在N=8、16时不能将三条谱线分开，在N=64和256时能将三条谱线分开。因为x6(n)的最高频率是f=10,由奈奎斯特抽样定理知抽样频率fs2f=20时，才能无失真的恢复原来的信号，即可以将原信号频谱分开。实验二、用双线性法设计IIR数字滤波器一、 实验目的1、 熟悉双线性法设计IIR数字滤波器的原理和方法。2、 掌握数字滤波器的计算机仿真方法3、 通过观察对实际心电图信号的滤波作用，获得数字滤波器的感性认识。二、 上机程序T=1;Wp=2/T*tan(0.2*pi/2);Ws=2/T*tan(0.3*pi/2);Rp=1;Rs=15;N,Wc=buttord(Wp,Ws,Rp,Rs,s);B,A=butter(N,Wc,s); C,D=bilinear(B,A,1/T);W=0:0.001*pi:0.5*pi;H=freqs(B,A,W);Hd=freqz(C,D,W);figureplot(W/pi,abs(H)grid ontitle(模拟巴特沃斯滤波器)xlabel(Frequency/Hz)ylabel(Magnitude)figure;plot(W/pi,abs(Hd)grid ontitle(数字巴特沃斯滤波器)xlabel(Didital Frequency/pi)ylabel(Magnitude)Hd_db=-20*log10(abs(Hd(1)./(abs(Hd)+eps);figure;plot(W/pi,Hd_db)grid ontitle(数字巴特沃斯滤波器波特图)xlabel(Didital Frequency/pi)ylabel(bd_Magnitude)x=-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,-38,60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,6,4,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0;M=512;Fx=fft(x,M);y=filter(C,D,x);Fy=fft(y,M);figure;subplot(1,2,1);plot(x);title(心电图序列);xlabel(n);ylabel(x);f=1/M*0:M-1*250;subplot(1,2,2);plot(f,abs(Fx);title(原序列的频谱分析);xlabel(f);ylabel(|Fx|)figure;subplot(1,2,1);plot(y);title(滤波后的信号);xlabel(n);ylabel(y);f=1/M*0:M-1*250;subplot(1,2,2);plot(f,abs(Fy);title(滤波后的频谱分析);xlabel(f);ylabel(|Fy|);三、思考题1、本实验中用双线性法设计数字滤波器的映射公式s=(2/T)(1-z-1)/(1+z-1)中的T对设计的模拟滤波器波形有无影响？对设计的数字滤波器波形有无影响？为什么？答：对设计的模拟滤波器波形有影响，对设计的数字滤波器波形无影响；因为由=T可知，当数字指标指定时，改变T对有影响，但对数字指标无影响，也即对设计的模拟滤波器波形有影响，对设计的数字滤波器波形无影响。2、滤波前后心电图的波形有何不同？为什么？答：滤波后信号的高于50Hz的分量都被滤掉了。因为滤波器的阻带开始频率fst=37.5Hz,所以信号的高于50Hz的分量都被滤掉了。实验三、用窗函数法设计FIR数字滤波器一、实验目的1、 熟悉用窗函数法设计FIR数字滤波器的原理和方法。2、 熟悉线性FIR数字滤波器特性。3、 了解各种窗函数对滤波器性能的影响。二、上机程序 function hd=ideal(N,wc)for n=0:N-1 if n=(N-1)/2 hd(n+1)=wc/pi; else hd(n+1)=sin(wc*(n-(N-1)/2)/(pi*(n-(N-1)/2); endendclear all;cla;close all;N=33;wc=pi/4;hd=ideal(N,wc);w1=boxcar(N);w2=hamming(N);w3=hann(N);w4=blackman(N);h1=hd.*w1;h2=hd.*w2;h3=hd.*w3;h4=hd.*w4;M=512;fh1=fft(h1,M);db1=-20*log10(abs(fh1(1)./(abs(fh1)+eps);fh2=fft(h2,M);db2=-20*log10(abs(fh2(1)./(abs(fh2)+eps);fh3=fft(h2,M);db3=-20*log10(abs(fh3(1)./(abs(fh3)+eps);fh4=fft(h4,M);db4=-20*log10(abs(fh4(1)./(abs(fh4)+eps);w=2/M*0:M-1;figuresubplot(2,2,1);stem(h1)subplot(2,2,2);plot(w,abs(fh1)subplot(2,2,3);plot(w,db1)subplot(2,2,4);plot(w,angle(fh1)figuresubplot(2,2,1);stem(h2)subplot(2,2,2);plot(w,abs(fh2)subplot(2,2,3);plot(w,db2)subplot(2,2,4);plot(w,angle(fh2)figuresubplot(2,2,1);stem(h3)subplot(2,2,2);plot(w,abs(fh3)subplot(2,2,3);plot(w,db3)subplot(2,2,4);plot(w,angle(fh3)figuresubplot(2,2,1);stem(h4)subplot(2,2,2);plot(w,abs(fh4)subplot(2,2,3);plot(w,db4)subplot(2,2,4);plot(w,angle(fh4)N1=15;wc=pi/4;hd1=ideal(N1,wc);w1=hamming(N1);hn1=hd1.*w1;M=512;fh1=fft(hn1,M);w=2/M*0:M-1;figuresubplot(2,2,1);stem(hn1);title(单位抽样响应h1(n);xlabel(n);ylabel(h1(n);subplot(2,2,2);plot(w,abs(fh1);title(频率幅度);xlabel(w);ylabel(|H(k)|);db1=-20*log10(abs(fh1(1)./(abs(fh1)+eps);subplot(2,2,3);plot(w,db1);title(db频率);xlabel(w);ylabel(db);subplot(2,2,4);plot(w,angle(fh1);title(相位特性);xlabel(w);ylabel(angle(fh1);N2=33;wc=pi/4;hd2=ideal(N2,wc);w2=hamming(N2);hn2=hd2.*w2;M=512;fh2=fft(hn2,M);w=2/M*0:M-1;figuresubplot(2,2,1);stem(hn2);title(单位抽样响应h2(n);xlabel(n);ylabel(h2(n);subplot(2,2,2);plot(w,abs(fh2);title(频率幅度);xlabel(w);ylabel(|H2(k)|);db2=-20*log10(abs(fh2(1)./(abs(fh2)+eps);subplot(2,2,3);plot(w,db2);title(db频率);xlabel(w);ylabel(db2);subplot(2,2,4);plot(w,angle(fh2);title(相位特性);xlabel(w);ylabel(angle(fh2);N=512;t=0:1/512:1/N*(N-1);x1=sin(100*pi*t)+sin(200*pi*t)+sin(300*pi*t);figure;plot(t,x1);title(原信号x(t)波形);xlabel(t);ylabel(x(t);figure;fh1=fft(x1,N);f=1/N*0:N-1*512;plot(f,abs(fh1);title(滤波前的信号频谱)xlabel(f);ylabel(|fh1|);N=40;Wp=100*pi/512;Wst=150*pi/512