数字信号处理仿真实验.docx
《数字信号处理仿真实验.docx》由会员分享，可在线阅读，更多相关《数字信号处理仿真实验.docx（28页珍藏版）》请在冰豆网上搜索。
数字信号处理仿真实验
数字信号处理上机实验仿真报告
学院电子工程学院
班级
学号
姓名
实验一：
信号、系统及系统响应
一、实验目的
（1）熟悉连续信号经理想采样前后的频谱变化关系，加深对时域采样定理的理解。
（2）熟悉时域离散系统的时域特性。
（3）利用卷积方法观察分析系统的时域特性。
（4）掌握序列傅里叶变换的计算机实现方法，利用序列的傅里叶变换对连续信号、离散信号及系统响应进行频域分析。
二、实验原理
采样是连续信号数字处理的第一个关键环节。
对一个连续信号
进行理想采样的过程可用（1.1）式表示。
（1.1）
其中
为
的理想采样，
为周期冲激脉冲，即
（1.2）
的傅里叶变换
为
（1.3）
将（1.2）式代入（1.1）式并进行傅里叶变换，
（1.4）
式中的
就是采样后得到的序列
，即
的傅里叶变换为
（1.5）
比较（1.5）和（1.4）可知
（1.6）
为了在数字计算机上观察分析各种序列的频域特性，通常对
在
上进行M点采样来观察分析。
对长度为N的有限长序列
，有
（1.7）
其中
一个时域离散线性时不变系统的输入/输出关系为
（1.8）
上述卷积运算也可以转到频域实现
（1.9）
三、实验内容及步骤
（1）认真复习采样理论、离散信号与系统、线性卷积、序列的傅里叶变换及性质等有关内容，阅读本实验原理与方法。
（2）编制实验用主程序及相应子程序。
①信号产生子程序，用于产生实验中要用到的下列信号序列：
xa（t）=Ae-atsin（Ω0t）u（t）
进行采样，可得到采样序列
xa（n）=xa（nT）=Ae-anTsin（Ω0nT）u（n）,0≤n<50
其中A为幅度因子，a为衰减因子，Ω0是模拟角频率，T为采样间隔。
这些参数都要在实验过程中由键盘输入，产生不同的xa（t）和xa（n）。
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）调通并运行实验程序，完成下述实验内容：
①分析采样序列的特性。
a.取采样频率fs=1kHz,即T=1ms。
b.改变采样频率，fs=300Hz，观察|X（ejω）|的变化，并做记录（打印曲线）；进一步降低采样频率，fs=200Hz，观察频谱混叠是否明显存在，说明原因，并记录（打印）这时的|X（ejω）|曲线。
②时域离散信号、系统和系统响应分析。
a.观察信号xb（n）和系统hb（n）的时域和频域特性；利用线性卷积求信号xb（n）通过系统hb（n）的响应y（n），比较所求响应y（n）和hb（n）的时域及频域特性，注意它们之间有无差别，绘图说明，并用所学理论解释所得结果。
b.观察系统ha（n）对信号xc（n）的响应特性。
③卷积定理的验证。
（4）主程序框图
①分析采样序列的特性
②时域离散信号、系统和系统响应分析
③卷积定理的验证
四．实验程序及对应波形
1.子程序
function[XN,n,k]=DFT（xn,N）
n=0:
N-1;
k=-200:
200;
XN=xn*exp（-j*2*pi/N）.^（n'*k）;%计算DFT[x（n）]
%产生矩形序列
functionx=juxing（n2）;
x=[1,ones（1,n2）];
function[x,n]=maichong（n0,n1,n2）
n=（n1:
n2）;
x=（n==n0）;
%产生信号Xa（n）
functionx=xn（A,a,w,fs）
n=0:
50-1;
x=A*exp（（-a）*n/fs）.*sin（w*n/fs）.*juxing（49）;
functionx=u（t）;
x=（t>=0）;
%产生脉冲信号
function[x,n]=maichong（n0,n1,n2）
n=（n1:
n2）;
x=（n==n0）;
2.主程序
①分析采样序列的特性
程序代码如下（见后附程序Ex1_1）:
A=100;a=200;w0=200;k=-200:
200;
T=0.001;t=0:
T:
0.06;N=50;k1=0:
1:
N;W1max=2*pi*500;W1=W1max*k1/N;w1=W1/pi;
xat=A*exp（-a*t）.*sin（w0*t）.*u（t）;
Xa=xat*exp（-j*t'*W1）;
figure
（1）;
subplot（2,1,1）;
plot（t,xat）;
xlabel（'t'）;
ylabel（'xa（t）'）;
title（'采样信号a（t）'）;
axis（[0,0.06,-5,35]）;
subplot（2,1,2）;
plot（w1,abs（Xa））;
xlabel（'w'）;
ylabel（'X（jw）'）;title（'xa（t）的频谱'）;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%fs=1000Hz
A=100;a=200;w0=200;k=-200:
200;
fs=1000;
w=k/50;
xan=xn（A,a,w0,fs）;%产生信号xa（n）
X=DFT（xan,50）;
figure
（2）
subplot（2,1,1）
n=0:
49;
stem（n,xan,'.'）;axis（[0,50,-20,50]）;
xlabel（'n'）;
ylabel（'xa（n）'）;
title（'采样信号fs=1000Hz'）;
subplot（2,1,2）;
plot（w,abs（X））;
xlabel（'w/pi'）;
ylabel（'X（e^jw）'）;title（'xa（n）的频谱（fs=1000Hz）'）;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%fs=300Hz
fs=300;
xan=xn（A,a,w0,fs）;%产生信号xa（n）
X=DFT（xan,50）;
figure（3）
subplot（2,1,1）
n=0:
49;
stem（n,xan,'.'）;axis（[0,50,-20,50]）;
xlabel（'n'）;
ylabel（'xa（n）'）;
title（'采样信号fs=300Hz'）;
subplot（2,1,2）;
plot（w,abs（X））;
xlabel（'w/pi'）;
ylabel（'X（e^jw）'）;title（'xa（n）的频谱（fs=300Hz）'）;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%fs=200Hz
fs=200;
xan=xn（A,a,w0,fs）;%产生信号xa（n）
X=DFT（xan,50）;
figure（4）
subplot（2,1,1）
n=0:
49;
stem（n,xan,'.'）;axis（[0,50,-20,50]）;
xlabel（'n'）;
ylabel（'xa（n）'）;
title（'采样信号fs=200Hz'）;
subplot（2,1,2）;
plot（w,abs（X））;
xlabel（'w/pi'）;
ylabel（'X（e^jw）'）;title（'xa（n）的频谱（fs=200Hz）'）;
仿真结果如下:
由图可见，在折叠频率w=π,即f=fs/2=500Hz处混叠很小。
当fs=300Hz时，存在较明显的混叠失真；当fs=200时，发生严重的混叠失真。
②时域离散信号、系统和系统响应分析
a.观察信号xb（n）和系统hb（n）的时域和频域特性；
实验代码如下：
（见后附程序Ex1_2）
%信号xb
N=64
xb=zeros（1,N）
xb
（1）=1;xn1=[0:
N-1]
figure
（1）;subplot（2,1,1）
stem（xn1,xb,'.'）
title（'信号xb时域特性'）
w=linspace（-5*pi,5*pi,1000）
Xb=xb*exp（-j*xn1'*w）
figure
（1）
subplot（2,1,2）
plot（w/pi,abs（Xb）,'.'）
title（'信号xb频域特性'）
%ϵͳ信号hb
hb=[1,2.5,2.5,1];nhb=[0:
3]
figure
（2）;subplot（2,1,1）
stem（nhb,hb,'.'）
title（'系统函数ýhb时域特性'）
w=linspace（-2.8*pi,2.8*pi,1000）
Hb=hb*exp（-j*nhb'*w）
figure
（2）
subplot（2,1,2）
plot（w/pi,abs（Hb）,'.'）
title（'系统函数ýhb频域特性'）
%响应
yb=conv（xb,hb）;
nyb=[0:
length（yb）-1]
figure（3）;subplot（2,1,1）
stem（yb,'.'）
title（'响应yb时域特性'）
w=linspace（-2.8*pi,2.8*pi,1000）
Yb=yb*exp（-j*nyb'*w）
figure（3）
subplot（2,1,2）
plot（w/pi,abs（Yb）,'.'）
title（'响应Yb频域特性'）
仿真结果如下:
b.观察系统ha（n）对信号xc（n）的响应特性。
实验代码如下（见后附程序Ex1_3）:
%信号xc
N=10
xc=ones（1,N）
xn1=[0:
N-1]
figure
（1）;subplot（2,1,1）
stem（xn1,xc,'.'）
title（'信号xc时域响应'）
w=linspace（-5*pi,5*pi,1000）
Xc=xc*exp（-j*xn1'*w）
figure
（1）
subplot（2,1,2）
plot（w/pi,abs（Xc））
%ϵͳ信号ha
N=10
ha=ones（1,N）
hn1=[0:
N-1]
figure
（2）;subplot（2,1,1）
stem（hn1,ha,'.'）
title（'系统函数ha时域响应'）
w=linspace（-2.8*pi,2.8*pi,1000）
Ha=ha*exp（-j*hn1'*w）
figure
（2）
subplot（2,1,2）
plot（w/pi,abs（Ha）,'.'）
title（'系统函数ha频域响应'）
%响应
y2=conv（xc,ha）;
ny2=[0:
length（y2）-1]
figure（3）;subplot（2,1,1）
stem（y2,'.'）
title（'响应y2时域响应'）
w=linspace（-2.8*pi,2.8*pi,1000）
Y2=y2*exp（-j*ny2'*w）
figure（3）
subplot（2,1,2）
plot（w/pi,abs（Y2）,'.'）
title（'响应Yb频域响应'）
仿真结果如下:
③卷积定理的验证。
实验代码如下：
（见后附程序Ex1_4）:
%时域卷积定理的验证
%信号xc
N=10
xc=ones（1,N）
xn1=[0:
N-1]
w=linspace（-5*pi,5*pi,1000）
Xc=xc*exp（-j*xn1'*w）
%系统信号ha
N=10
ha=ones（1,N）
hn1=[0:
N-1]
w=linspace（-5*pi,5*pi,1000）
Ha=ha*exp（-j*hn1'*w）
%响应
y2=conv（xc,ha）;
ny2=[0:
length（y2）-1]
w=linspace（-5*pi,5*pi,1000）
Y2=y2*exp（-j*ny2'*w）
figure
（1）
subplot（2,1,1）
plot（w/pi,abs（Y2））
title（'响应Yb=DTFT（conv（xc,ha））频域特性'）
Y=Xc.*Ha
subplot（2,1,2）
plot（w/pi,abs（Y））
title（'响应Y=Xc.*Ha频域特性'）
仿真结果如下：
两种方法得到的谱是一样的，即验证了卷积定理。
五、思考题
（1）在分析理想采样序列特性的实验中，采样频率不同时，相应理想采样序列的傅里叶变换频谱的数字频率度量是否都相同?
它们所对应的模拟频率是否相同?
为什么?
答：
由
可知，若采样频率不同，则其周期T不同，相应的数字频率
也不相同；而因为是同一信号，故其模拟频率
保持不变。
（2）在卷积定理验证的实验中，如果选用不同的频域采样点数M值，例如，选M=10和M=20，分别做序列的傅里叶变换，求得
所得结果之间有无差异?
为什么?
答：
有差异。
因为所得
图形由其采样点数唯一确定，由频域采样定理可知，若M小于采样序列的长度N，则恢复原序列时会发生时域混叠现象。
实验二：
用FFT作谱分析
一、实验目的：
（1）进一步加深DFT算法原理和基本性质的理解（因为FFT只是DFT的一种快速算法，所以FFT的运算结果必然满足DFT的基本性质）。
（2）熟悉FFT算法原理和FFT子程序的应用。
（3）学习用FFT对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差及其原因，以便在实际中正确应用FFT。
二、实验步骤：
（1）复习DFT的定义、性质和用DFT作谱分析的有关内容。
（2）复习FFT算法原理与编程思想，并对照DIT—FFT运算流图和程序框图，读懂本实验提供的FFT子程序。
（3）编制信号产生子程序，产生以下典型信号供谱分析用：
=R4（n）
=cos（pi/4*n）
=sin（pi/8*n）
=cos（pi*8*t）+sin（pi*16*t）+cos（20*pi*t）
应当注意，如果给出的是连续信号xa（t）,则首先要根据其最高频率确定采样速率fs以及由频率选择采样点数N，然后对其进行软件采样（即计算x（n）=xa（nT）,0<=n<=N-1）,产生对应序列x（n）。
对信号x6（t），频率分辨率的选择要以能分辨开其中的三个频率对应的谱线为准则。
对周期序列，最好截取周期的整数倍进行分析，否则有可能产生较大的分析误差。
请实验者根据DFT的隐含周期性思考这个问题。
（4）编写主程序。
下图给出了主程序框图，供参考。
本实验提供FFT子程序和通用绘图子程序。
三．实验程序及对应波形
（1）对2中所给出的信号逐个进行谱分析。
实验代码如下（见后附程序Ex2_1）
%信号x1
x1=[1,1,1,1]
n1=[0:
3]
figure
（1）;subplot（3,1,1）
stem（n1,x1）
title（'信号x1时域特性'）
w=linspace（-5*pi,5*pi,1000）
X1=x1*exp（-j*n1'*w）
figure
（1）
subplot（3,1,2）
plot（w/pi,abs（X1））
title（'DTFT分析信号x1频域特性'）
X11=fftshift（fft（x1,128））
figure
（1）
subplot（3,1,3）
stem（abs（X11））
title（'fft分析信号x1频域特性'）
%信号x2
x2=[1,2,3,4,4,3,2,1]
n2=[0:
7]
figure
（2）;subplot（3,1,1）
stem（n2,x2）
title（'信号x2时域特性'）
w=linspace（-5*pi,5*pi,1000）
X2=x2*exp（-j*n2'*w）
figure
（2）
subplot（3,1,2）
plot（w/pi,abs（X2））
title（'DTFT分析信号x2频域特性'）
X22=fftshift（fft（x2,128））
figure
（2）
subplot（3,1,3）
stem（abs（X22））
title（'fft分析信号x2频域特性'）
%信号x3
x3=[43211234]
n3=[0:
7]
figure（3）;subplot（3,1,1）
stem（n3,x3）
title（'信号x3时域特性'）
w=linspace（-5*pi,5*pi,1000）
X3=x3*exp（-j*n3'*w）
figure（3）
subplot（3,1,2）
plot（w/pi,abs（X3））
title（'DTFT信号x3频域特性'）
X33=fftshift（fft（x3,128））
figure（3）
subplot（3,1,3）
stem（abs（X3））
title（'fft分析信号x3频域特性'）
%信号x4=cos（pi/4*n）
n4=[0:
15]
x4=cos（pi/4*n4）
figure（4）;subplot（3,1,1）
stem（n4,x4）
title（'信号x4时域特性'）
w=linspace（-5*pi,5*pi,1000）
X4=x4*exp（-j*n4'*w）
figure（4）
subplot（3,1,2）
plot（w/pi,abs（X4））
title（'DTFT信号x4频域特性'）
X44=fftshift（fft（x4,128））
figure（4）
subplot（3,1,3）
stem（abs（X44））
title（'fft分析信号x4频域特性'）
%信号x5=sin（pi/8*n）
n5=[0:
15]
x5=sin（pi/8*n5）
figure（5）;subplot（3,1,1）
stem（n5,x5）
title（'信号x5时域特性'）
w=linspace（-5*pi,5*pi,1000）
X5=x5*exp（-j*n5'*w）
figure（5）
subplot（3,1,2）
plot（w/pi,abs（X5））
X55=fftshift（fft（x5,128））
figure（5）
subplot（3,1,3）
stem（abs（X55））
title（'fft分析信号x4频域特性'）
%信号x6=cos（8*pi*t）+cos（16*pi*t）+cos（20*pi*t）
t=0:
0.0001:
2.5
xa=cos（8*pi*t）+cos（16*pi*t）+cos（20*pi*t）%原始的连续时间信号xa
figure（6）;subplot（3,1,1）;
plot（t,xa）;title（'连续时间信号曲线xa'）;;axis（[02.5-55]）%画连续时间信号曲线
T=0.01;n6=0:
2.5/T;%建立离散自变量向量
x6=cos（8*pi*n6*T）+cos（16*pi*n6*T）+cos（20*pi*n6*T）%采样周期为T的离散时间信号x
figure（6）;subplot（3,1,2）;
stem（n6*T,x6）;title（'取样信号'）;axis（[02.5-55]）%画离散时间信号曲线
w=linspace（-5*pi,5*pi,2500）;%均匀取点
X6=x6*exp（-j*n6'*w）%计算DTFT
figure（6）;subplot（3,1,3）;
plot（w/pi,abs（X6））;title（'DTFT分析取样信号频谱'）;
仿真结果如下：
（2）令x（n）=x4（n）+x5（n），用FFT计算8点和16点离散傅里叶变换，X（k）=DFT［x（n）］
（3）令x（n）=x4（n）+jx5（n），重复
（2）。
代码如下：
（见后附程序Ex2_2）
%信号x4=cos（pi/4*n）
n4=[0:
15]
x4=cos（pi/4*n4）
figure
（1）;subplot（2,1,1）
stem（n4,x4）
title（'信号x4时域特性'）
w=linspace（-5*pi,5*pi,1000）
X4=x4*exp（-j*n4'*w）
figure
（1）
subplot（2,1,2）
plot（w/pi,abs（X4））
title（'DTFT分析信号x4频域特性'）
%信号x5=sin（pi/8*n）
n5=[0:
15]
x5=sin（pi/8*n5）
figure
（2）;subplot（2,1,1）
stem（n5,x5）
title（'DTFT分析信号x5时域特性'）
w=linspace（-5*pi,5*pi,1000）
X5=x5*exp（-j*n5'*w）
figure
（2）
subplot（2,1,2）
plot（w/pi,abs（X5））
%x（n）=x4（n）+x5（n）
xn1=x4+x5
nx1=0:
length（xn1）-1
figure（3）;subplot（3,1,1）
stem（nx1,xn1）
title（'xn1=x4+x5时域特性'）
Xn1_8=fftshift（fft（xn1,8））
n1_8=0:
7
figure（3）;subplot（3,1,2）
stem（n1_8,abs（Xn1_8））
title（'xn1的8点fft变换'）
Xn1_16=fftshift（fft（xn1,16））
n1_16=0:
15
figure（3）;subplot（3,1,3）
stem（n1_16,abs（Xn1_16））
title（'xn1的16点fft变换'）
%x（n）=x4（n）+jx5（n）
xn2=x4+j*x5
nx2=0:
length（xn2）-1
figure（4）;subplot（3,1,1）
stem（nx2,xn2）
title（'xn2=x4+j*x5时域特性'）
Xn2_8=fftshift（fft（xn2,8））
n2_8=0:
7
figure（4）;subplot（3,1,2）
stem（n2_8,abs（Xn2_8））
title（'xn2的8点fft变换'）
Xn2_16=fftshift（fft（xn2,16））
n2_16=0:
15
figure（4）;subplot（3,1,3）
stem（n2_16,abs（Xn2_16））
title（'xn2的16点fft变换'）
仿真结果如下：
四、思考题
（1）在N=8时，x2（n）和x3（n）的幅频特性会相同吗?
为什么?
N=16呢?
答：
在N=8时，x2（n）和x3（n）的幅频特性会相同在N=8时x2（n）和x
升级会员 