数字信号处理MATLAB实验程序

1、2.1clcclose all;n=0:15;p=8;q=2;x=exp(-(n-p).2/q);figure(1);subplot(3,1,1);stem(n,x);title('exp(-(n-p)2/q),p=8,q=2');xk1=fft(x,16);q=4;x=exp(-(n-p).2/q);subplot(3,1,2);xk2=fft(x,16);stem(n,x);title('exp(-(n-p)2/q),p=8,q=4');q=8;x=exp(-(n-p).2/q);xk3=fft(x,16);subplot(3,1,3);stem(n,x);

2、title('exp(-(n-p)2/q),p=8,q=8');%时域特性figure(2);subplot(3,1,1);stem(n,abs(xk1);title('exp(-(n-p)2/q),p=8,q=2');subplot(3,1,2);stem(n,abs(xk2);title('exp(-(n-p)2/q),p=8,q=4');subplot(3,1,3);stem(n,abs(xk3);title('exp(-(n-p)2/q),p=8,q=8');%频域特性%p=8;q=8;figure(3);subplot(

3、3,1,1);stem(n,x);title('exp(-(n-p)2/q),p=8,q=8');xk1=fft(x,16);p=13;x=exp(-(n-p).2/q);subplot(3,1,2);xk2=fft(x,16);stem(n,x);title('exp(-(n-p)2/q),p=13,q=8');p=14;x=exp(-(n-p).2/q);xk3=fft(x,16);subplot(3,1,3);stem(n,x);title('exp(-(n-p)2/q),p=14,q=8');%时域特性figure(4);subplot(

4、3,1,1);stem(n,abs(xk1);title('exp(-(n-p)2/q),p=8,q=8');subplot(3,1,2);stem(n,abs(xk2);title('exp(-(n-p)2/q),p=13,q=8');subplot(3,1,3);stem(n,abs(xk3);title('exp(-(n-p)2/q),p=14,q=8');%频域特性2.2clcclose alln=0:15;f=0.0625;a=0.1;xb1=exp(-a*n).*sin(2*pi*f*n);xk1=fft(xb1,16);f=0.43

5、75;xb2=exp(-a*n).*sin(2*pi*f*n);xk2=fft(xb2,16);f=0.5625;xb3=exp(-a*n).*sin(2*pi*f*n);xk3=fft(xb3,16);figure(1);subplot(3,1,1);stem(n,xb1);title('f=0.0625');subplot(3,1,2);stem(n,xb2);title('f=0.4375');subplot(3,1,3);stem(n,xb3);title('f=0.5625');figure(2);subplot(3,1,1);stem

6、(n,abs(xk1);title('f=0.0625');subplot(3,1,2);stem(n,abs(xk2);title('f=0.4375');subplot(3,1,3);stem(n,abs(xk3);title('f=0.5625');2.3clcclose all;N=8;n=0:N-1;xc=0:7;xd=0:7;for m=0:7; if(0<=m & m<=3) xc(m+1)=m; xd(m+1)=4-m; elseif(4<=m & m<=7) xc(m+1)=8-m; xd

7、(m+1)=m-4; endendfigure(1);subplot(2,1,1);stem(n,xc);title('三角波序列');subplot(2,1,2);stem(n,xd);title('反三角波序列');xck=fft(xc,N);xdk=fft(xd,N);figure(2);subplot(2,1,1);stem(n,abs(xck);title('三角波序列8点fft');subplot(2,1,2);stem(n,abs(xdk);title('反三角波序列8点fft');figure(3);xck1=ff

8、t(xc,32);xdk1=fft(xd,32);subplot(2,1,1);stem(0:4*N-1,abs(xck1);title('三角波序列32点fft');subplot(2,1,2);stem(0:4*N-1,abs(xdk1);title('反三角波序列32点fft');2.4clcclose allN=16;detaf=1/16;n=0:N-1;x=sin(2*pi*0.125*n)+cos(2*pi*(0.125+detaf)*n);xk=fft(x,N);figure(1);subplot(2,1,1);stem(n,abs(xk);tit

9、le('N=16,f=1/16');detaf=1/64;x=sin(2*pi*0.125*n)+cos(2*pi*(0.125+detaf)*n);xk1=fft(x,N);subplot(2,1,2);stem(n,abs(xk1);title('N=16,f=1/64');N=128;detaf=1/16;n=0:N-1;x=sin(2*pi*0.125*n)+cos(2*pi*(0.125+detaf)*n);xk=fft(x,N);figure(2);subplot(2,1,1);stem(n,abs(xk);title('N=128,f=1/

10、16');detaf=1/64;x=sin(2*pi*0.125*n)+cos(2*pi*(0.125+detaf)*n);xk1=fft(x,N);subplot(2,1,2);stem(n,abs(xk1);title('N=128,f=1/64');2.5clcclose allN=16;p=8;q=2;a=.1;f=.0625;n=0:N-1;xa=exp(-(n-p).2/q);xb=exp(-a*n).*sin(2*pi*f*n);xak=fft(xa,N);xbk=fft(xb,N);f=ifft(xak.*xbk,N);figure(1);subplot

11、(4,1,1);stem(0:N-1,xa);title('xa(n)=exp(-(n-p)2/q)');subplot(4,1,2);stem(0:N-1,xb);title('xb(n)=exp(-a*n)*sin(2*pi*f*n)');subplot(4,1,3);stem(0:N-1,f);title('xa(n)xb(n)');xak1=fft(xa,2*N);xbk1=fft(xb,2*N);h=ifft(xak1.*xbk1,2*N);subplot(4,1,4);stem(0:2*N-2,h(1:2*N-1);title(

12、9;xa*xb');2.6clcclose allN=512;xe=rand(1,N)-0.5;for m=0:7; if(0<=m & m<=3) xc(m+1)=m; elseif(4<=m & m<=7) xc(m+1)=8-m; endendxck=fft(xc,71);f=0:518;fm=0:141;for m=1:8; xm=xe(m-1)*64+1:m*64); xmk=fft(xm,71); fm(1:71)=fm(72:142); fm(72:142)=ifft(xmk.*xck,71); if m=2 f(1:7)=fm(1

13、:7); end if m>1 f(m-2)*64+8:(m-1)*64)=fm(8:64); f(m-1)*64+1:(m-1)*64+7)=fm(65:71)+fm(72:78); end if m=8 f(m-1)*64+8:m*64+7)=fm(79:142); end end stem(0:512+6,f);%重叠相加法title('重叠相加法');for m=1:8 xm1(8:71)=xe(m-1)*64+1:m*64); if m=1 xm1(1:7)=0; else xm1(1:7)=xe(m-1)*64-6:(m-1)*64); end xmk1=ff

14、t(xm1,71); fm1=ifft(xmk1.*xck,71); f1(m-1)*64+1:m*64)=fm1(8:71);end xmk2=fft(xe(506:512),71);h=ifft(xmk2.*xck,71);f1(513:519)=h(8:14);figure(2)stem(0:518,f1);title('重叠保留法');fa=conv(xc,xe);detaf=f1-f;figure(3);stem(0:518,detaf);2.7clcclose allN=16;p=8;q=2;a=.1;n=0:N-1;f=0.0625;xa=exp(-(n-p).2

15、/q);xb=exp(-a*n).*sin(2*pi*f*n);xak=fft(xa,N*2);xbk=fft(xb,2*N);rm=real(ifft(conj(xak).*xbk);rm=rm(N+2:2*N) rm(1:N);m=(-N+1):N-1;stem(m,rm);title('线性相关');figure(2);xak1=fft(xa,N);xbk1=fft(xb,N);rm1=real(ifft(xak1.*xbk1,N);stem(n,rm1);3.1clearclose allwc=2*1000*tan(2*pi*300/(2*1000);wt=2*1000

16、*tan(2*pi*200/(2*1000);N,wn=cheb1ord(wc,wt,0.8,20,'s');B,A=cheby1(N,0.8,wn,'high','s');num,den=bilinear(B,A,1000);h,w=freqz(num,den);f=w/pi*500;plot(f,20*log10(abs(h);axis(0,500,-80,10);grid;xlabel('频率/Hz');ylabel('幅度/dB');3.2clear;close all;wc=2*pi*200;wr=2*p

17、i*300;n1,wn1 = buttord(wc,wr,1,25,'s');B,A=butter(n1,wn1,'low','s');num1,den1=impinvar(B,A,1000);h1,w=freqz(num1,den1);wc=2*1000*tan(2*pi*200/(2*1000);wr=2*1000*tan(2*pi*300/(2*1000);n,wn = buttord(wc,wr,1,25,'s');B,A=butter(n,wn,'low','s');num2,den2=b

18、ilinear(B,A,1000);h2,w=freqz(num2,den2);f=w/pi*500;plot(f,20*log10(abs(h1),'-.',f,20*log10(abs(h2),'-');legend('脉冲响应不变法','双线性变换法');axis(0,500,-80,10)grid;xlabel('频率/Hz');ylabel('幅度/dB');3.3clearclose allwc=2*8000*tan(2*pi*1200/(2*8000);wr=2*8000*tan(2*

19、pi*2000/(2*8000);n,wn=buttord(wc,wr,0.5,40,'s');B,A=butter(n,wn,'low','s');num,den=bilinear(B,A,8000);h1,w=freqz(num,den);n,wn=cheb1ord(wc,wr,0.5,40,'s');B,A=cheby1(n,0.5,wn,'low','s');num,den=bilinear(B,A,8000);h2,w=freqz(num,den);n,wn=ellipord(wc,wr,

20、0.5,40,'s');B,A=ellip(n,0.5,40,wn,'low','s');num,den=bilinear(B,A,8000);h3,w=freqz(num,den);f=w/pi*4000;plot(f,20*log10(abs(h1),'-.',f,20*log10(abs(h2),'.-',f,20*log10(abs(h3),'-');axis(0,4000,-100,10);grid;xlabel('频率/Hz');ylabel('幅度/dB'

21、;);legend('Butterworth','Chebyshev','Ellicptic')hold onx=1200;plot(x)3.4clearclose all;n,wc=buttord(2000 3000,1500 6000,3,20,'s');B,A=butter(n,wc,'bandpass','s');num,den=bilinear(B,A,30000);h1,w=freqz(num,den);num,den=impinvar(B,A,30000);h2,w=freqz(num

22、,den);f=w/pi*15000;plot(f,20*log10(abs(h1),'-.',f,20*log10(abs(abs(h2);grid on;axis(0,3000,-100,0);legend('脉冲响应不变法','双线性变换');xlabel('频率/Hz');ylabel('幅度/dB');3.5clearclose allw1=2*10000*tan(2*pi*1000/(2*10000);w2=2*10000*tan(2*pi*2000/(2*10000);w3=2*10000*tan(2

23、*pi*500/(2*10000);w4=2*10000*tan(2*pi*3000/(2*10000);n,wn=cheb1ord(w3 w4,w1 w2,3,18,'s')B,A=cheby1(n,3,wn,'stop','s');num,den=bilinear(B,A,10000);h,w=freqz(num,den);f=w/pi*5000;plot(f,20*log10(abs(h);axis(0,5000,-100,0);gridxlabel('频率/Hz');ylabel('幅度/dB');4.1%

24、N=45,计算并画出矩形框、汉明窗、布莱克曼窗的归一化的幅度谱，并比较各自的主要特点%（1）矩形窗（Rectangle Window） 调用格式：w=boxcar(n)，根据长度 n 产生一个矩形窗 w。%（2）三角窗（Triangular Window） 调用格式：w=triang(n)，根据长度 n 产生一个三角窗 w。%（3）汉宁窗（Hanning Window） 调用格式：w=hanning(n)，根据长度 n 产生一个汉宁窗 w。%（4）海明窗（Hamming Window） 调用格式：w=hamming(n)，根据长度 n 产生一个海明窗 w。%（5）布拉克曼窗（Blackman

25、Window）调用格式：w=blackman(n)，根据长度 n 产生一个布拉克曼窗 w。%（6）恺撒窗（Kaiser Window） 调用格式：w=kaiser(n,beta)，根据长度 n 和影响窗函数旁瓣的参数产生一个恺撒窗wclear allclose allN=45;w1=boxcar(N);w2=hamming(N);w3=blackman(N);h,w=freqz(w1,N);figure(1)plot(w/pi,20*log10(abs(h);axis(0,1,-80,10);grid onxlabel('归一化频率/');ylabel('幅度/dB&#

26、39;);title('矩形窗');figure(2)h,w=freqz(w2,N);plot(w/pi,20*log10(abs(h);axis(0,1,-80,10);grid onxlabel('归一化频率/');ylabel('幅度/dB');title('汉明窗');figure(3)h,w=freqz(w3,N);plot(w/pi,20*log10(abs(h);axis(0,1,-150,10);grid onxlabel('归一化频率/');ylabel('幅度/dB');titl

27、e('布莱克曼窗');4.2%N=15,带通滤波器的两个通带边界分别是w1=0.3,w2=0.5。用汉宁窗设计此线性相位滤波器，观察%它的实际3dB和20dB带宽。N=45，重复这一设计，观察幅频和相位特性的变化，注意N变化的影响。close allclear allN=15;w1=0.3;w2=0.5;w=hanning(N);n=0:N-1;alfa=(N-1)/2;h=fir1(N-1,w1 w2,w);h1,w3=freqz(h,1);figure(1)subplot(2,1,1);plot(w3/pi,20*log10(abs(h1);grid on;axis(0,1

28、,-80,10);xlabel('归一化频率/');ylabel('幅度/dB');subplot(2,1,2);plot(w3/pi,angle(h1);grid on;axis(0,1,-4,4);xlabel('归一化频率/');ylabel('角度/rad');N=45;w=hanning(N);n=0:N-1;alfa=(N-1)/2;h=fir1(N-1,w1 w2,w);h1,w3=freqz(h,1);figure(2)subplot(2,1,1);plot(w3/pi,20*log10(abs(h1);grid

29、on;axis(0,1,-80,10);xlabel('归一化频率/');ylabel('幅度/dB');subplot(2,1,2);plot(w3/pi,angle(h1);grid on;axis(0,1,-4,4);xlabel('归一化频率/');ylabel('角度/rad');4.3close allclear allN=15;w1=0.3;w2=0.5;wn1=boxcar(N);wn2=blackman(N);hn1=fir1(N-1,w1 w2,wn1);hn2=fir1(N-1,w1 w2,wn2);h1,w

30、3=freqz(hn1,1);figure(1)plot(w3/pi,20*log10(abs(h1);grid on;axis(0,1,-80,10);xlabel('归一化频率/');ylabel('幅度/dB');title('矩形窗,N=15');h1,w3=freqz(hn2,1);figure(2)plot(w3/pi,20*log10(abs(h1);grid on;axis(0,1,-80,10);xlabel('归一化频率/');ylabel('幅度/dB');title('布莱克曼窗,N

31、=15');N=45;wn1=boxcar(N);wn2=blackman(N);hn1=fir1(N-1,w1 w2,wn1);hn2=fir1(N-1,w1 w2,wn2);h1,w3=freqz(hn1,1);figure(3)plot(w3/pi,20*log10(abs(h1);grid on;axis(0,1,-80,10);xlabel('归一化频率/');ylabel('幅度/dB');title('矩形窗,N=45');h1,w3=freqz(hn2,1);figure(4)plot(w3/pi,20*log10(abs

32、(h1);grid on;axis(0,1,-110,10);xlabel('归一化频率/');ylabel('幅度/dB');title('布莱克曼窗,N=45');4.4close allclear allN=40;%beta=4for n=1:3 if n=1 beta=4; elseif n=2 beta=6; else beta=10; end;w=kaiser(N,beta);h=fir1(N-1,0.2 0.4 0.6 0.8,w);h1,w1=freqz(h,1);figure(n)subplot(2,1,1);plot(w1/p

33、i,20*log10(abs(h1);grid on;axis(0,1,-80,10);xlabel('归一化频率/');ylabel('幅度/dB');if n=1 title('beta=4'); elseif n=2 title('beta=6'); else title('beta=10'); end;subplot(2,1,2);plot(w1/pi,angle(h1);grid on;axis(0,1,-4,4);xlabel('归一化频率/');ylabel('角度/rad&#

34、39;);end4.5clear allclose allN=45;k=0:N-1;for k=0:N-1 w=2*pi/N*k; hk(1,k+1)=0; if (w>=0.2*pi) && (w<=0.4*pi)|(w>=0.6*pi && w<=0.8*pi)|(w>=1.2*pi && w<=1.4*pi)|(w>=1.6*pi && w<=1.8*pi) hk(1,k+1)=1; endendk=0:N-1;hk(1,5)=0.5;hk(1,11)=0.5;hk(1,14)=0.5;hk(1,20)=0.5;hk(1,27)=0.5;hk(1,33)=0.5;hk(1,36)=0.5;hk(1,42)=0.5;thetak=-k*2*pi/N*(N-1)/2);hk1=hk.*exp(j*thetak);hn=ifft(hk1);h1,w1=freqz(h