实验一时域离散信号与系统变换域分析报告205.docx

实验一时域离散信号与系统变换域分析报告205.docx

《实验一时域离散信号与系统变换域分析报告205.docx》由会员分享，可在线阅读，更多相关《实验一时域离散信号与系统变换域分析报告205.docx（17页珍藏版）》请在冰豆网上搜索。

实验一时域离散信号与系统变换域分析报告205

实验一时域离散信号与系统变换域分析

一、实验目的

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]）

程序流程图：

求零状态响应

求全响应