数字信号处理实验程序2

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-p2/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-p2/q,p=8,q=4'q=8;x=exp(-(n-p.2/q;xk3=fft(x,16;subplot(3,1,3;stem(n,x;title('exp(-(n-p2/q,

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

3、'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-p2/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-p2/q,p=14,q=8'%时域特性figure(4;subplot(3,1,1;stem(n,abs(xk1;title('exp(-(n-p2/q,p=8,q=8'subplot(3,1,2;

4、stem(n,abs(xk2;title('exp(-(n-p2/q,p=13,q=8'subplot(3,1,3;stem(n,abs(xk3;title('exp(-(n-p2/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.4375;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

5、,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(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(

6、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<=3xc(m+1=m;xd(m+1=4-m;elseif(4<=m & m<=7xc(m+1=8-m;xd(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

7、(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=fft(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.4cl

8、cclose 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;title('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

9、;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/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

10、(-(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(4,1,1;stem(0:N-1,xa;title('xa(n=exp(-(n-p2/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(nxb(n'xak1=fft(xa,2*N;xbk1=fft(x

11、b,2*N;h=ifft(xak1.*xbk1,2*N;subplot(4,1,4;stem(0:2*N-2,h(1:2*N-1;title('xa*xb'2.6clcclose allN=512;xe=rand(1,N-0.5;for m=0:7;if(0<=m & m<=3xc(m+1=m;elseif(4<=m & m<=7xc(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:

12、142;fm(72:142=ifft(xmk.*xck,71;if m=2f(1:7=fm(1:7;endif m>1f(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;endif m=8f(m-1*64+8:m*64+7=fm(79:142;endendstem(0:512+6,f;%重叠相加法title('重叠相加法'for m=1:8xm1(8:71=xe(m-1*64+1:m*64;if m=1xm1(1:7=0;else xm1(1:7=xe(m-1*64-6:(m-1*64;en

13、dxmk1=fft(xm1,71;fm1=ifft(xmk1.*xck,71;f1(m-1*64+1:m*64=fm1(8:71;endxmk2=fft(xe(506:512,71;h=ifft(xmk2.*xck,71;f1(513:519=h(8:14;figure(2stem(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/q;xb=exp(-a*

14、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*tan(2*pi*200/(2*1000; N,wn=che

15、b1ord(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*pi*300;n1,wn1 = buttord(wc,wr,1,25,'s

16、' 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=bilinear(B,A,1000; h2,w=freqz(num2,den2;f=w/pi*

17、500;plot(f,20*log10(abs(h1,'-.',f,20*log10(abs(h2 ,'-'legend('脉冲响应不变法','双线性变换法' axis(0,500,-80,10grid;xlabel('频率/Hz'ylabel('幅度/dB'3.3clearclose allwc=2*8000*tan(2*pi*1200/(2*8000;wr=2*8000*tan(2*pi*2000/(2*8000;n,wn=buttord(wc,wr,0.5,40,'s'B,A=

18、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,0.5,40,'s'B,A=ellip(n,0.5,40,wn,'low','s'num,d

19、en=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'legend('Butterworth','Chebyshev','Ellicptic' hold onx=1200

20、;plot(x3.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,den;f=w/pi*15000;plot(f,20*log10(abs(h1,'-.',f,20*log10(abs(ab s(h2;grid on;axis(0,3

21、000,-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*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=cheb

22、y1(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%N=45,计算并画出矩形框、汉明窗、布莱克曼窗的归一化的幅度谱,并比较各自的主要特点%(1矩形窗(Rectangle Window调用格式:w=boxcar(n,根据长度n 产生一个矩形窗w。%(2三角窗(Triangular

23、Window调用格式:w=triang(n,根据长度n 产生一个三角窗w。%(3汉宁窗(Hanning Window调用格式:w=hanning(n,根据长度n 产生一个汉宁窗w。%(4海明窗(Hamming Window调用格式:w=hamming(n,根据长度n 产生一个海明窗w。%(5布拉克曼窗(Blackman Window调用格式:w=blackman(n,根据长度n 产生一个布拉克曼窗w。%(6恺撒窗(Kaiser Window调用格式:w=kaiser(n,beta,根据长度n 和影响窗函数旁瓣的参数产生一个恺撒窗wclear allclose allN=45;w1=boxcar

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

25、figure(3h,w=freqz(w3,N;plot(w/pi,20*log10(abs(h;axis(0,1,-150,10;grid onxlabel('归一化频率/'ylabel('幅度/dB'title('布莱克曼窗'4.2%N=15,带通滤波器的两个通带边界分别是w1=0.3,w2=0.5。用汉宁窗设计此线性 相位滤波器，观察 %它的实际 3dB 和 20dB 带宽。N=45，重复 这一设计，观察幅频和相位特性的变化，注 意 N 变化的影响。 close all clear all N=15; w1=0.3; w2=0.5; w=ha

26、nning(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,-80,10; xlabel('归一化频率/' ylabel('幅度/dB' subplot(2,1,2; plot(w3/pi,angle(h1; grid on; axis(0,1,-4,4; xlabel('归一化频率/' ylabel('角度/rad'

27、 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 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('

28、角度/rad' 4.3 close all clear all N=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,w3=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=freq

29、z(hn2,1; figure(2 plot(w3/pi,20*log10(abs(h1; grid on; axis(0,1,-80,10; xlabel('归一化频率/' ylabel('幅度/dB' title('布莱克曼窗,N=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(

30、0,1,-80,10; xlabel('归一化频率/' ylabel('幅度/dB' title('矩形窗,N=45' h1,w3=freqz(hn2,1; figure(4 plot(w3/pi,20*log10(abs(h1; grid on; axis(0,1,-110,10; xlabel('归一化频率/' ylabel('幅度/dB' title('布莱克曼窗,N=45' 4.4 close all clear all N=40; %beta=4 for n=1:3 if n=1 beta

31、=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/pi,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

32、title('beta=10' end; subplot(2,1,2; plot(w1/pi,angle(h1; grid on; axis(0,1,-4,4; xlabel('归一化频率/' ylabel('角度/rad' end 4.5 clear all close all N=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

33、|(w>=1.6*pi && w<=1.8*pi hk(1,k+1=1; end end k=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(hn,1; plot(w1/pi,20*log10(abs(h1; grid on; axis(0,1,-80,10; xlabel('归一化频率/' ylabel