实验一时域离散信号与及系统变换域分析.docx
《实验一时域离散信号与及系统变换域分析.docx》由会员分享，可在线阅读，更多相关《实验一时域离散信号与及系统变换域分析.docx（19页珍藏版）》请在冰豆网上搜索。
实验一时域离散信号与及系统变换域分析
实验一时域离散信号与及系统变换域分析
一、实验目的
1．了解时域离散信号的产生及基本运算实现。
2．掌握离散时间傅里叶变换实现及系统分析方法。
3.熟悉离散时间傅里叶变换性质。
4.掌握系统Z域分析方法。
5.培养学生运用软件分析、处理数字信号的能力。
二、实验设备
1、计算机
2、Matlab7.0以上版本
三、实验内容
1、对于给定的时域离散信号会进行频谱分析，即序列的傅里叶变换及其性质分析。
2、对于离散系统会进行频域分析及Z域分析。
包括频谱特性、零极点画图、稳定性分析。
3、对于差分方程会用程序求解，包括求单位冲击序列响应，零输入响应、零状态响应、全响应，求其系统函数，及其分析。
4、信号时域采样及其频谱分析，序列恢复。
5、扩展部分主要是关于语音信号的读取及其播放。
四、实验原理
1、序列的产生及运算
在Matlab中自带了cos、sin、exp（指数）等函数，利用这些函数可以产生实验所需序列。
序列的运算包括序列的加法、乘法，序列
的移位
，翻褶
等。
序列的加法或乘法指同序号的序列值逐项对应相加或相乘，但Matlab中“+”“.*”运算是对序列的值直接进行加或乘，不考虑两序列的序号是否相同，因此编程时考虑其序号的对应。
2、序列的傅里叶变换及其性质
序列的傅里叶变换定义：
，其幅度特性为
，在Matlab中采用abs函数；相位特性为
，在Matlab中采用angle函数。
序列傅里叶变换的性质：
（1）FT的周期性
，实序列傅里叶变换的对称性
。
对实序列和复序列分别进行傅里叶变换，通过图形结果观察周期性即对称性。
（2）FT的频移特性
，对序列在时域乘以
，然后进傅里叶变换，比较其结果和直接对序列进行傅里叶变换的不同。
（3）时域卷积定理：
若
，对序列
和
进行线性卷积得到
，分别对它们进行傅里叶变换，应满足
。
3、离散时间系统的Z域分析
已知离散时间系统的差分方程为
，对等号两边进行Z变换，得到其系统函数
及系统零极点，对系统函数进行反变换得到单位取样响应
，根据单位取样响应或系统函数的系数可以得到频率响应
，根据极点位置判断系统稳定性。
4、信号时域采样及恢复
给定连续信号
，对其用不同的采样频率进行采样，根据时域采样定理，采样信号的频谱是原模拟信号频谱沿频率轴以
为周期延拓而成的，并且要不失真地还原出模拟信号时，要满足
，因此当采样频率满足和不满足采样定理时，所得到的频谱是不同的。
根据采样信号进行信号恢复时，采用内插公式
实现。
五、实验步骤
1、序列的基本运算
1.1产生余弦信号
及带噪信号
0<=n<=50（噪声采用randn函数）
1.2已知
，
求两个序列的和、乘积、序列x1的移位序列（右移2位）,序列x2的翻褶序列，画出原序列及运算结果图。
2、序列的傅里叶变换
2.1已知序列
。
试求它的傅里叶变换，并且画出其幅度、相角、实部和虚部的波形，并分析其含有的频率分量主要位于高频区还是低频区。
2.2令
，求其傅立叶变换
。
分别用
和
对其进行采样，求出离散时间傅立叶变换
，画出相应频谱，分析结果的不同及原因。
3、序列的傅里叶变换性质分析
3.1已知序列
，
，求其傅里叶变换，并讨论其傅里叶变换的周期性和对称性。
3.2已知序列
，
，求其傅里叶变换，并讨论其傅里叶变换的周期性和对称性。
为了方便，考虑在两个周期，例如[
]中2M+1个均匀频率点上计算FT，并且观察其周期性和对称性。
为此给出function文件如下，求解FT变换：
function[X,w]=ft1（x,n,k）
w=（pi/abs（max（k）/2））*k
X=x*（exp（-j*pi/abs（max（k）/2）））.^（n'*k）
3.3编写程序验证序列傅里叶变换频移性质，时域卷积定理（时域卷积后的频域特性）。
（所需信号自行选择）
4、时域差分方程的求解
4.1求解差分方程y（n）＋a1y（n-1）＋a2y（n-2）=b0x（n）＋b1x（n-1）的零状态响应和全响应。
已知X（n）为单位取样序列，y（-1）=1,y（-2）=2,a1=0.5,a2=0.06,b0=2,b1＝3。
5、离散系统的Z域分析
5.1利用系统函数
分析系统的稳定性。
假设系统函数如下式：
，试判断系统是否稳定。
5.2已知线性时不变系统的系统函数
，编写程序求其单位取样响应，频率响应及系统零极点，并画出相应图形。
6、创新训练拓展内容
6.1利用Matlab自带的录音功能，或利用Goldwave等音频编辑软件，对语音或其他音频信号进行采集并保存为*.wav文件。
要求：
（1）采用不同的采样频率（2000Hz，4000Hz，8000Hz，16000Hz等）。
（2）对采集得到的信号进行播放，并画图。
（3）分析在不同采样频率下得到的信号有何不同。
6.2设定一个连续时间信号，进行抽样和恢复，要求分析不同采样频率对恢复结果的影响，给出实验程序及各关键步骤图形结果。
6.3设计内容
设计一个离散系统，给定系统函数或差分方程，设定激励及初始条件。
要求：
（1）绘制系统函数零极点图，判断稳定性；
（2）求单位脉冲响应h（n）；
（3）求系统零输入响应及零状态响应，要求零状态响应采样三种方法求解（卷积的方法、迭代解法、变换域求解方法），激励自定；
（4）分析系统频响特性，画出频响函数幅频曲线和相频曲线。
六、实验要求
第一部分：
验证实验内容
根据给定的实验内容，部分实验给出了参考程序段，见下面各段程序。
请基于Matlab环境进行验证实验。
第二部分：
编程实验内容
对于给定的实验内容中，没有参考程序段的部分，进行编程，给出实验结果，并进行相应的分析。
第三部分：
创新训练拓展内容
此部分内容，要求给出程序设计流程图（画法见附录3），给出程序内容的解释，并对结果进行分析。
7、思考题
下面四个二阶网络的系统函数具有一样的极点分布：
1）
2）
3）
4）
请分析研究零点分布对于单位脉冲响应的影响。
要求：
（1）分别画出各系统的零、极点分布图；
（2）分别求出各系统的单位脉冲响应，并画出其波形；
（3）分析零点分布对于单位脉冲响应的影响。
八、实验参考资料
1、高西全,丁玉美.数字信号处理[M].西安：
西安电子科技大学出版社,2008
2、张德丰.详解MATLAB数字信号处理[M].北京：
电子工业出版社,2010
3、王月明,张宝华.MATLAB基础与应用教程[M].北京:
北京大学出版社,2012
附：
实验报告要求：
实验名称：
-------
班级：
组号：
姓名1（学号）、姓名2（学号）、姓名3（学号）
一、实验目的（手写）
二、实验主要内容（要根据自己组所做内容写，做了的写，没做的不要写）
例如：
1.对序列的产生和运算方法进行实现
2.序列的傅里叶变换实现、性质及分析
等等
三、实验主要仪器、设备及软件（手写）
四、实验步骤、结果与分析
例如：
1.序列的运算
序列为……，进行加法、乘法、……运算
运算结果为……
2.序列的傅里叶变换实现及分析
（1）已知序列
。
试求它的傅里叶变换，并且画出其幅度、相角、实部和虚部的波形，并分析其含有的频率分量主要位于高频区还是低频区。
程序
结果
分析
（2）序列的傅里叶变换性质分析
……
注1：
（包括程序框图及代码、图形、数据等），其中程序框图、代码、图形可以直接打印，结果分析手写。
注2：
对已给出（程序、结果及分析）的验证性实验，自己运行即可，可以不用写在报告中。
对已给出（程序）的验证性实验，程序可以不用写在实验报告中，只写出结果和分析。
五、实验结论与总结（手写）
六、思考题（分析手写）
七、实验参考资料
附：
实验所需部分函数及验证性程序:
1、序列的基本运算
%1.单位取样序列x（n）=delta（n-n0）要求n1<=n0<=n2
function[x,n]=impseq（n0,n1,n2）
n=[n1:
n2];x=[（n-n0）==0];==是逻辑判断
%2.单位阶跃序列x（n）=u（n-n0）要求n1<=n0<=n2
function[x,n]=stepseq（n0,n1,n2）
n=[n1:
n2];x=[（n-n0）>=0];
%3.信号加y（n）=x1（n）+x2（n）
%find函数：
找出非零元素的索引号
%x1:
第一个序列的值，n1:
序列x1的索引号
%x2:
第二个序列的值，n2:
序列x2的索引号
function[y,n]=sigadd（x1,n1,x2,n2）
n=min（min（n1）,min（n2））:
max（max（n1）,max（n2））;
y1=zeros（1,length（n））;y2=y1;
y1（find（（n>=min（n1））&（n<=max（n1））==1））=x1;
y2（find（（n>=min（n2））&（n<=max（n2））==1））=x2;
y=y1+y2;
%4.信号乘y（n）=x1（n）*x2（n）
function[y,n]=sigmult（x1,n1,x2,n2）
n=min（min（n1）,min（n2））:
max（max（n1）,max（n2））;
y1=zeros（1,length（n））;y2=y1;
y1（find（（n>=min（n1））&（n<=max（n1））==1））=x1;
y2（find（（n>=min（n2））&（n<=max（n2））==1））=x2;
y=y1.*y2;
%5.移位y（n）=x（n-n0）
function[y,n]=sigshift（x,m,n0）
n=m+n0;y=x;
%6.翻褶y（n）=x（-n）
function[y,n]=sigfold（x,n）
y=fliplr（x）;n=-fliplr（n）;
2、序列的傅里叶变换
%7.求序列
的傅里叶变换
w=[0:
1:
500]*pi/500
X=exp（j*w）./（exp（j*w）-0.5*ones（1,501））
magX=abs（X）
angX=angle（X）
realX=real（X）
imagX=imag（X）
subplot（2,2,1）
plot（w/pi,magX）
grid
xlabel（'frequencyinpiunits'）
title（'MagnitudePart'）
ylabel（'Magnitude'）
subplot（2,2,3）
plot（w/pi,angX）
grid
xlabel（'frequencyinpiunits'）
title（'AnglePart'）
ylabel（'Radians'）
subplot（2,2,2）
plot（w/pi,realX）
grid
xlabel（'frequencyinpiunits'）
title（'RealPart'）
ylabel（'Real'）
subplot（2,2,4）
plot（w/pi,imagX）
grid
xlabel（'frequencyinpiunits'）
title（'ImaginaryPart'）
ylabel（'Imaginary'）
程序执行结果：
%8令
，绘制其傅立叶变换
。
用不同频率对其进行采样，分别画出
。
Dt=0.00005;%步长为0.00005s
t=-0.005:
Dt:
0.005;
xa=exp（-1000*abs（t））;%取时间从-0.005s到0.005s这段模拟信号
Wmax=2*pi*2000;%信号最高频率为2
*2000
K=500;%频域正半轴取500个点进行计算
k=0:
1:
K;
W=k*Wmax/K;%
求模拟角频率
Xa=xa*exp（-j*t'*W）*Dt;%计算连续时间傅立叶变换（利用矩阵运算实现）
Xa=real（Xa）;%取实部
W=[-fliplr（W）,W（2:
501）];%将角频率范围扩展为从-到+
Xa=[fliplr（Xa）,Xa（2:
501）];
subplot（2,2,1）;
plot（t*1000,xa）;%画出模拟信号，横坐标为时间（毫秒），纵坐标为幅度
xlabel（'time（millisecond）'）;ylabel（'xa（t）'）;
title（'anologsignal'）;
subplot（2,2,2）;
plot（W/（2*pi*1000）,Xa*1000）;%画出连续时间傅立叶变换
xlabel（'frequency（kHZ）'）;%横坐标为频率（kHz）
ylabel（'xa（jw）'）;%纵坐标为幅度
title（'FT'）;
%下面为采样频率5kHz时的程序
T=0.0002;%采样间隔为
n=-25:
1:
25;
x=exp（-1000*abs（n*T））;%离散时间信号
K=500;k=0:
1:
K;w=pi*k/K;%w为数字频率
X=x*exp（-j*n'*w）;%计算离散时间傅立叶变换（序列的傅立叶变换）
X=real（X）;
w=[-fliplr（w）,w（2:
K+1）];
X=[fliplr（X）,X（2:
K+1）];
subplot（2,2,3）;
stem（n*T*1000,x）;%画出采样信号（离散时间信号）
xlabel（'time（millisecond）'）;
ylabel（'x1（n）'）;
title（'discretesignal'）;
subplot（2,2,4）;
plot（w/pi,X）;%画出离散时间傅立叶变换
xlabel（'frequency（radian）'）;%横坐标为弧度
ylabel（'x1（jw）'）;title（'DTFT'）;
3、序列的傅里叶变换性质分析
%9已知序列
，
，求其傅里叶变换，并讨论其傅里叶变换的周期性和对称性。
n=0:
10
x=（0.9*exp（j*pi/3））.^n
k=-200:
200
[X,w]=ft1（x,n,k）
magX=abs（X）
angX=angle（X）
subplot（2,1,1）
plot（w/pi,magX）
grid
xlabel（'frequencyinpiunits'）
ylabel（'/X/'）
title（'MagnitudePart'）
subplot（2,1,2）
plot（w/pi,angX/pi）
grid
xlabel（'frequencyinpiunits'）
ylabel（'Radians/pi'）
title（'AnglePart'）
由图可见，序列
的傅里叶变换对
是周期的，但不是共轭对称的。
%10、已知序列
，
，求其傅里叶变换，并讨论其傅里叶变换的周期性和对称性。
n=-5:
5
x=（-0.9）.^n
k=-200:
200
[X,w]=ft1（x,n,k）
magX=abs（X）
angX=angle（X）
subplot（2,1,1）
plot（w/pi,magX）
grid
xlabel（'frequencyinpiunits'）
ylabel（'/X/'）
title（'MagnitudePart'）
subplot（2,1,2）
plot（w/pi,angX/pi）
grid
xlabel（'frequencyinpiunits'）
ylabel（'Radians/pi'）
title（'AnglePart'）
由图可见，序列
的傅里叶变换对
是周期的，是共轭对称的。
4、时域差分方程的求解
采用filter函数实现线性常系数差分方程的递推求解，函数调用格式如下：
●yn=filter（B,A,xn）计算输入信号xn的零状态响应yn
●yn=filter（B,A,xn,xi）计算输入信号xn的全响应yn，xi为等效初始条件的输入序列
●xi=filtic（B,A,ys,xs）由初始条件计算xi的函数
4.1求解差分方程y（n）＋a1y（n-1）＋a2y（n-2）=b0x（n）＋b1x（n-1）的零状态响应和全响应。
已知X（n）为单位取样序列，y（-1）=1,y（-2）=2,a1=0.5,a2=0.06,b0=2,b1＝3。
程序：
xn=[1zeros（1,20）]
B=[2,3]
A=[1,0.5,0.06]
ys=[1,2]
xi=filtic（B,A,ys）
yn1=filter（B,A,xn）
yn2=filter（B,A,xn,xi）
subplot（2,1,1）
n1=0:
length（yn1）-1
stem（n1,yn1,'.'）
axis（[0,21,-3,3]）
subplot（2,1,2）
n2=0:
length（yn2）-1
stem（n2,yn2,'.'）
结果图形：
上图为零状态响应、下图为全响应。
5、离散系统的Z域分析
%11利用系统函数
分析系统的稳定性。
假设系统函数如下式：
，试判断系统是否稳定。
解：
%调用roots函数求极点，并判断系统的稳定性
A=［3，－3.98，1.17，2.3418，－1.5147］；
%H（z）的分母多项式系数
p=roots（A）%求H（z）的极点
pm=abs（p）；%求H（z）的极点的模
ifmax（pm）<1disp（′系统因果稳定′），else，
disp（′系统不因果稳定′），end
程序运行结果如下：
极点：
－0.74860.6996-0.7129i 0.6996+0.7129i0.6760
pm=0.74860.99880.99880.6760
由极点分布可知系统因果稳定。
附录3：
例：
求解差分方程y（n）＋a1y（n-1）＋a2y（n-2）=b0x（n）＋b1x（n-1）的零状态响应或全响应。
已知X（n）为单位取样序列，y（-1）=1,y（-2）=2,a1=0.5,a2=0.06,b0=2,b1＝3。
程序：
B=[2,3];
A=[1,0.5,0.06];
ys=[1,2];
xn=[1zeros（1,20）];
aa=input（‘aa=’）;%1对应求零状态响应，其他值对应求全响应
ifaa==1
yn=filter（B,A,xn）;
else
xi=filtic（B,A,ys）;
yn=filter（B,A,xn,xi）;
end
n=0:
length（yn）-1;
stem（n,yn,'.'）
axis（[0,21,-3,3]）
程序流程图：
求零状态响应
求全响应
升级会员 