西电数字信号处理上机实验报告
数字信号处理上机实验报告
14020710021 张吉凯
第一次上机
实验一：
设给定模拟信号()1000t a x t e -=，t 的单位是ms 。
(1) 利用MATLAB 绘制出其时域波形和频谱图（傅里叶变换），估计其等效带宽（忽略谱分量降低到峰值的3%以下的频谱）。 (2) 用两个不同的采样频率对给定的()a x t 进行采样。 ○1()()15000s a f x t x n =以样本秒采样得到。
()()11j x n X e ω画出及其频谱。 ○2()()11000s a f x t x n =以样本秒采样得到。
()()
11j x n X e ω画出及其频谱。
比较两种采样率下的信号频谱，并解释。
（1）MATLAB 程序：
N=10; Fs=5; T s=1/Fs;
n=[-N:T s: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)时域波形');
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^(j\omega))');
ind = find(X >=0.03*max(X))*0.01; eband = (max(ind) -min(ind)); fprintf('等效带宽为%fKHZ\n',eband); 运行结果：
等效带宽为12.110000KHZ
（2）MATLAB程序：
N=10;
omega=-3*pi:0.01:3*pi;
%Fs=5000
Fs=5;
T s=1/Fs;
n=-N:T s: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));
title('频谱图(f_s=5000)');
xlabel('\omega/\pi');ylabel('X_1(f)');
grid on;
%Fs=1000
Fs=1;
T s=1/Fs;
n=-N:T s: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 on;
subplot(2,2,4);
plot(omega/pi,abs(X));
title('频谱图(f_s=1000)');
xlabel('\omega/\pi');
ylabel('X_2(f)');
grid on;
运行结果：
实验二：
给定一指数型衰减信号()()0cos 2at x t e f t π-=，采样率1
s f T
=，T 为采样周期。为方便起见，重写成复指数形式()02j f t at x t e e π-=。
采样后的信号为()02j f nT anT x nT e e π-=，加窗后长度为L 的形式为：
()(),0,1,,1L x nT x nT n L ==-L
这3个信号()x t ，()x nT ，()L x nT 的幅度谱平方分别为： 模拟信号：
()()()
2
2
2
01
2X f a f f π=
+-
采样信号：
()()()2
201
?12cos 2aT aT
X
f e f f T e π--=--+
加窗（取有限个采样点）信号：
()()()()()22
02012cos 2?12cos 2aTL aTL L
aT
aT
e f f TL e X f e
f f T e
ππ------+=--+
且满足如下关系： ()()()()???lim ,lim s L
L f X f X f TX f X f →∞
→∞
== 实验内容
100.2sec ,0.5Hz,1Hz 2Hz =10s s a f f f L -====取采样频率分别取和，。
（1） 在同一张图上画出：模型号幅度谱平方()2
X f ；
()()2
?1Hz 2Hz 0Hz 3Hz s s f f TX
f f ==≤≤和时，采样信号幅度谱平方
（2） 在同一张图上画出：模型号幅度谱平方()2
X f ；
()()2
?2Hz 0Hz 3Hz s f TX
f f =≤≤时，采样信号幅度谱平方；改变L 值，
结果又如何？
（1）MATLAB 程序：
f=0:0.01:3; alpha=0.2; f0=0.5; L=10; T1=1; T2=0.5;
Xa=1./(alpha^2+(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(-alpha*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');
运行结果：
（2）MATLAB程序：
f=0:0.01:3;
alpha=0.2;
f0=0.5;
L1=5;
L2=10;
L3=20;
T1 = 0.5
Xa=1./(alpha^2+(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*L2))./( 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,'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)');
运行结果：
实验三：
设(){}11,2,2x n =，(){}21,2,3,4x n =，编写MATLAB 程序，计算： （１） ５点圆周卷积()1y n ； （２） ６点圆周卷积()2y n ； （３） 线性卷积()3y n ；
（４） 画出的()1y n ，()2y n 和()3y n 时间轴对齐。
MATLAB 程序：
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 on
title('五点圆周卷积y1(n)');
xlabel('n'),ylabel('y1(n)');axis([0 6 0 15])
subplot(312)
stem(y2);
grid on
title('六点圆周卷积y2(n)');
xlabel('n'),ylabel('y2(n)');axis([0 6 0 15])
subplot(313)
stem(y3);
grid on
title('线性卷积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]);
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]);
运行结果：
实验四：
给定因果系统：()()()0.91y n y n x n =-+ （１） 求系统函数()H z 并画出零极点示意图。 （２） 画出系统的幅频特性()
j H e ω和相频特性()?ω。 （３） 求脉冲响应()h n 并画序列图。
提示：在ＭＡＴＬＡＢ中，zplane(b,a) 函数可画零极点图；Freqz(b,a,N)可给出[]0,π范围内均匀间隔的N 点频率响应的复振幅；Impz(b,a,N)可求()H z 的逆变换（即脉冲响应）。
MATLAB 程序：
a = [1,0]
b = [1,-0.9] figure(1) zplane(b,a);
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); impz(b,a);
运行结果：
第二次上机
1. 给定模拟信号()()()2sin 45cos 8x t t t ππ=+，对其进行采样，用DFT （FFT ）进行信号频谱分析。
（１） 确定最小采样频率和最小采样点数。
（２） 若以()0.010:1t n n N ==-秒进行采样，至少需要取多少采样点？ （３） 用DFT 的点数50,100N =画出信号的N 点DFT 的幅度谱，讨论幅度谱结果。
（４） N 分别为64N =和60N =，能否分辨出信号的所有频率分量。 （５） 在（３）和（４）的条件下做补0 FFT ，分析结果。
（６） 在不满足最小采样点数的情况下做补0DFT ，观察是否可以分辨出两
个频率分量。
（1）最小采样频率：8 ；最小采样点数：4 （2） 最小采样点数：50 （3）（4）MATLAB 程序：
N1=50;N2=100;N3=64;N4=60;
n1=0:N1-1; n2=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,'.');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'); subplot(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后的幅度值');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,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');ylabel('X9'); 运行结果：
2. 设雷达发射线性调频信号()()2
exp 2h t j t πμ=，13510μ=?，采样率
9210s f =?，采
样点数20000N =。回波信号()()()1
2
s t h t h t ττ=-+-，6110τ-=，62 1.110τ-=?。 （１） 画出()h t 的频谱。
（２） 利用DFT 的时延性质产生()s t ，比较直接在时域产生和在频域产生（再
变换到时域）的结果是否相同。
（３） 匹配滤波的结果是()()()y t s t h t *
=*-，（“*”表示线性卷积）。分别用直
接线性卷积和DFT 的卷积定理求解()y t 。比较二者结果，并记录两种方法的运行时间（用tic ，toc 指令）。 （４） 画出()y t 的频谱。
（1）MATLAB 程序：
mu=5e13;
fs=2e9;
T s=1/fs;
mu=5e13;
fs=2e9;
T s=1/fs;
N=20000;
tao1=1e-6;
tao2=1.1e-6;
delay1=ceil(tao1/T s); delay2=ceil(tao2/T s); 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;
运行结果：
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;
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;
运行结果：
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-N*T s; 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.
figure(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');
[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('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);
h=20*log10(abs(h));
hold on;
plot(w,h,'g-');
legend('butter','cheby1','cheby2','ellip')
grid on;
xlabel('w');
ylabel('h');