数字信号处理malab上机实验含子程序.docx
《数字信号处理malab上机实验含子程序.docx》由会员分享，可在线阅读，更多相关《数字信号处理malab上机实验含子程序.docx（26页珍藏版）》请在冰点文库上搜索。
数字信号处理malab上机实验含子程序
数字信号处理上机内容
实验一、求线性时不变系统的输出
一、实验目的：
学习用线性卷积法求网络输出的方法。
二、实验原理：
一般网络或系统用线性常系数差分方程描述，如果已知差分方程和输入信号，用递推法求解差分方程或者求网络输出，最适合用计算计求解。
一、实验结果及程序清单（含分析及函数用法）
1、%-----
（2）调用filter解差分方程以及单位脉冲响应------
closeall;clearall
A=[1,-0.9];B=[0.05,0.05];%系统差分方程系数向量B和A
x1n=[ones（1,8）zeros（1,25）];%产生信号x1（n）=R8（n），用zeros用来加点的个数
x2n=ones（1,30）;%产生信号x2（n）=u（n）
hn=impz（B,A,25）;%求系统单位脉冲响应h（n）
subplot（3,1,1）;stem（hn）;%调用函数stem绘图
title（'（a）系统单位脉冲响应h（n）'）;
y1n=filter（B,A,x1n）;%求系统对x1（n）的响应y1（n）
subplot（3,1,2）;stem（y1n）;
title（'（b）系统对R8（n）的响应y1（n）'）;
y2n=filter（B,A,x2n）;%求系统对x2（n）的响应y2（n）
subplot（3,1,3）;stem（y2n）;
title（'（c）系统对u（n）的响应y2（n）'）;
分析：
1、在时域求系统响应的方法有两种，第一种是通过解差分方程求得系统输出；第二种是已知系统的单位脉冲响应，通过求输入信号和系统单位脉冲响应的线性卷积求得系统输出。
2、检验系统的稳定性，其方法是在输入端加入单位阶跃序列，观察输出波形，如果波形稳定在一个常数值上，系统稳定，否则不稳定。
（a）25个点数和程序所写一致。
Filter函数实现线性常系数差分方程的递推求解。
调用格式如下：
Y=[filter（B,A,x）]计算系统对输入信号x的零状态响应输出信号向量Y，
B,A是差分方程的系数向量。
即
B=[a1,a2……am]A=[b1,b2……bn]
2、%-----（3）调用conv函数计算卷积-------
x1n=[ones（1,8）];%产生信号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）;
subplot（2,2,1）;stem（h1n）;
title（'（d）系统单位脉冲响应h1（n）'）;
subplot（2,2,3）;stem（y21n）;
title（'（e）h1（n）与R8（n）的卷积y21（n）'）;
subplot（2,2,2）;stem（h2n）;
title（'（f）系统单位脉冲响应h2（n）'）;
subplot（2,2,4）;stem（y22n）;
title（'（g）h2（n）与R8（n）的卷积y22（n）'）;
分析：
（d）（f）单位脉冲响应点数与程序要求一致
（e）（g）卷积点数满足M+N-1的要求，图形也满足要求。
Conv函数用于计算两个有限长序列的卷积
C=conv（A,B）计算两个有限长序列向量A和B的卷积
3、%-----（4）实验方法检查系统是否稳定-------
closeall;clearall
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
y1n=filter（B,A,un）;%谐振器对u（n）的响应y31（n）
y2n=filter（B,A,xsin）;%谐振器对u（n）的响应y31（n）
subplot（2,1,1）;stem（y1n）;
title（'（h）谐振器对u（n）的响应y31（n）'）;
subplot（2,1,2）;stem（y2n）;
title（'（i）谐振器对正弦信号的响应y32（n）'）;
分析：
（h）中由：
在系统的输入端加入单位阶跃序列，如果系统的输出趋近一个常数（包括零），就可以断定系统是稳定的
输出显然趋近于零，所以是稳定的
（i）中谐振器具有对某个频率进行谐振的性质，本实验中的谐振器的谐振频率是0.4rad,因此稳定波形为sin（0.4n）。
二、思考题简要分析
（1）如果输入信号为无限长序列，系统的单位脉冲响应是有限长序列，可否用线性卷积法求系统的响应?
如何求？
答：
如果输入信号为无限长序列，系统的单位脉冲响应是有限长序列，可用分段线性卷积法求系统的响应
（2）如果信号经过低通滤波器，把信号的高频分量滤掉，时域信号会有何变化，用前面第一个实验结果进行分析说明。
答：
如果信号经过低通滤波器，则信号的高频分量被滤掉，时域信号的变化减缓，在有阶跃处附近产生过渡带。
因此，当输入矩形序列时，输出序列的开始和终了都产生了明显的过渡带，见第一个实验结果的波形。
三、总结
心里总结：
起初没有看书，什么都不懂，做实验的时候也是什么都不懂，混时间，非常的痛苦。
说要认真的写实验报告，也是笑话，还要分析，更是痛苦。
这次发时间把书本的第一章好好的看了几遍，包括原理性的东西，以及其中涉及的matlab的代码全看了。
再则就是参考了matlab基础书籍以及网上的一些资料，最重要的是matlab的帮助文档作用很大，函数怎么用的都写好了。
把实验一完成后，觉得也没有什么呀，当初那么的害怕做什么呢！
！
！
！
总结即在实验原理中说明的两点：
1、在时域求系统响应的方法有两种，第一种是通过解差分方程求得系统输出；第二种是已知系统的单位脉冲响应，通过求输入信号和系统单位脉冲响应的线性卷积求得系统输出。
2、检验系统的稳定性，其方法是在输入端加入单位阶跃序列，观察输出波形，如果波形稳定在一个常数值上，系统稳定，否则不稳定。
实验二：
时域采样与频域采样
一、实验原理与方法
1、时域采样定理：
a）对模拟信号
以间隔T进行时域等间隔理想采样，形成的采样信号的频谱
是原模拟信号频谱
以采样角频率
（
）为周期进行周期延拓。
公式为：
b）采样频率
必须大于等于模拟信号最高频率的两倍以上，才能使采样信号的
频谱不产生频谱混叠。
2、频域采样定理：
公式为：
由公式可知，频域采样点数N必须大于等于时域离散信号的长度M（即N≥M），才能使时域不产生混叠，则N点IDFT[
]得到的序列
就是原序列x（n）,即
=x（n）。
二、实验结果及程序清单（含分析及函数用法）
1、%---------
（1）时域采样理论验证程序---------
Tp=64/1000;%观察时间Tp=64毫秒
%产生M长采样序列x（n）
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）;%调用自编绘图函数stem绘制序列图
boxon;title（'（a）Fs=1000Hz'）;
k=0:
M-1;fk=k/Tp;
subplot（3,2,2）;stem（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完全相同。
Tp=64/1000;%观察时间Tp=64毫秒
%产生M长采样序列x（n）
Fs=300;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,3）;
stem（xnt）;%调用自编绘图函数stem绘制序列图
boxon;title（'（b）Fs=300Hz'）;
k=0:
M-1;fk=k/Tp;
subplot（3,2,4）;stem（fk,abs（Xk））;title（'（b）T*FT[xa（nT）],Fs=300Hz'）;
xlabel（'f（Hz）'）;ylabel（'幅度'）;axis（[0,Fs,0,1.2*max（abs（Xk））]）;
Tp=64/1000;%观察时间Tp=64毫秒
%产生M长采样序列x（n）
Fs=200;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,5）;
stem（xnt）;%调用自编绘图函数stem绘制序列图
boxon;title（'（c）Fs=200Hz'）;
k=0:
M-1;fk=k/Tp;
subplot（3,2,6）;stem（fk,abs（Xk））;title（'（c）T*FT[xa（nT）],Fs=200Hz'）;
xlabel（'f（Hz）'）;ylabel（'幅度'）;axis（[0,Fs,0,1.2*max（abs（Xk））]）;
分析：
，当采样频率为1000Hz时，频谱混叠很小：
当采样频率为300Hz时，频谱混叠很严重；当采样频率为200Hz时，频谱混叠更很严重。
Fft函数的调用格式：
Xk=fft（xn,N）
调用参数xn为被交换的时域序列向量，N是DFT变换的区间长度，当N大于xn的长度时，fft函数自动在xn后面补零。
当N小于xn的长度时，fft函数计算xn的前面N个元素构成的N长序列的N点DFT,忽略xn后面的元素。
2、%-------------
（2）频域采样理论验证-----------
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]）
分析;该图验证了频域采样理论和频域采样定理。
对信号x（n）的频谱函数X（ejω）在［0，2π］上等间隔采样N=16时，N点IDFT［XN（k）］得到的序列正是原序列x（n）以16为周期进行周期延拓后的主值区序列
当N=16时，由于N当N=32时，由于N>M，频域采样定理，所以不存在时域混叠失真
Ifft函数用法同fft函数
三、思考题简要分析
1、如果序列x（n）的长度为M，希望得到其频谱
在
上的N点等间隔采样，当N答：
先对原序列x（n）以N为周期进行周期延拓后取主值区序列，
再计算N点DFT则得到N点频域采样：
实验三：
用FFT对信号作频谱分析
一、实验原理与方法
1、周期信号的频谱是离散谱，只有用整数倍周期的长度作FFT，得到的离散谱才能代表周期信号的频谱。
如果不知道信号周期，可以尽量选择信号的观察时间长一些。
2、对模拟信号进行谱分析时，首先要按照采样定理将其变成时域离散信号。
如果是模拟周期信号，也应该选取整数倍周期的长度，经过采样后形成周期序列，按照周期序列的谱分析进行。
3、频谱分辨率直接和FFT的变换区间N有关，因为FFT能够实现的频率分辨率是
，因此要求
。
二、实验结果及程序清单（含分析及函数用法）
1、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（3,2,1）;
mstem（X1k8）;%绘制8点DFT的幅频特性图
title（'（1a）8点DFT[x_1（n）]'）;xlabel（'ω/π'）;ylabel（'幅度'）;
axis（[0,2,0,1.2*max（abs（X1k8））]）
subplot（3,2,2）;mstem（X1k16）;%绘制16点DFT的幅频特性图
title（'（1b）16点DFT[x_1（n）]'）;xlabel（'ω/π'）;ylabel（'幅度'）;
axis（[0,2,0,1.2*max（abs（X1k16））]）
subplot（3,2,3）;mstem（X2k8）;%绘制8点DFT的幅频特性图
title（'（2a）8点DFT[x_2（n）]'）;xlabel（'ω/π'）;ylabel（'幅度'）;
axis（[0,2,0,1.2*max（abs（X2k8））]）
subplot（3,2,4）;mstem（X2k16）;%绘制16点DFT的幅频特性图
title（'（2b）16点DFT[x_2（n）]'）;xlabel（'ω/π'）;ylabel（'幅度'）;
axis（[0,2,0,1.2*max（abs（X2k16））]）
subplot（3,2,5）;mstem（X3k8）;%绘制8点DFT的幅频特性图
title（'（3a）8点DFT[x_3（n）]'）;xlabel（'ω/π'）;ylabel（'幅度'）;
axis（[0,2,0,1.2*max（abs（X3k8））]）
subplot（3,2,6）;mstem（X3k16）;%绘制16点DFT的幅频特性图
title（'（3b）16点DFT[x_3（n）]'）;xlabel（'ω/π'）;ylabel（'幅度'）;
axis（[0,2,0,1.2*max（abs（X3k16））]）
2、%-------------实验内容
（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
subplot（2,2,1）;mstem（X4k8）;%绘制8点DFT的幅频特性图
title（'（a）8点DFT[x_4（n）]'）;xlabel（'ω/π'）;ylabel（'幅度'）;
axis（[0,2,0,1.2*max（abs（X4k8））]）
subplot（2,2,3）;mstem（X4k16）;%绘制16点DFT的幅频特性图
title（'（b）16点DFT[x_4（n）]'）;xlabel（'ω/π'）;ylabel（'幅度'）;
axis（[0,2,0,1.2*max（abs（X4k16））]）
subplot（2,2,2）;mstem（X5k8）;%绘制8点DFT的幅频特性图
title（'（a）8点DFT[x_5（n）]'）;xlabel（'ω/π'）;ylabel（'幅度'）;
axis（[0,2,0,1.2*max（abs（X5k8））]）
subplot（2,2,4）;mstem（X5k16）;%绘制16点DFT的幅频特性图
title（'（b）16点DFT[x_5（n）]'）;xlabel（'ω/π'）;ylabel（'幅度'）;
axis（[0,2,0,1.2*max（abs（X5k16））]）
分析：
的周期为8，所以N=8和N=16均是其周期的整数倍，得到正确的单一频率正弦波的频谱，仅在0.25π处有1根单一谱线。
如图（a）和（b）所示。
的周期为16，所以N=8不是其周期的整数倍，得到的频谱不正确，如图（a）所示。
N=16是其一个周期，得到正确的频谱，仅在0.25π和0.125π处有2根单一谱线,如图（b）所示。
3、%----------实验内容（3）---------
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））]）
分析：
x6（t）有3个频率成分，f1=4Hz,f2=8Hz,f3=10Hz。
所以x6（t）的周期为0.5s。
采样频率Fs=64Hz=16f1=8f2=6.4f3。
变换区间N=16时，观察时间Tp=16T=0.25s，不是x6（t）的整数倍周期，所以所得频谱不正确，如图8.3.1（6a）所示。
变换区间N=32，64时，观察时间Tp=0.5s，1s，
升级会员 