数字信号处理实验二,三

1、数字信号处理上机实验二实验报告班级1302013 学号姓名：贾昊禛1. 给定模拟信号，对其进行采样，用DFT（FFT）进行信号频谱分析。（） 确定最小采样频率和最小采样点数。（） 若以秒进行采样，至少需要取多少采样点？（） 用DFT的点数画出信号的点DFT的幅度谱，讨论幅度谱结果。（） 分别为和，能否分辨出信号的所有频率分量。（） 在（）和（）的条件下做补0 FFT，分析结果。（） 在不满足最小采样点数的情况下做补0DFT，观察是否可以分辨出两个频率分量。(1)w1=4*pi;w2=8*pi;f1=w1/2/pi;f2=w2/2/pi;disp('最小采样频率

2、：')fs1=2*max(f1,f2)f=f2-f1;disp('最小采样点数：')N= ceil(fs1/abs(f);运行结果：最小采样频率：fs1 = 8最小采样点数：N = 4(2) T=0.01;fs2=1/T;disp('以T=0.01n采样，最少的采样点个数为：')N0=ceil(fs2/abs(f)运行结果：以T=0.01n采样，最少的采样点个数为：N0 =50（3）（4）N1=50;N2=100;N3=64;N4=60;n1=0:N1-1;n2=0:N2-1;n3=0:N3-1;n4=0:N4-1;x1=2*cos(w1*n1*T)+5

3、*cos(w2*n1*T);x2=2*cos(w1*n2*T)+5*cos(w2*n2*T);x3=2*cos(w1*n3*T)+5*cos(w2*n3*T);x4=2*cos(w1*n4*T)+5*cos(w2*n4*T);X1=abs(fft(x1,N1);X2=abs(fft(x2,N2);X3=abs(fft(x3,N3);X4=abs(fft(x4,N4);figure(1)subplot(2,2,1),stem(n1,X1,'.');title('N=50');xlabel('n1');ylabel('X1');sub

4、plot(2,2,2),stem(n2,X2,'.');title('N=100');xlabel('n2');ylabel('X2');subplot(2,2,3),stem(n3,X3,'.');title('N=64');xlabel('n3');ylabel('X3');subplot(2,2,4),stem(n4,X4,'.');title('N=60');xlabel('n4');ylabel('X4

5、');运行结果：（5）Nb=200;nb=0:Nb-1;X5=abs(fft(x1,Nb);X6=abs(fft(x2,Nb);X7=abs(fft(x3,Nb);X8=abs(fft(x4,Nb);figure(2)subplot(2,2,1),stem(nb,X5,'.');title('N=50补零到200后的幅度值');xlabel('nb');ylabel('X5');subplot(2,2,2),stem(nb,X6,'.');title('N=100补零到200后的幅度值');

6、xlabel('nb');ylabel('X6');subplot(2,2,3),stem(nb,X7,'.');title('N=64补零到200后的幅度值');xlabel('nb');ylabel('X7');subplot(2,2,4),stem(nb,X8,'.');title('N=60补零到200后的幅度值');xlabel('nb');ylabel('X8');运行结果：（6）N9=3;n9=0:N9-1;x9=2*co

7、s(w1*n9*T)+5*cos(w2*n9*T);X9=abs(fft(x9,Nb);figure(3)stem(nb,X9,'.');title('N=3不满足最小采样点时补零到200后的幅度值');xlabel('nb');ylabel('X9');2. 设雷达发射线性调频信号，采样率，采样点数。回波信号，。（） 画出的频谱。（） 利用DFT的时延性质产生，比较直接在时域产生和在频域产生（再变换到时域）的结果是否相同。（） 匹配滤波的结果是，（“”表示线性卷积）。分别用直接线性卷积和DFT的卷积定理求解。比较二者结果，并记录

8、两种方法的运行时间（用tic，toc指令）。（） 画出的频谱。（1）mu=5e13;fs=2e9;Ts=1/fs;N=20000;tao1=1e-6;tao2=1.1e-6;delay1=ceil(tao1/Ts);delay2=ceil(tao2/Ts);n=0:N-1;t=n*Ts;h=exp(j*2*pi*mu*t.2);H=abs(fft(h,N);figure(1);stem(n,H,'.');运行结果：（2）Ns=N+max(delay1,delay2);s1=zeros(1,Ns);s2=zeros(1,Ns);s1(delay1+1:delay1+N)=exp(

9、2*j*pi*mu*t.2);s2(delay2+1:delay2+N)=exp(2*j*pi*mu*t.2);s=s1+s2;f=(0:1/Ns:1-1/Ns)*fs;x=zeros(1,Ns);x(1:N)=h;x1=fft(x).*exp(-j*2*pi*tao1*f);x2=fft(x).*exp(-j*2*pi*tao2*f);xsum=ifft(x1+x2);figure(2)plot(abs(s)hold onplot(abs(xsum),'r')（3）hh=zeros(1,N);hh(1)=h(1);hh(2:end)=fliplr(h(2:end);hh=co

10、nj(hh);ticy=conv(s,hh);toctaob=(0:length(y)-1)*Ts-N*Ts;figure(3);plot(taob,abs(y)hold onNy=N+Ns-1;tichhf=conj(fft(h,Ny); hf=fft(s,Ny);yf=(ifft(hf.*hhf);tocstem(taob,circshift(abs(yf.'),N),'r')xlabel('10');ylabel('a');运行结果：Elapsed time is 3.462685 seconds.Elapsed time is 0

11、.133051 seconds.（4）figure(4)plot(abs(fft(yf)实验结果：数字信号处理实验三：数字滤波器设计及其实现1. IIR滤波器设计（1） 用matlab确定一个数字IIR低通滤波器所有四种类型的最低阶数。指标如下：40kHz的采样率，4kHz的通带边界频率，8kHz的阻带边界频率，0.5dB的通带波纹，40dB的最小阻带衰减。并在同一张图中画出每种滤波器的频率响应。（1）程序：fc=40;fp=4;fs=8;rp=0.5;rs=40;wp=2*pi*fp/fc;ws=2*pi*fs/fc;disp('buttord')n,wc=buttord(w

12、p,ws,rp,rs,'s')b,a=butter(n,wc,'low','s');w=0:0.002:5;h,w=freqs(b,a,w);h=20*log10(abs(h);plot(w,h,'b-')disp('cheb1ord')n,wpo=cheb1ord(wp,ws,rp,rs,'s')b,a=cheby1(n,rp,wpo,'low','s');h,w=freqs(b,a,w);h=20*log10(abs(h);hold onplot(w,h,'

13、;r-')disp('cheb2ord')n,wso=cheb2ord(wp,ws,rp,rs,'s')b,a=cheby2(n,rs,wso,'low','s');h,w=freqs(b,a,w);h=20*log10(abs(h);hold onplot(w,h,'k-')disp('ellipord')n,wc=ellipord(wp,ws,rp,rs,'s')b,a=ellip(n,rp,rs,wc,'low','s');h,w=fre

14、qs(b,a,w);h=20*log10(abs(h);hold onplot(w,h,'g-')legend('butter','cheby1','cheby2','ellip')grid onxlabel('w')ylabel('h')运行结果：buttordn = 9wc = 0.7533cheb1ordn = 5wpo = 0.6283cheb2ordn = 5wso = 1.2069ellipordn = 4wc = 0.6283（2） 用matlab确定一个数字IIR高通

15、滤波器所有四种类型的最低阶数。指标如下：3500Hz的采样率，1050Hz的通带边界频率，600Hz的阻带边界频率，1dB的通带波纹，50dB的最小阻带衰减。并在同一张图中画出每种滤波器的频率响应。程序：fc=3500;fp=1050;fs=600;rp=1;rs=50;wp=2*pi*fp/fc;ws=2*pi*fs/fc;disp('buttord')n,wc=buttord(wp,ws,rp,rs,'s')b,a=butter(n,wc,'high','s');w=0:0.001:6;h,w=freqs(b,a,w);h=2

16、0*log10(abs(h);plot(w,h,'b-')disp('cheb1ord')n,wpo=cheb1ord(wp,ws,rp,rs,'s')b,a=cheby1(n,rp,wpo,'high','s');h,w=freqs(b,a,w);h=20*log10(abs(h);hold onplot(w,h,'r-')disp('cheb2ord')n,wso=cheb2ord(wp,ws,rp,rs,'s')b,a=cheby2(n,rs,wso,'

17、high','s');h,w=freqs(b,a,w);h=20*log10(abs(h);hold onplot(w,h,'k-')disp('ellipord')n,wc=ellipord(wp,ws,rp,rs,'s')b,a=ellip(n,rp,rs,wc,'high','s');h,w=freqs(b,a,w);h=20*log10(abs(h);hold onplot(w,h,'g-')legend('butter','cheby1

18、9;,'cheby2','ellip')axis(0,5,-120,0)grid onxlabel('w')ylabel('h')运行结果：buttordn = 12wc = 1.7402cheb1ordn = 7wpo = 1.8850cheb2ordn = 7wso = 1.2049ellipordn = 5wc = 1.8850（3） 用matlab确定一个数字IIR带通滤波器所有四种类型的最低阶数。指标如下：7kHz的采样率，1.4kHz和2.1kHz的通带边界频率，1.05kHz和2.45kHz的阻带边界频率，0.4dB

19、的通带波纹，50dB的最小阻带衰减。并在同一张图中画出每种滤波器的频率响应。程序：fc=7;fp=1.4,2.1;fs=1.05,2.45;rp=0.4;rs=50;wp=2*pi*fp/fc;ws=2*pi*fs/fc;disp('buttord')n,wc=buttord(wp,ws,rp,rs,'s')b,a=butter(n,wc ,'s');w=0:0.002:5;h,w=freqs(b,a,w);h=20*log10(abs(h);plot(w,h,'b-')disp('cheb1ord')n,wpo=

20、cheb1ord(wp,ws,rp,rs,'s')b,a=cheby1(n,rp,wpo ,'s');h,w=freqs(b,a,w);h=20*log10(abs(h);hold onplot(w,h,'r-')disp('cheb2ord')n,wso=cheb2ord(wp,ws,rp,rs,'s')b,a=cheby2(n,rs,wso ,'s');h,w=freqs(b,a,w);h=20*log10(abs(h);hold onplot(w,h,'k-')disp(

21、9;ellipord')n,wc=ellipord(wp,ws,rp,rs,'s')b,a=ellip(n,rp,rs,wc ,'s');h,w=freqs(b,a,w);h=20*log10(abs(h);hold onplot(w,h,'g-')legend('butter','cheby1','cheby2','ellip')axis(0,5,-120,20)grid onxlabel('w')ylabel('h')buttord运行结果

22、：n = 12wc = 1.2305 1.9250cheb1ordn = 7wpo = 1.2566 1.8850cheb2ordn = 7wso = 1.1050 2.1437ellipordn = 5wc = 1.2566 1.8850（4） 用matlab确定一个数字IIR带阻滤波器所有四种类型的最低阶数。指标如下：12kHz的采样率，2.1kHz和4.5kHz的通带边界频率，2.7kHz和3.9kHz的阻带边界频率，0.6dB的通带波纹，45dB的最小阻带衰减。并在同一张图中画出每种滤波器的频率响应。程序：fc=12;fp=2.1,4.5;fs=2.7,3.9;rp=0.5;rs=40

23、;wp=2*pi*fp/fc;ws=2*pi*fs/fc;disp('buttord')n,wc=buttord(wp,ws,rp,rs,'s')b,a=butter(n,wc,'stop','s');w=0:0.002:5;h,w=freqs(b,a,w);h=20*log10(abs(h);plot(w,h,'b-')disp('cheb1ord')n,wpo=cheb1ord(wp,ws,rp,rs,'s')b,a=cheby1(n,rp,wpo,'stop',

24、's');h,w=freqs(b,a,w);h=20*log10(abs(h);hold onplot(w,h,'r-')disp('cheb2ord')n,wso=cheb2ord(wp,ws,rp,rs,'s')b,a=cheby2(n,rs,wso,'stop','s');h,w=freqs(b,a,w);h=20*log10(abs(h);hold onplot(w,h,'k-')disp('ellipord')n,wc=ellipord(wp,ws,rp,r

25、s,'s')b,a=ellip(n,rp,rs,wc,'stop','s');h,w=freqs(b,a,w);h=20*log10(abs(h);hold onplot(w,h,'g-')legend('butter','cheby1','cheby2','ellip',4)grid onaxis(0,5,-120,20)xlabel('w')ylabel('h')运行结果：buttordn = 10wc = 1.2726 2.268

26、4cheb1ordn = 6wpo = 1.2252 2.3561cheb2ordn = 6wso = 1.3845 2.0851ellipordn = 5wc =1.996 2.3562用到的函数： butter，buttord，cheb2ord，chebl1，cheby2，ellip，ellipord.2. FIR滤波器设计分别用矩形窗，Blackman窗，Hamming窗，Hanning窗和Bartlett窗设计截止频率为0.3Pi，窗长为M（M=11，41，81，121）的FIR低通滤波器。在图中画出：（1） 理想低通滤波器的冲激响应；程序：wc=0.3*pi; M=11,41,81,

27、121; for i=1:1:4tao=(M(i)-1)/2; h_d1n=zeros(1,N(i); for n=1:1:N(i) if n=tao h_d1n(n)=wc/pi; else h_d1n(n)=(sin(wc*(n-tao)/(pi*(n-tao); endendsubplot(2,2,i)n1=1:1:M(i);stem(n1,h_d1n,'b');title( 'a = ', num2str(i) );end运行结果：（2） 所加窗函数；程序：M=11,41,81,121; for i=1:1:4n=0:1:M(i)-1;wn_boxcar

28、=boxcar(M(i); wn_bartlett=bartlett(M(i); wn_hanning=hanning(M(i); wn_hamming=hamming(M(i); wn_blackman=blackman(M(i); figure(i)subplot(3,2,1),stem(n,wn_boxcar,'r.');title('M=', num2str(i) ,'矩形窗函数'); subplot(3,2,2),stem(n,wn_bartlett,'r.');title('M=', num2str(i

29、) ,' 三角窗函数'); subplot(3,2,3),stem(n,wn_hanning,'r.');title('M=', num2str(i) ,' 汉宁窗函数'); subplot(3,2,4),stem(n,wn_hamming,'r.');title('M=', num2str(i) ,'海明窗函数'); subplot(3,2,5),stem(n,wn_blackman,'r.');title('M=', num2str(i) ,

30、9; 布莱克曼窗函数'); end运行结果：（3） 加窗后的滤波器冲激响应（4） 滤波器的幅频特性；程序：clear all; clc; wc=0.3; N=11,41,81,121; k=1;for i=1:1:4M=N(i)-1;h_boxcar=fir1(M,wc,boxcar(N(i);h_bartlett=fir1(M,wc,bartlett(N(i);h_hanning=fir1(M,wc,hanning(N(i);h_hamming=fir1(M,wc,hamming(N(i);h_blackman=fir1(M,wc,blackman(N(i);n=0:1:N(i)-1

31、; figure(k)k=k+1;subplot(3,2,1),stem(n,h_boxcar,'.');title('N=',num2str(N(i),' 矩形窗函数冲击响应'); subplot(3,2,2),stem(n,h_bartlett,'.');title('N=',num2str(N(i), '三角窗函数冲击响应');subplot(3,2,3),stem(n,h_hanning,'.');title('N=',num2str(N(i), '汉

32、宁窗函数冲击响应'); subplot(3,2,4),stem(n,h_hamming,'.');title('N=',num2str(N(i), '海明窗函数冲击响应'); subplot(3,2,5),stem(n,h_blackman,'.');title('N=',num2str(N(i), '布莱克曼窗函数冲击响应'); Hx_boxcar=freqz(h_boxcar);dBHx_boxcar=20*log10(abs(Hx_boxcar)/max(abs(Hx_boxcar);

33、 Hx_bartlett=freqz(h_bartlett);dBHx_bartlett=20*log10(abs(Hx_bartlett)/max(abs( Hx_bartlett);Hx_hanning=freqz(h_hanning);dBHx_hanning=20*log10(abs(Hx_hanning)/max(abs(Hx_hanning); Hx_hamming=freqz(h_hamming);dBHx_hamming=20*log10(abs(Hx_hamming)/max(abs(Hx_hamming);Hx_blackman=freqz(h_blackman);dBHx

34、_blackman=20*log10(abs(Hx_blackman)/max(abs( Hx_blackman); w=0:pi/512:511*pi/512; figure(k)k=k+1;subplot(3,2,1),plot(w/pi,dBHx_boxcar);title('N=',num2str(N(i),'矩形窗幅度响应 (dB)');xlabel('omega/pi');ylabel('20log|H(ejomega)|(dB)'); grid on;subplot(3,2,2),plot(w/pi,dBHx_bar

35、tlett); title('N=',num2str(N(i),'三角窗幅度响应 (dB)');xlabel('omega/pi');ylabel('20log|H(ejomega)|(dB)'); grid on; subplot(3,2,3),plot(w/pi,dBHx_hanning); title('N=',num2str(N(i),'汉宁窗幅度响应 (dB)');xlabel('omega/pi');ylabel('20log|H(ejomega)|(dB)');grid on; subplo