数字信号处理上机实验答案第三版.docx

数字信号处理上机实验答案第三版.docx

《数字信号处理上机实验答案第三版.docx》由会员分享，可在线阅读，更多相关《数字信号处理上机实验答案第三版.docx（27页珍藏版）》请在冰豆网上搜索。

数字信号处理上机实验答案第三版

实验1：

系统响应及系统稳定性

实验程序清单:

closeall;clearall

%======内容1：

调用filter解差分方程，由系统对u（n）的响应判断稳定性======

A=[1,-0.9];B=[0.05,0.05];%系统差分方程系数向量B和A

x1n=[11111111zeros（1,50）];%产生信号x1（n）=R8（n）

x2n=ones（1,128）;%产生信号x2（n）=u（n）

hn=impz（B,A,58）;%求系统单位脉冲响应h（n）

subplot（2,2,1）;y='h（n）';stem（hn,'y'）;%调用函数tstem绘图

title（'（a）系统单位脉冲响应h（n）'）;

y1n=filter（B,A,x1n）;%求系统对x1（n）的响应y1（n）

subplot（2,2,2）;y='y1（n）';stem（y1n,'y'）;

title（'（b）系统对R8（n）的响应y1（n）'）;

y2n=filter（B,A,x2n）;%求系统对x2（n）的响应y2（n）

subplot（2,2,4）;y='y2（n）';stem（y2n,'y'）;

title（'（c）系统对u（n）的响应y2（n）'）;

%===内容2：

调用conv函数计算卷积============================

x1n=[11111111];%产生信号x1（n）=R8（n）

h1n=[ones（1,10）zeros（1,10）];

h2n=[12.52.51zeros（1,10）];

y21n=conv（h1n,x1n）;

y22n=conv（h2n,x1n）;

figure

（2）

subplot（2,2,1）;y='h1（n）';stem（h1n,'y'）;%调用函数tstem绘图

title（'（d）系统单位脉冲响应h1（n）'）;

subplot（2,2,2）;y='y21（n）';stem（y21n,'y'）;

title（'（e）h1（n）与R8（n）的卷积y21（n）'）;

subplot（2,2,3）;y='h2（n）';stem（h2n,'y'）;%调用函数tstem绘图

title（'（f）系统单位脉冲响应h2（n）'）;

subplot（2,2,4）;y='y22（n）';stem（y22n,'y'）;

title（'（g）h2（n）与R8（n）的卷积y22（n）'）;

%=========内容3：

谐振器分析========================

un=ones（1,256）;%产生信号u（n）

n=0:

255;

xsin=sin（0.014*n）+sin（0.4*n）;%产生正弦信号

A=[1,-1.8237,0.9801];B=[1/100.49,0,-1/100.49];%系统差分方程系数向量B和A

y31n=filter（B,A,un）;%谐振器对u（n）的响应y31（n）

y32n=filter（B,A,xsin）;%谐振器对u（n）的响应y31（n）

figure（3）

subplot（2,1,1）;y='y31（n）';stem（y31n,'y'）;

title（'（h）谐振器对u（n）的响应y31（n）'）;

subplot（2,1,2）;y='y32（n）';stem（y32n,'y'）;

title（'（i）谐振器对正弦信号的响应y32（n）'）;

实验程序运行结果及分析讨论

程序运行结果如图10.1.1所示。

实验内容

（2）系统的单位冲响应、系统对和的响应序列分别如图（a）、（b）和（c）所示；

实验内容（3）系统h1（n）和h2（n）对的输出响应分别如图（e）和（g）所示；

实验内容（4）系统对和的响应序列分别如图（h）和（i）所示。

由图（h）可见，系统对的响应逐渐衰减到零，所以系统稳定。

由图（i）可见，系统对的稳态响应近似为正弦序列，这一结论验证了该系统的谐振频率是0.4rad。

简答思考题

（1）如果输入信号为无限长序列，系统的单位脉冲响应是有限长序列，可否用线性卷积法求系统的响应。

①对输入信号序列分段；②求单位脉冲响应h（n）与各段的卷积；③将各段卷积结果相加。

具体实现方法有第三章介绍的重叠相加法和重叠保留法。

（2）如果信号经过低通滤波器，把信号的高频分量滤掉，时域信号的剧烈变化将被平滑，由实验内容

（1）结果图10.1.1（a）、（b）和（c）可见，经过系统低通滤波使输入信号、和的阶跃变化变得缓慢上升与下降。

实验二时域采样与频域采样（注：

本实验程序来自互联网，前半部分运行有误，请同学们自行检察，运行截图是正确的，可作参考）

实验程序清单:

1时域采样理论的验证程序清单

%时域采样理论验证程序exp2a.m

Tp=64/1000;%观察时间Tp=64微秒

%产生M长采样序列x（n）

%Fs=1000;T=1/Fs;

Fs=1000;T=1/Fs;

M=Tp*Fs;n=0:

M-1;

A=444.128;alph=pi*50*2^0.5;omega=pi*50*2^0.5;

xnt=A*exp（-alph*n*T）.*sin（omega*n*T）;

Xk=T*fft（xnt,M）;%M点FFT[xnt）]

yn='xa（nT）';subplot（3,2,1）;

stem（xnt,yn）;%调用自编绘图函数stem绘制序列图

boxon;title（'（a）Fs=1000Hz'）;

k=0:

M-1;fk=k/Tp;

subplot（3,2,2）;plot（fk,abs（Xk））;title（'（a）T*FT[xa（nT）],Fs=1000Hz'）;

xlabel（'f（Hz）'）;ylabel（'幅度'）;axis（[0,Fs,0,1.2*max（abs（Xk））]）

%=================================================

%Fs=300Hz和Fs=200Hz的程序与上面Fs=1000Hz完全相同。

2频域采样理论的验证程序清单

%频域采样理论验证程序exp2b.m

M=27;N=32;n=0:

M;

%产生M长三角波序列x（n）

xa=0:

floor（M/2）;xb=ceil（M/2）-1:

-1:

0;xn=[xa,xb];

Xk=fft（xn,1024）;%1024点FFT[x（n）],用于近似序列x（n）的TF

X32k=fft（xn,32）;%32点FFT[x（n）]

x32n=ifft（X32k）;%32点IFFT[X32（k）]得到x32（n）

X16k=X32k（1:

2:

N）;%隔点抽取X32k得到X16（K）

x16n=ifft（X16k,N/2）;%16点IFFT[X16（k）]得到x16（n）

subplot（3,2,2）;stem（n,xn,'.'）;boxon

title（'（b）三角波序列x（n）'）;xlabel（'n'）;ylabel（'x（n）'）;axis（[0,32,0,20]）

k=0:

1023;wk=2*k/1024;%

subplot（3,2,1）;plot（wk,abs（Xk））;title（'（a）FT[x（n）]'）;

xlabel（'\omega/\pi'）;ylabel（'|X（e^j^\omega）|'）;axis（[0,1,0,200]）

k=0:

N/2-1;

subplot（3,2,3）;stem（k,abs（X16k）,'.'）;boxon

title（'（c）16点频域采样'）;xlabel（'k'）;ylabel（'|X_1_6（k）|'）;axis（[0,8,0,200]）

n1=0:

N/2-1;

subplot（3,2,4）;stem（n1,x16n,'.'）;boxon

title（'（d）16点IDFT[X_1_6（k）]'）;xlabel（'n'）;ylabel（'x_1_6（n）'）;axis（[0,32,0,20]）

k=0:

N-1;

subplot（3,2,5）;stem（k,abs（X32k）,'.'）;boxon

title（'（e）32点频域采样'）;xlabel（'k'）;ylabel（'|X_3_2（k）|'）;axis（[0,16,0,200]）

n1=0:

N-1;

subplot（3,2,6）;stem（n1,x32n,'.'）;boxon

title（'（f）32点IDFT[X_3_2（k）]'）;xlabel（'n'）;ylabel（'x_3_2（n）'）;axis（[0,32,0,20]）

实验程序运行结果

1时域采样理论的验证程序运行结果exp2a.m如图10.3.2所示。

由图可见，采样序列的频谱的确是以采样频率为周期对模拟信号频谱的周期延拓。

当采样频率为1000Hz时频谱混叠很小；当采样频率为300Hz时，在折叠频率150Hz附近频谱混叠很严重；当采样频率为200Hz时，在折叠频率110Hz附近频谱混叠更很严重。

2时域采样理论的验证程序exp2b.m运行结果如图10.3.3所示。

该图验证了频域采样理论和频域采样定理。

对信号x（n）的频谱函数X（ejω）在[0，2π]上等间隔采样N=16时，N点IDFT[

]得到的序列正是原序列x（n）以16为周期进行周期延拓后的主值区序列：

由于N

与x（n）不相同，如图图10.3.3（c）和（d）所示。

当N=32时，如图图10.3.3（c）和（d）所示，由于N>M，频域采样定理，所以不存在时域混叠失真，因此。

与x（n）相同。

简答思考题

先对原序列x（n）以N为周期进行周期延拓后取主值区序列，

再计算N点DFT则得到N点频域采样：

实验三：

用FFT对信号作频谱分析

10.3.2实验程序清单

%第10章实验3程序exp3.m

%用FFT对信号作频谱分析

clearall;closeall

%实验内容

（1）===================================================

x1n=[ones（1,4）];%产生序列向量x1（n）=R4（n）

M=8;xa=1:

（M/2）;xb=（M/2）:

-1:

1;x2n=[xa,xb];%产生长度为8的三角波序列x2（n）

x3n=[xb,xa];

X1k8=fft（x1n,8）;%计算x1n的8点DFT

X1k16=fft（x1n,16）;%计算x1n的16点DFT

X2k8=fft（x2n,8）;%计算x1n的8点DFT

X2k16=fft（x2n,16）;%计算x1n的16点DFT

X3k8=fft（x3n,8）;%计算x1n的8点DFT

X3k16=fft（x3n,16）;%计算x1n的16点DFT

%以下绘制幅频特性曲线

subplot（2,2,1）;stem（X1k8）;%绘制8点DFT的幅频特性图

title（'（1a）8点DFT[x_1（n）]'）;xlabel（'ω/π'）;ylabel（'幅度'）;

axis（[0,2,0,1.2*max（abs（X1k8））]）

subplot（2,2,3）;stem（X1k16）;%绘制16点DFT的幅频特性图

title（'（1b）16点DFT[x_1（n）]'）;xlabel（'ω/π'）;ylabel（'幅度'）;

axis（[0,2,0,1.2*max（abs（X1k16））]）

figure

（2）

subplot（2,2,1）;stem（X2k8）;%绘制8点DFT的幅频特性图

title（'（2a）8点DFT[x_2（n）]'）;xlabel（'ω/π'）;ylabel（'幅度'）;

axis（[0,2,0,1.2*max（abs（X2k8））]）

subplot（2,2,2）;stem（X2k16）;%绘制16点DFT的幅频特性图

title（'（2b）16点DFT[x_2（n）]'）;xlabel（'ω/π'）;ylabel（'幅度'）;

axis（[0,2,0,1.2*max（abs（X2k16））]）

subplot（2,2,3）;stem（X3k8）;%绘制8点DFT的幅频特性图

title（'（3a）8点DFT[x_3（n）]'）;xlabel（'ω/π'）;ylabel（'幅度'）;

axis（[0,2,0,1.2*max（abs（X3k8））]）

subplot（2,2,4）;stem（X3k16）;%绘制16点DFT的幅频特性图

title（'（3b）16点DFT[x_3（n）]'）;xlabel（'ω/π'）;ylabel（'幅度'）;

axis（[0,2,0,1.2*max（abs（X3k16））]）

%实验内容

（2）周期序列谱分析==================================

N=8;n=0:

N-1;%FFT的变换区间N=8

x4n=cos（pi*n/4）;

x5n=cos（pi*n/4）+cos（pi*n/8）;

X4k8=fft（x4n）;%计算x4n的8点DFT

X5k8=fft（x5n）;%计算x5n的8点DFT

N=16;n=0:

N-1;%FFT的变换区间N=16

x4n=cos（pi*n/4）;

x5n=cos（pi*n/4）+cos（pi*n/8）;

X4k16=fft（x4n）;%计算x4n的16点DFT

X5k16=fft（x5n）;%计算x5n的16点DFT

figure（3）

subplot（2,2,1）;stem（X4k8）;%绘制8点DFT的幅频特性图

title（'（4a）8点DFT[x_4（n）]'）;xlabel（'ω/π'）;ylabel（'幅度'）;

axis（[0,2,0,1.2*max（abs（X4k8））]）

subplot（2,2,3）;stem（X4k16）;%绘制16点DFT的幅频特性图

title（'（4b）16点DFT[x_4（n）]'）;xlabel（'ω/π'）;ylabel（'幅度'）;

axis（[0,2,0,1.2*max（abs（X4k16））]）

subplot（2,2,2）;stem（X5k8）;%绘制8点DFT的幅频特性图

title（'（5a）8点DFT[x_5（n）]'）;xlabel（'ω/π'）;ylabel（'幅度'）;

axis（[0,2,0,1.2*max（abs（X5k8））]）

subplot（2,2,4）;stem（X5k16）;%绘制16点DFT的幅频特性图

title（'（5b）16点DFT[x_5（n）]'）;xlabel（'ω/π'）;ylabel（'幅度'）;

axis（[0,2,0,1.2*max（abs（X5k16））]）

%实验内容（3）模拟周期信号谱分析===============================

figure（4）

Fs=64;T=1/Fs;

N=16;n=0:

N-1;%FFT的变换区间N=16

x6nT=cos（8*pi*n*T）+cos（16*pi*n*T）+cos（20*pi*n*T）;%对x6（t）16点采样

X6k16=fft（x6nT）;%计算x6nT的16点DFT

X6k16=fftshift（X6k16）;%将零频率移到频谱中心

Tp=N*T;F=1/Tp;%频率分辨率F

k=-N/2:

N/2-1;fk=k*F;%产生16点DFT对应的采样点频率（以零频率为中心）

subplot（3,1,1）;stem（fk,abs（X6k16）,'.'）;boxon%绘制8点DFT的幅频特性图

title（'（6a）16点|DFT[x_6（nT）]|'）;xlabel（'f（Hz）'）;ylabel（'幅度'）;

axis（[-N*F/2-1,N*F/2-1,0,1.2*max（abs（X6k16））]）

N=32;n=0:

N-1;%FFT的变换区间N=16

x6nT=cos（8*pi*n*T）+cos（16*pi*n*T）+cos（20*pi*n*T）;%对x6（t）32点采样

X6k32=fft（x6nT）;%计算x6nT的32点DFT

X6k32=fftshift（X6k32）;%将零频率移到频谱中心

Tp=N*T;F=1/Tp;%频率分辨率F

k=-N/2:

N/2-1;fk=k*F;%产生16点DFT对应的采样点频率（以零频率为中心）

subplot（3,1,2）;stem（fk,abs（X6k32）,'.'）;boxon%绘制8点DFT的幅频特性图

title（'（6b）32点|DFT[x_6（nT）]|'）;xlabel（'f（Hz）'）;ylabel（'幅度'）;

axis（[-N*F/2-1,N*F/2-1,0,1.2*max（abs（X6k32））]）

N=64;n=0:

N-1;%FFT的变换区间N=16

x6nT=cos（8*pi*n*T）+cos（16*pi*n*T）+cos（20*pi*n*T）;%对x6（t）64点采样

X6k64=fft（x6nT）;%计算x6nT的64点DFT

X6k64=fftshift（X6k64）;%将零频率移到频谱中心

Tp=N*T;F=1/Tp;%频率分辨率F

k=-N/2:

N/2-1;fk=k*F;%产生16点DFT对应的采样点频率（以零频率为中心）

subplot（3,1,3）;stem（fk,abs（X6k64）,'.'）;boxon%绘制8点DFT的幅频特性图

title（'（6a）64点|DFT[x_6（nT）]|'）;xlabel（'f（Hz）'）;ylabel（'幅度'）;

axis（[-N*F/2-1,N*F/2-1,0,1.2*max（abs（X6k64））]）

程序运行结果分析讨论：

请读者注意，用DFT（或FFT）分析频谱，绘制频谱图时，最好将X（k）的自变量k换算成对应的频率，作为横坐标便于观察频谱。

为了便于读取频率值，最好关于π归一化，即以

作为横坐标。

1、实验内容

（1）

图（1a）和（1b）说明

的8点DFT和16点DFT分别是

的频谱函数的8点和16点采样；

因为

，所以，

与

的8点DFT的模相等，如图（2a）和（3a）。

但是，当N=16时，

与

不满足循环移位关系，所以图（2b）和（3b）的模不同。

2、实验内容

（2），对周期序列谱分析

的周期为8，所以N=8和N=16均是其周期的整数倍，得到正确的单一频率正弦波的频谱，仅在0.25π处有1根单一谱线。

如图（4b）和（4b）所示。

的周期为16，所以N=8不是其周期的整数倍，得到的频谱不正确，如图（5a）所示。

N=16是其一个周期，得到正确的频谱，仅在0.25π和0.125π处有2根单一谱线,如图（5b）所示。

3、实验内容（3），对模拟周期信号谱分析

有3个频率成分，

。

所以

的周期为0.5s。

采样频率

。

变换区间N=16时，观察时间Tp=16T=0.25s，不是

的整数倍周期，所以所得频谱不正确，如图（6a）所示。

变换区间N=32,64时，观察时间Tp=0.5s，1s，是

的整数周期，所以所得频谱正确，如图（6b）和（6c）所示。

图中3根谱线正好位于

处。

变换区间N=64时频谱幅度是变换区间N=32时2倍，这种结果正好验证了用DFT对中期序列谱分析的理论。

注意：

（1）用DFT（或FFT）对模拟信号分析频谱时，最好将X（k）的自变量k换算成对应的模拟频率fk，作为横坐标绘图，便于观察频谱。

这样，不管变换区间N取信号周期的几倍，画出的频谱图中有效离散谐波谱线所在的频率值不变，如图（6b）和（6c）所示。

（2）本程序直接画出采样序列N点DFT的模值，实际上分析频谱时最好画出归一化幅度谱，这样就避免了幅度值随变换区间N变化的缺点。

本实验程序这样绘图只要是为了验证了用DFT对中期序列谱分析的理论。

简答思考题

思考题

（1）和

（2）的答案请读者在教材3.？

节找，思考题（3）的答案在程序运行结果分析讨论已经详细回答。

实验四IIR数字滤波器设计及软件实现

2、实验程序清单

%实验4程序exp4.m

%IIR数字滤波器设计及软件实现

clearall;closeall

Fs=10000;T=1/Fs;%采样频率

%调用信号产生函数mstg产生由三路抑制载波调幅信号相加构成的复合信号st

st=mstg;

%低通滤波器设计与实现=========================================

fp=280;fs=450;

wp=2*fp/Fs;ws=2*fs/Fs;rp=0.1;rs=60;%DF指标（低通滤波器的通、阻带边界频）

[N,wp]=ellipord（wp,ws,rp,rs）;%调用ellipord计算椭圆DF阶数N和通带截止频率wp

[B,A]=ellip（N,rp,rs,wp）;%调用ellip计算椭圆带通DF系统函数系数向量B和A

y1t=filter（B,A,st）;%滤波器软件实现

%低通滤波器设计与实现绘图部分

figure

（2）;subplot（3,1,1）;

myplot（B,A）;%调用绘图函数myplot绘制损耗函数曲线

yt='y_1（t）';

subplot（3,1,2）;tplot（y1t,T,yt）;%调用绘图函数tplot绘制滤波器输出波形

%带通滤波器设计与实现====================================================

fpl=440;fpu=560;fsl=275;fsu=900;

wp=[2*fpl/Fs,2*fpu/Fs];ws=[2*fsl/Fs,2*fsu/Fs];rp=0.1;rs=60;

[N,wp]=ellipord（wp,ws,rp,rs）;%调用ellipord计算椭圆DF阶数N和通带截止频率wp

[B,A]=ellip（N,rp,rs,wp）;%调用ellip计算椭圆带通DF系统函数系数向量B和A

y2t=filter（B,A,st）;%滤波器软件实现

%带通滤波器设计与实现绘图部分（省略）

%高通滤波器设计与实现================================================

fp=890;fs=600;

wp=2*fp/Fs;ws=2*fs/Fs;rp=0.1;rs=60;%DF指标（低通滤波器的通、阻带边界频）

[N,wp]=ellipord（wp,ws,rp,rs）;%调用ellipord计算椭圆DF阶数N和通带截止频率wp

[B,A]=ellip（N,rp,rs,wp,'high'）;%调用ellip计算椭圆带通DF系统函数系数向量B和A

y3t=filter（B,A,st）;%滤波器软件实现

%高低通滤波器设计与实现绘图部分（省略）

10.4.3实验程序