西电数字信号处理上机实验报告

1、数字信号处理上机实验报张吉凯第一次上机实验一：设给定模拟信号，的单位是ms。(1) 利用MATLAB绘制出其时域波形和频谱图（傅里叶变换），估计其等效带宽（忽略谱分量降低到峰值的3%以下的频谱）。(2) 用两个不同的采样频率对给定的进行采样。比较两种采样率下的信号频谱，并解释。（1）MATLAB程序：N=10; Fs=5;Ts=1/Fs;n=-N:Ts:N; xn=exp(-abs(n); w=-4*pi:0.01:4*pi;X=xn*exp(-j*(n'*w);subplot(211)plot(n,xn);title('x_a(t)时域波形'

2、);xlabel('t/ms');ylabel('x_a(t)');axis(-10, 10, 0, 1);subplot(212);plot(w/pi,abs(X);title('x_a(t)频谱图');xlabel('omega/pi');ylabel('X_a(e(jomega)');ind = find(X >=0.03*max(X)*0.01;eband = (max(ind) -min(ind);fprintf('等效带宽为 %fKHZn',eband);运行结果：等效带宽为 12

3、.110000KHZ（2）MATLAB程序：N=10; omega=-3*pi:0.01:3*pi; %Fs=5000Fs=5;Ts=1/Fs;n=-N:Ts:N; xn=exp(-abs(n); X=xn*exp(-j*(n'*omega);subplot(2,2,1);stem(n,xn);grid on;axis(-10, 10, 0, 1.25);title('时域波形(f_s=5000)');xlabel('n');ylabel('x_1(n)');subplot(2,2,2);plot(omega/pi,abs(X);titl

4、e('频谱图(f_s=5000)');xlabel('omega/pi');ylabel('X_1(f)');grid on;%Fs=1000Fs=1;Ts=1/Fs;n=-N:Ts:N; xn=exp(-abs(n); X=xn*exp(-j*(n'*omega);subplot(2,2,3);stem(n,xn);grid on;axis(-10, 10, 0, 1.25);title('时域波形(f_s=1000)');xlabel('n');ylabel('x_2(n)');grid

5、 on;subplot(2,2,4);plot(omega/pi,abs(X);title('频谱图(f_s=1000)');xlabel('omega/pi');ylabel('X_2(f)');grid on;运行结果：实验二：给定一指数型衰减信号，采样率，为采样周期。为方便起见，重写成复指数形式。采样后的信号为，加窗后长度为的形式为：这3个信号，的幅度谱平方分别为：模拟信号：采样信号：加窗（取有限个采样点）信号：且满足如下关系：实验内容（1） 在同一张图上画出：模型号幅度谱平方；（2） 在同一张图上画出：模型号幅度谱平方；改变值，结果又如何

6、？（1）MATLAB程序：f=0:0.01:3;alpha=0.2;f0=0.5;L=10;T1=1;T2=0.5;Xa=1./(alpha2+(2*pi*(f-f0).2);Xs1=T1*(1-2*exp(-alpha*T1*L)*cos(2*pi*(f-f0)*T1*L)+exp(-2*alpha*T1*L)./(1-2*exp(-alpha*T1)*cos(2*pi*(f-f0)*T1)+exp(-2*alpha*T1);Xs2=T2*(1-2*exp(-alpha*T2*L)*cos(2*pi*(f-f0)*T2*L)+exp(-2*alpha*T2*L)./(1-2*exp(-alp

7、ha*T2)*cos(2*pi*(f-f0)*T2)+exp(-2*alpha*T2);plot(f,Xa,'b');hold on;plot(f,Xs1,'g');hold on;plot(f,Xs2,'r');xlabel('f/Hz');ylabel('|X(f)|2');grid on;legend('模拟信号幅度谱平方|X(f)|2', 'f_s=1Hz时，采样信号幅度谱平方|TX(f)|2', 'f_s=2Hz时，采样信号幅度谱平方|TX(f)|2');运

8、行结果：（2）MATLAB程序：f=0:0.01:3;alpha=0.2;f0=0.5;L1=5;L2=10;L3=20;T1 = 0.5Xa=1./(alpha2+(2*pi*(f-f0).2);Xs1=T1*(1-2*exp(-alpha*T1*L1)*cos(2*pi*(f-f0)*T1*L1)+exp(-2*alpha*T1*L1)./(1-2*exp(-alpha*T1)*cos(2*pi*(f-f0)*T1)+exp(-2*alpha*T1);Xs2=T1*(1-2*exp(-alpha*T1*L2)*cos(2*pi*(f-f0)*T1*L2)+exp(-2*alpha*T1*L

9、2)./(1-2*exp(-alpha*T1)*cos(2*pi*(f-f0)*T1)+exp(-2*alpha*T1);Xs3=T1*(1-2*exp(-alpha*T1*L3)*cos(2*pi*(f-f0)*T1*L3)+exp(-2*alpha*T1*L3)./(1-2*exp(-alpha*T1)*cos(2*pi*(f-f0)*T1)+exp(-2*alpha*T1);plot(f,Xa,'b');hold on;plot(f,Xs1,'g');hold on;plot(f,Xs2,'r');hold on;plot(f,Xs3,&#

10、39;y')xlabel('f/Hz');ylabel('|X(f)|2');grid on;legend('模拟信号幅度谱平方|X(f)|2', 'f_s=2Hz时，采样信号幅度谱平方|TX(f)|2(L=5)','f_s=2Hz时，采样信号幅度谱平方|TX(f)|2(L=10)','f_s=2Hz时，采样信号幅度谱平方|TX(f)|2(L=20)');运行结果：实验三：设，编写MATLAB程序，计算：（） 点圆周卷积；（） 点圆周卷积；（） 线性卷积；（） 画出的，和时间轴对齐。MATL

11、AB程序：a = 1,2,2;b = 1,2,3,4;y1 = cconv(a,b,5);y2 = cconv(a,b,6);y3 = conv(a,b);figure(1);subplot(311)stem(y1);grid ontitle('五点圆周卷积y1(n)');xlabel('n'),ylabel('y1(n)');axis(0 6 0 15)subplot(312)stem(y2);grid ontitle('六点圆周卷积y2(n)');xlabel('n'),ylabel('y2(n)

12、9;);axis(0 6 0 15)subplot(313)stem(y3);grid ontitle('线性卷积y3(n)');xlabel('n'),ylabel('y3(n)');axis(0 6 0 15);运行结果：x1=1,2,2;x2=1,2,3,4;n1=0:4;y1=cconv(x1,x2,5);n2=0:5;y2=cconv(x1,x2,6);n3=0:length(x1)+length(x2)-2;y3=conv(x1,x2);subplot(3,1,1);stem(n1,y1);grid on;axis(-1,6,0,16

13、);subplot(3,1,2);stem(n2,y2);grid on;axis(-1,6,0,16);subplot(3,1,3);stem(n3,y3);grid on;axis(-1,6,0,16);运行结果：实验四：给定因果系统：（） 求系统函数并画出零极点示意图。（） 画出系统的幅频特性和相频特性。（） 求脉冲响应并画序列图。提示：在中，zplane(b,a) 函数可画零极点图；Freqz(b,a,N)可给出范围内均匀间隔的点频率响应的复振幅；Impz(b,a,N)可求的逆变换（即脉冲响应）。MATLAB程序：a = 1,0b = 1,-0.9figure(1)zplane(b,a

14、);title('零极点分布图')w=-3*pi:0.01:3*pi;h,phi=freqz(b,a,w);figure(2);subplot(3,1,1);plot(w, abs(h);grid on;title('幅频特性');xlabel('f/Hz'),ylabel('H(f)');subplot(3,1,2);plot(w, phi);grid on;title('相频特性');xlabel('f/Hz'),ylabel('W(f)');subplot(3,1,3);imp

15、z(b,a);运行结果：第二次上机1. 给定模拟信号，对其进行采样，用DFT（FFT）进行信号频谱分析。（） 确定最小采样频率和最小采样点数。（） 若以秒进行采样，至少需要取多少采样点？（） 用DFT的点数画出信号的点DFT的幅度谱，讨论幅度谱结果。（） 分别为和，能否分辨出信号的所有频率分量。（） 在（）和（）的条件下做补0 FFT，分析结果。（） 在不满足最小采样点数的情况下做补0DFT，观察是否可以分辨出两个频率分量。（1）最小采样频率：8 ；最小采样点数：4（2） 最小采样点数：50（3）（4）MATLAB程序：N1=50;N2=100;N3=64;N4=60;n1=0:N1-1;n2

16、=0:N2-1;n3=0:N3-1;n4=0:N4-1;w1=4*pi;w2=8*pi;T=0.01;x1=2*cos(w1*n1*T)+5*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,'.'

17、;);grid on;title('N=50');xlabel('n1');ylabel('X1');subplot(2,2,2),stem(n2,X2,'.');grid on;title('N=100');xlabel('n2');ylabel('X2');subplot(2,2,3),stem(n3,X3,'.');grid on;title('N=64');xlabel('n3');ylabel('X3');su

18、bplot(2,2,4),stem(n4,X4,'.');grid on;title('N=60');xlabel('n4');ylabel('X4');运行结果：（5）MATLAB程序：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,'.');grid on;title('N=50补零到200后的幅度值

19、9;);xlabel('nb');ylabel('X5');subplot(2,2,2),stem(nb,X6,'.');grid on;title('N=100补零到200后的幅度值');xlabel('nb');ylabel('X6');subplot(2,2,3),stem(nb,X7,'.');grid on;title('N=64补零到200后的幅度值');xlabel('nb');ylabel('X7');subplot(2

20、,2,4),stem(nb,X8,'.');grid on;title('N=60补零到200后的幅度值');xlabel('nb');ylabel('X8');运行结果：（6）MATLAB程序：N9=3;n9=0:N9-1;x9=2*cos(w1*n9*T)+5*cos(w2*n9*T);X9=abs(fft(x9,Nb);figure(3)stem(nb,X9,'.'); grid on;title('N=3不满足最小采样点时补零到200后的幅度值');xlabel('nb');

21、ylabel('X9');运行结果：2. 设雷达发射线性调频信号，采样率，采样点数。回波信号，。（） 画出的频谱。（） 利用DFT的时延性质产生，比较直接在时域产生和在频域产生（再变换到时域）的结果是否相同。（） 匹配滤波的结果是，（“”表示线性卷积）。分别用直接线性卷积和DFT的卷积定理求解。比较二者结果，并记录两种方法的运行时间（用tic，toc指令）。（） 画出的频谱。（1）MATLAB程序：mu=5e13;fs=2e9;Ts=1/fs;mu=5e13;fs=2e9;Ts=1/fs;N=20000;tao1=1e-6;tao2=1.1e-6;delay1=ceil(tao

22、1/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,'.');grid on;运行结果：（2）MATLAB程序：Ns=N+max(delay1,delay2);s1=zeros(1,Ns);s2=zeros(1,Ns);s1(delay1+1:delay1+N)=exp(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;

23、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);grid on;hold on;plot(abs(xsum),'r');grid on;运行结果：（3）MATLAB程序：hh=zeros(1,N);hh(1)=h(1);hh(2:end)=fliplr(h(2:end);hh=conj(hh);tic;y=conv(s,hh);toc;taob=(0:length(y)-1)*Ts-

24、N*Ts;figure(3);plot(taob,abs(y);grid on;hold on;Ny=N+Ns-1;tic;hhf=conj(fft(h,Ny); hf=fft(s,Ny);yf=(ifft(hf.*hhf);toc;stem(taob,circshift(abs(yf.'),N),'r');grid on;xlabel('10');ylabel('a');运行结果：Elapsed time is 3.462685 seconds.Elapsed time is 0.133051 seconds.（4）MATLAB程序：f

25、igure(4);plot(abs(fft(yf);grid on;实验结果：第三次上机1. IIR滤波器设计（1）用matlab确定一个数字IIR低通滤波器所有四种类型的最低阶数。指标如下：40kHz的采样率，4kHz的通带边界频率，8kHz的阻带边界频率，0.5dB的通带波纹，40dB的最小阻带衰减。并在同一张图中画出每种滤波器的频率响应。（1）MATLAB程序：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(wp,ws,rp,rs,'s

26、9;);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 on;plot(w,h,'r-');disp(&

27、#39;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 on;plot(w,h,'k-');disp('ellipord');n,wc=ellipord(wp,ws,rp,rs,'s');b,a=ellip(n,rp,rs,wc,'low','s');h,w=freqs(b,a,w)

28、;h=20*log10(abs(h);hold on;plot(w,h,'g-');legend('butter','cheby1','cheby2','ellip')grid on;xlabel('w');ylabel('h');运行结果：buttordn =9wc =0.7533cheb1ordn =5wpo = 0.6283cheb2ordn = 5wso = 1.2069ellipordn =4wc =0.6283（2）用matlab确定一个数字IIR高通滤波器所有四种类型的

29、最低阶数。指标如下：3500Hz的采样率，1050Hz的通带边界频率，600Hz的阻带边界频率，1dB的通带波纹，50dB的最小阻带衰减。并在同一张图中画出每种滤波器的频率响应。MATLAB程序：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=20*

30、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 on;plot(w,h,'r-');disp('cheb2ord');n,wso=cheb2ord(wp,ws,rp,rs,'s');b,a=cheby2(n,rs,wso,

31、'high','s');h,w=freqs(b,a,w);h=20*log10(abs(h);hold on;plot(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 on;plot(w,h,'g-');legend('butter',

32、9;cheby1','cheby2','ellip');axis(0,5,-120,0);grid on;xlabel('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的阻带边

33、界频率，0.4dB的通带波纹，50dB的最小阻带衰减。并在同一张图中画出每种滤波器的频率响应。MATLAB程序：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('ch

34、eb1ord');n,wpo=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 on;plot(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 on;plot(

35、w,h,'k-');disp('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 on;plot(w,h,'g-');legend('butter','cheby1','cheby2','ellip');axis(0,5,-120,20);grid on;xlabel('w

36、9;);ylabel('h');运行结果：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的最小阻带衰减。并在同一张图中画出每种滤波器的频率响应。MATLAB程序：fc=12;fp=2.1

37、,4.5;fs=2.7,3.9;rp=0.5;rs=40;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=che

38、by1(n,rp,wpo,'stop','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&#

39、39;)n,wc=ellipord(wp,ws,rp,rs,'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')运行结果：bu

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

41、MATLAB程序：wc=0.3*pi; M=11,41,81,121; for i=1:1:4tao=(M(i)-1)/2; h_d1n=zeros(1,M(i); for n=1:1:M(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');grid on;title( 'a = ', num2str(i) );end运行结果：（2）所加窗函数；MATLAB程序：M=11

42、,41,81,121; for i=1:1:4n=0:1:M(i)-1;wn_boxcar=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_bartlet

43、t,'r.');title('M=', num2str(i) ,' 三角窗函数'); 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,'

44、;r.');title('M=', num2str(i) ,' 布莱克曼窗函数'); end运行结果： （2）加窗后的滤波器冲激响应（2） 滤波器的幅频特性；MATLAB程序：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

45、(i);h_blackman=fir1(M,wc,blackman(N(i);n=0:1:N(i)-1; 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,

46、'.');title('N=',num2str(N(i), '汉宁窗函数冲击响应'); 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);dBH

47、x_boxcar=20*log10(abs(Hx_boxcar)/max(abs(Hx_boxcar); 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_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)');xlab