数字信号处理上机实验.docx
《数字信号处理上机实验.docx》由会员分享，可在线阅读，更多相关《数字信号处理上机实验.docx（24页珍藏版）》请在冰豆网上搜索。
数字信号处理上机实验
《数字信号处理》上机实验报告
院系:
班级:
学号：
姓名：
实验一：
信号、系统及系统响应····················1
实验二：
用FFT作谱分析··························7
实验三：
用双线性变换法设计IIR数字滤波器········15
四：
总结·······································18
实验一：
信号、系统及系统响应
1.实验目的
（1）熟悉联系信号经理想采样后的频谱变化关系，加深对时域采样定理的理解
（2）熟悉时域离散系统的时域特性。
（3）利用卷积方法观察分析系统的时域特性。
（4）掌握序列傅里叶变换的计算机实现方法，利用序列的傅里叶变换对连续信号、离散信号及系统响应进行频域分析。
2.实验原理与方法
（1）时域采样。
（2）LTI系统的输入输出关系。
3.实验内容及步骤
（1）认真复习采样理论、离散信号与系统、线性卷积、序列的傅里叶变换及性质等有关内容，阅读本实验原理与方法。
（2）编制实验用主程序及相应子程序。
①信号产生子程序，用于产生实验中要用到的下列信号序列：
a.xa（t）=A*e^-at*sin（Ω0t）u（t）
b.单位脉冲序列：
xb（n）=δ（n）
c.矩形序列：
xc（n）=RN（n）,N=10
②系统单位脉冲响应序列产生子程序。
本实验要用到两种FIR系统。
a.ha（n）=R10（n）;
b.hb（n）=δ（n）+2.5δ（n-1）+2.5δ（n-2）+δ（n-3）
③有限长序列线性卷积子程序
用于完成两个给定长度的序列的卷积。
可以直接调用MATLAB语言中的卷积函数conv。
conv用于两个有限长度序列的卷积，它假定两个序列都从n=0开始。
调用格式如下：
y=conv（x,h）
（3）调通并运行实验程序，完成下述实验内容
①分析采样序列的特性。
采样信号xa（n）的参数为A=444.128,a=50
=50
。
a.取采样频率fs=1kHz,，即T=1ms。
b.改变采样频率，fs=300Hz，观察|X（e^jω）|的变化，并做记录（打印曲线）；进一步降低采样频率，fs=200Hz，观察频谱混叠是否明显存在，说明原因，并记录（打印）这时的|X（e^jω）|曲线。
程序代码如下：
function[y]=x（t,A,a,omiga0）
y=A*exp（-a*t）.*sin（omiga0*t）.*u（t）;
end
%分析采样序列特性
A=444.128;a=50*2^0.5*pi;omiga0=50*2^0.5*pi;
fs1=1000;fs2=300;fs3=200;
w=linspace（-2*pi,2*pi,10000）;
n=0:
49;
x1=x（n/fs1,A,a,omiga0）;
x2=x（n/fs2,A,a,omiga0）;
x3=x（n/fs3,A,a,omiga0）;
X1=x1*exp（-j*n'*w）;
X2=x2*exp（-j*n'*w）;
X3=x3*exp（-j*n'*w）;
subplot（3,2,1）;stem（n,x1,'.'）;ylabel（'y'）;xlabel（'n'）;
title（'时间函数'）;
subplot（3,2,2）;plot（w/pi,abs（X1）,'r'）;ylabel（'|X（jf）|'）;xlabel（'\omega/\pi'）;
title（'频谱图'）;text（1.5,1200,'f=1000Hz'）;
subplot（3,2,3）;stem（n,x2,'.'）;ylabel（'y'）;xlabel（'n'）;
title（'时间函数'）;
subplot（3,2,4）;plot（w/pi,abs（X2）,'r'）;ylabel（'|X（jf）|'）;xlabel（'\omega/\pi'）;
title（'频谱图'）;text（1.5,500,'f=300Hz'）;
subplot（3,2,5）;stem（n,x3,'.'）;ylabel（'y'）;xlabel（'n'）;
title（'时间函数'）;
subplot（3,2,6）;plot（w/pi,abs（X3）,'r'）;ylabel（'|X（jf）|'）;xlabel（'\omega/\pi'）;
title（'频谱图'）;text（1.5,250,'f=200Hz'）;
图1对xa不同取样频率得到的取样图和频谱图
由图可知，采样频率不同，所得到的采样后信号和其傅里叶变换都不同。
fs=1000Hz时，频谱的混叠效应很小，fs=300Hz时，混叠效应加大，fs=200Hz时，混叠效应进一步加大。
这是因为采样频率越来越接近临界采样频率，所以会造成频谱的混叠。
②时域离散信号、系统和系统响应分析。
a.观察信号xb（n）和系统hb（n）的时域和频域特性；利用线性卷积求信号xb（n）通过系统hb（n）的响应y（n），比较所求响应y（n）和hb（n）的时域及频域特性，注意它们之间有无差别，绘图说明，并用所学理论解释所得结果。
b.观察系统ha（n）对信号xc（n）的响应特性。
程序代码如下：
a.
clearall;
xb=double（impact（0,0,1））;
xc=single（R（1,1,10））;
ha=single（R（1,1,10））;
hb=impact（0,0,3）+2.5*impact（1,0,3）+2.5*impact（2,0,3）+impact（3,0,3）;
[Xb,w]=DFT（xb,2）;
[Hb,w]=DFT（hb,4）;
y=conv（xb,hb）;
[Y,w]=DFT（y,5）;
subplot（3,2,1）;stem（0:
1,xb,'.'）;ylabel（'xb'）;xlabel（'n'）;
title（'时间函数'）;
subplot（3,2,2）;plot（w/pi,Xb,'r'）;ylabel（'|Xb（jf）|'）;xlabel（'\omega/\pi'）;
title（'频谱图'）;text（1.5,700,'f=1000Hz'）;
subplot（3,2,3）;stem（hb,'.'）;ylabel（'hb'）;xlabel（'n'）;
title（'时间函数'）;
subplot（3,2,4）;plot（w/pi,Hb,'r'）;ylabel（'|Hb（jf）|'）;xlabel（'\omega/\pi'）;
title（'频谱图'）;text（1.5,700,'f=300Hz'）;
subplot（3,2,5）;stem（y,'.'）;ylabel（'y'）;xlabel（'n'）;
title（'时间函数'）;
subplot（3,2,6）;plot（w/pi,Y,'r'）;ylabel（'|Y（jf）|'）;xlabel（'\omega/\pi'）;
title（'频谱图'）;text（1.5,700,'f=300Hz'）;
图2xb、hb的时域频域特性
δ（n）的傅里叶变换恒为1。
由图可知，xb（n）*hb（n）=hb（n），即y（n）=hb（n），所以y（n）和hb（n）的频谱也是完全一样的。
b.
clearall;
xc=single（R（1,1,10））;%xc=R10（n）
ha=single（R（1,1,10））;
yn=conv（ha,xc）;
[Y,w]=DFT（yn,19）;
subplot（2,2,1）;stem（yn,'.'）;ylabel（'yn'）;xlabel（'n'）;
title（'时间函数'）;
subplot（2,2,2）;plot（w/pi,abs（Y）,'r'）;ylabel（'|Y（jf）|'）;xlabel（'\omega/\pi'）;
title（'频谱图'）;
xc=single（R（1,1,5））;%xc=R5（n）
ha=single（R（1,1,10））;
yn=conv（ha,xc）;
[Y,w]=DFT（yn,14）;
subplot（2,2,3）;stem（yn,'.'）;ylabel（'yn'）;xlabel（'n'）;
title（'时间函数'）;
subplot（2,2,4）;plot（w/pi,abs（Y）,'r'）;ylabel（'|Y（jf）|'）;xlabel（'\omega/\pi'）;
title（'频谱图'）;
图3xc长度为10和5时与ha的卷积时序图和频域图
y（n）的长度与理论计算的相同。
两序列卷积后的长度为n1+n2-1，用此方法可以快速验证两个序列的卷积结果是否正确。
两次卷积的长度分别为10+10-1=19和10+5-1=14由图可以判断，卷积的结果正确。
③卷积定理的验证。
A=1,a=0.4,
=2.0734，T=1。
程序代码如下：
n=0:
49;
xa=x（n,1,0.4,2.0374）;
hb=impact（0,0,3）+2.5*impact（1,0,3）+2.5*impact（2,0,3）+impact（3,0,3）;
[Xa,w]=DFT（xa,50）;
yn=conv（xa,hb）;
Y=DFT（yn,53）;
figure
（1）;
subplot（2,2,1）;stem（yn,'.'）;ylabel（'yn'）;xlabel（'n'）;
title（'时间函数'）;
subplot（2,2,3）;plot（w/pi,abs（Y）,'r'）;ylabel（'|Y（jf）|'）;xlabel（'\omega/\pi'）;
title（'频谱图'）;
subplot（2,2,4）;plot（w/pi,angle（Y）,'r'）;ylabel（'相位'）;xlabel（'\omega/\pi'）;
title（'Y（jf）的相位'）;
Hb=DFT（hb,4）;
Yn=Hb.*Xa;
figure
（2）;
subplot（1,2,1）;plot（w/pi,abs（Yn）,'r'）;ylabel（'|Yn（jf）|'）;xlabel（'\omega/\pi'）;
title（'频谱图'）;
subplot（1,2,2）;plot（w/pi,angle（Yn）,'r'）;ylabel（'相位'）;xlabel（'\omega/\pi'）;
title（'Yn（jf）的相位'）;
先对xa,hb求卷积再转换为频谱图如下：
图4xa时域图、xa与hb卷积频域图和对应相位图
Hb.*Xa的频谱图如下
图5Hb.*Xa的频谱图
对比图4图5即可验证卷积定理的正确性。
4、思考题
（1）在分析理想采样序列特性的实验中，采样频率不同时，相应理想采样序列的傅里叶变换频谱的数字频率度量是否都相同？
它们所对应的模拟频率是否相同？
为什么？
答：
数字频率度量不相同，但他们所对应的模拟频率相同。
由公式w=Ω*Ts可知，采样间隔Ts的变化会引起数字频率w的变化，但是不会引起模拟频率Ω的变化。
（2）在卷积定理验证的实验中，如果选用不同的频域采样点数M值，例如，选M=10和M=20，分别做序列的傅里叶变换，求得的结果有无差异？
答：
有差异，采样点数不一样，得到的傅里叶变换的点数也会不一样。
实验二：
用FFT作谱分析
1、实验目的
（1）进一步加深DFT算法原理和基本性质的理解（因为FFT只是DFT的一种快速算法，所以FFT的运算结果必然满足DFT的基本性质）。
（2）熟悉FFT算法原理和FFT子程序的应用。
（3）学习用FFT对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差及其原因，以便在实际中正确应用FFT。
2、实验步骤
（1）复习DFT的定义、性质和用DFT作谱分析的有关内容。
（2）复习FFT算法原理与编程思想，并对照DIT-FFT运算流图和程序框图，读懂本实验提供的FFT子程序。
（3）编制信号产生子程序，产生以下典型信号供谱分析用：
（4）编写主程序。
（5）按实验内容要求，上机实验，并写出实验报告。
3、实验内容
（1）对2中所给出的信号逐个进行谱分析。
程序代码如下：
x1=[11110000];
x2=[12344321];
x3=[43211234];
fs=64;
N1=8;N2=16;N3=32;N4=64;
%x4=cos（0.25*pi*n）;
%x5=sin（（pi/8）*n）;
%x7=cos（0.25*pi*n）+sin（（pi/8）*n）;
%x8=cos（0.25*pi*n）+j*sin（（pi/8）*n）;
%x1作图
i=1;
f1=fft（x1,N1）;
figure（i）;i=i+1;
subplot（1,3,1）;stem（0:
N1-1,x1,'.'）;title（'x1n的波形'）;
ylabel（'x1n'）;xlabel（'n'）;
subplot（1,3,2）;stem（0:
N1-1,abs（f1）,'.'）;title（'x1n的8点FFT'）;
ylabel（'|X（k）|'）;xlabel（'k'）;
f1=fft（x1,N2）;
subplot（1,3,3）;stem（0:
N2-1,abs（f1）,'.'）;title（'x1n的16点FFT'）;
ylabel（'|X（k）|'）;xlabel（'k'）;
%x2作图
f2=fft（x2,N1）;
figure（i）;i=i+1;
subplot（1,3,1）;stem（0:
N1-1,x2,'.'）;title（'x2n的波形'）;
ylabel（'x2n'）;xlabel（'n'）;
subplot（1,3,2）;stem（0:
N1-1,abs（f2）,'.'）;title（'x2n的8点FFT'）;
ylabel（'|X（k）|'）;xlabel（'k'）;
f2=fft（x2,N2）;
subplot（1,3,3）;stem（0:
N2-1,abs（f2）,'.'）;title（'x2n的16点FFT'）;
ylabel（'|X（k）|'）;xlabel（'k'）;
%x3作图
f3=fft（x3,N1）;
figure（i）;i=i+1;
subplot（1,3,1）;stem（0:
N1-1,x3,'.'）;title（'x3n的波形'）;
ylabel（'x3n'）;xlabel（'n'）;
subplot（1,3,2）;stem（0:
N1-1,abs（f3）,'.'）;title（'x3n的8点FFT'）;
ylabel（'|X（k）|'）;xlabel（'k'）;
f3=fft（x3,N2）;
subplot（1,3,3）;stem（0:
N2-1,abs（f3）,'.'）;title（'x3n的16点FFT'）;
ylabel（'|X（k）|'）;xlabel（'k'）;
%x4作图
n=0:
N1-1;
x4=sin（（pi/8）*n）;
f4=fft（x4,N1）;
figure（i）;i=i+1;
subplot（2,2,1）;stem（0:
N1-1,x4,'.'）;title（'x4n的8点波形'）;
ylabel（'x4n'）;xlabel（'n'）;
subplot（2,2,2）;stem（0:
N1-1,abs（f4）,'.'）;title（'x4n的8点FFT'）;
ylabel（'|X（k）|'）;xlabel（'k'）;
n=0:
N2-1;
x4=cos（0.25*pi*n）;
f4=fft（x4,N2）;
subplot（2,2,3）;stem（0:
N2-1,x4,'.'）;title（'x4n的16点波形'）;
ylabel（'x4n'）;xlabel（'n'）;
subplot（2,2,4）;stem（0:
N2-1,abs（f4）,'.'）;title（'x4n的16点FFT'）;
ylabel（'|X（k）|'）;xlabel（'k'）;
%x5作图
figure（i）;i=i+1;
n=0:
N1-1;
x5=cos（0.25*pi*n）;
f5=fft（x5,N1）;
subplot（2,2,1）;stem（0:
N1-1,x5,'.'）;title（'x5n的8点波形'）;
ylabel（'x5n'）;xlabel（'n'）;
subplot（2,2,2）;stem（0:
N1-1,abs（f5）,'.'）;title（'x5n的8点FFT'）;
ylabel（'|X（k）|'）;xlabel（'k'）;
n=0:
N2-1;
x5=cos（0.25*pi*n）;
f5=fft（x5,N2）;
subplot（2,2,3）;stem（0:
N2-1,x5,'.'）;title（'x5n的16点波形'）;
ylabel（'x5n'）;xlabel（'n'）;
subplot（2,2,4）;stem（0:
N2-1,abs（f5）,'.'）;title（'x5n的16点FFT'）;
ylabel（'|X（k）|'）;xlabel（'k'）;
%x6作图
figure（i）;i=i+1;
n=0:
N2-1;
x=cos（8*pi*n/fs）+cos（16*pi*n/fs）+cos（20*pi*n/fs）;
f6=fft（x,N2）;
subplot（3,2,1）;stem（0:
N2-1,x,'.'）;title（'x6n的16点波形'）;
ylabel（'x6n'）;xlabel（'n'）;
subplot（3,2,2）;stem（0:
N2-1,abs（f6）,'.'）;title（'x6n的16点FFT'）;
ylabel（'|X（k）|'）;xlabel（'k'）;
n=0:
N3-1;
x=cos（8*pi*n/fs）+cos（16*pi*n/fs）+cos（20*pi*n/fs）;
f6=fft（x,N3）;
subplot（3,2,3）;stem（0:
N3-1,x,'.'）;title（'x6n的32点波形'）;
ylabel（'x6n'）;xlabel（'n'）;
subplot（3,2,4）;stem（0:
N3-1,abs（f6）,'.'）;title（'x6n的32点FFT'）;
ylabel（'|X（k）|'）;xlabel（'k'）;
n=0:
N4-1;
x=cos（8*pi*n/fs）+cos（16*pi*n/fs）+cos（20*pi*n/fs）;
f6=fft（x,N4）;
subplot（3,2,5）;stem（0:
N4-1,x,'.'）;title（'x6n的64点波形'）;
ylabel（'x6n'）;xlabel（'n'）;
subplot（3,2,6）;stem（0:
N4-1,abs（f6）,'.'）;title（'x6n的64点FFT'）;
ylabel（'|X（k）|'）;xlabel（'k'）;
运行结果如下:
图6x1的8、16点谱分析
图7x2的8、16点谱分析
图8x3的8、16点谱分析
图9x4的8、16点谱分析
图10x5的8、16点谱分析
图11x6的16、32、64点谱分析
（2）令x7（n）=x4（n）+x5（n），用FFT计算8点和16点离散傅里叶变换。
程序代码如下：
%x7作图
figure（i）;i=i+1;
n=0:
N1-1;
x7=sin（（pi/8）*n）+cos（0.25*pi*n）;
f7=fft（x7,N1）;
t1=max（x7）;
t2=max（f7）;
subplot（3,2,1）;stem（0:
N1-1,x7,'.'）;title（'x7n的8点波形'）;
ylabel（'x7n'）;xlabel（'n'）;
subplot（3,2,2）;stem（0:
N1-1,abs（f7）,'.'）;title（'x7n的8点FFT'）;
ylabel（'|X（k）|'）;xlabel（'k'）;
n=0:
N2-1;
x7=sin（（pi/8）*n）+cos（0.25*pi*n）;
f7=fft（x7,N2）;
t1=max（x7）;
t2=max（f7）;
subplot（3,2,3）;stem（0:
N2-1,x7,'.'）;title（'x7n的16点波形'）;
ylabel（'x7n'）;xlabel（'n'）;
subplot（3,2,4）;stem（0:
N2-1,abs（f7）,'.'）;title（'x7n的16点FFT'）;
ylabel（'|X（k）|'）;xlabel（'k'）;
k=conj（f7）;
x4=（k+f7）/2;
subplot（3,2,5）;stem（0:
N2-1,abs（x4）,'.'）;title（'恢复的X4（K）'）;
ylabel（'|Re（X7（k））|'）;xlabel（'k'）;
x5=（k-f7）/2;
subplot（3,2,6）;stem（0:
N2-1,abs（x5）,'.'）;title（'恢复的X5（K）'）;
ylabel（'|jIm（X7（k））|'）;xlabel（'k'）;
运行结果
图12左1：
x7的8点时域图右1：
x7的8点频谱图
左2：
x7的16点时域图右2：
x7的16点频谱图
左3：
对ifft（Re（X7））的到的X4（k）频谱
右3：
对ifft（Im（X7））的到的X5（k）频谱
将上图左3和右3与图4图5中的16点频谱图对比可知两者对应完全相等，验证了DFT的如下性质：
xnepRe（Xn）xnopjIm（Xn）。
（3）令x8（n）=x4（n）+j*x5（n），重复
（2）。
程序代码如下：
%x8作图
figure（i）;i=i+1;
n=0:
N1-1;
x8=cos（0.25*pi*n）+j*sin（（pi/8）*n）;
f8=fft（x8,N1）;
subplot（3,2,1）;stem（0:
N1-1,x8,'.'）;title（'x8n的8点波形'）;
ylabel（'x8n'）;xlabel（'n'）;
subplot（3,2,2）;stem（0:
N1-1,abs（f8）,'.'）;title（'x8n的8点FFT'）;
ylabel（'|X（k）|'）;xlabel（'k'）;
n=0:
N2-1;
x8=cos（0.25*pi*n）+j*sin（（pi/8）*n）;
f8=fft（x8,N2）;
subplot（3,2,3）;stem（0:
N2-1,x8,'.'）;title（'x8n的16点波形'）;
ylabel（'x8n'）;xlabel（'n'）;
subplot（3,2,4）;stem（0:
N2-1,abs（f8）,'.'）;title（'x8n的16点FFT'）;
ylabel（'|X（k）|'）;xlabel（'k'）;
k
（1）=conj（f8
（1））;
form=2:
N2
k（
升级会员 