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

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

一、实验目的

1．了解时域离散信号的产生及基本运算实现。 2．掌握离散时间傅里叶变换实现及系统分析方法。 3. 熟悉离散时间傅里叶变换性质。 4. 掌握系统Z域分析方法。

5. 培养学生运用软件分析、处理数字信号的能力。 二、实验设备

1、计算机

2、Matlab7.0以上版本 三、实验内容

1、对于给定的时域离散信号会进行频谱分析，即序列的傅里叶变换及其性质分析。 2、对于离散系统会进行频域分析及Z域分析。包括频谱特性、零极点画图、稳定性分析。 3、对于差分方程会用程序求解，包括求单位冲击序列响应，零输入响应、零状态响应、全响应，求其系统函数，及其分析。

4、信号时域采样及其频谱分析，序列恢复。 5、扩展部分主要是关于语音信号的读取及其播放。 四、实验原理

1、序列的产生及运算

在Matlab中自带了cos、sin、exp（指数）等函数，利用这些函数可以产生实验所需序列。

序列的运算包括序列的加法、乘法，序列x(n)的移位x(nn0)，翻褶x(n)等。序列的加法或乘法指同序号的序列值逐项对应相加或相乘，但Matlab中“+”“.*”运算是对序列的值直接进行加或乘，不考虑两序列的序号是否相同，因此编程时考虑其序号的对应。

2、序列的傅里叶变换及其性质 序列的傅里叶变换定义：X(e)jnx(n)ejn|X(ej)|ej()，其幅度特性为|X(ej)|，

在Matlab中采用abs函数；相位特性为()，在Matlab中采用angle函数。

序列傅里叶变换的性质：

标准文案

实用文档

（1）FT的周期性X(ej(2M))X(ej)，实序列傅里叶变换的对称性X(ej)X(ej)。对实序列和复序列分别进行傅里叶变换，通过图形结果观察周期性即对称性。

（2）FT的频移特性FT[ej0nx(n)]X(ej(0))，对序列在时域乘以e换，比较其结果和直接对序列进行傅里叶变换的不同。

（3）时域卷积定理：若y(n)x(n)*h(n)，对序列x(n)和h(n)进行线性卷积得到y(n)，分别对它们进行傅里叶变换，应满足Y(ej)X(ej)H(ej)。

3、离散时间系统的Z域分析

已知离散时间系统的差分方程为aky(nk)bkx(nk)，对等号两边进行Z变换，得

k0k0NMj0n，然后进傅里叶变

到其系统函数H(z)及系统零极点，对系统函数进行反变换得到单位取样响应h(n)，根据单位取样响应或系统函数的系数可以得到频率响应H(ej)，根据极点位置判断系统稳定性。

4、信号时域采样及恢复

给定连续信号xa(t)，对其用不同的采样频率进行采样，根据时域采样定理，采样信号的频谱是原模拟信号频谱沿频率轴以s为周期延拓而成的，并且要不失真地还原出模拟信号时，要满足s2c，因此当采样频率满足和不满足采样定理时，所得到的频谱是不同的。

根据采样信号进行信号恢复时，采用内插公式xa(t)五、实验步骤

1、序列的基本运算

1.1 产生余弦信号x(n)cos(0.04n)及带噪信号y(n)cos(0.04n)0.2w(n) 0<=n<=50（噪声采用randn函数）

1.2 已知x1(n)2n1 1n5，x2(n)2n2 2n6,求两个序列的和、乘积、序列x1的移位序列（右移2位）,序列x2的翻褶序列，画出原序列及运算结果图。

2、序列的傅里叶变换

2.1 已知序列x(n)(0.5)nu(n)。试求它的傅里叶变换，并且画出其幅度、相角、实部和虚部的波形，并分析其含有的频率分量主要位于高频区还是低频区。

标准文案

nxa(nT)sin((tnT)/T)实现。

(tnT)/T实用文档

2.2 令xa(t)e1000|t|，求其傅立叶变换Xa(j)。分别用fs1kHz和fs5kHz对其进行采样，求出离散时间傅立叶变换X(ej)，画出相应频谱，分析结果的不同及原因。

3、序列的傅里叶变换性质分析

3.1 已知序列x(n)(0.9ej/3)n，0n10，求其傅里叶变换，并讨论其傅里叶变换的周期性和对称性。

3.2 已知序列x(n)(0.9)n，5n5，求其傅里叶变换，并讨论其傅里叶变换的周期性和对称性。

为了方便，考虑在两个周期，例如[2,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 利用系统函数H(z)分析系统的稳定性。假设系统函数如下式：

H(z)(z9)(z3)，试判断系统是否稳定。 4323z3.98z1.17z2.3418z1.51470.10.3z15.2 已知线性时不变系统的系统函数H(z)，编写程序求其单位取样响1210.8z0.12z应，频率响应及系统零极点，并画出相应图形。

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），给出程序内容的解释，并对结果进行分析。

七、思考题

下面四个二阶网络的系统函数具有一样的极点分布：

10.3z110.8z11）H1(z) 2）H2(z) 121211.6z0.9425z11.6z0.9425z10.8z111.6z10.8z23）H3(z) 4）H4(z)

11.6z10.9425z211.6z10.9425z2请分析研究零点分布对于单位脉冲响应的影响。 要求： （1） 分别画出各系统的零、 极点分布图；

（2） 分别求出各系统的单位脉冲响应，并画出其波形；

（3） 分析零点分布对于单位脉冲响应的影响。 八、实验参考资料

标准文案

实用文档

1、高西全,丁玉美.数字信号处理[M].西安：西安电子科技大学出版社,2008 2、张德丰.详解MATLAB 数字信号处理[M].北京：电子工业出版社,2010 3、王月明,张宝华.MATLAB基础与应用教程[M].北京:北京大学出版社,2012

附：实验报告要求：

实验名称：-------

班级： 组号： 姓名1（学号）、姓名2（学号）、姓名3（学号） 一、实验目的(手写)

二、实验主要内容(要根据自己组所做内容写，做了的写，没做的不要写) 例如：

1. 对序列的产生和运算方法进行实现 2. 序列的傅里叶变换实现、性质及分析 等等

三、实验主要仪器、设备及软件(手写) 四、实验步骤、结果与分析 例如： 1. 序列的运算

序列为……，进行加法、乘法、……运算 运算结果为……

2. 序列的傅里叶变换实现及分析

（1）已知序列x(n)(0.5)nu(n)。试求它的傅里叶变换，并且画出其幅度、相角、实部和虚部的波形，并分析其含有的频率分量主要位于高频区还是低频区。 程序 结果 分析

（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. 求序列x(n)(0.5)nu(n)的傅里叶变换

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('frequency in pi units') title('Magnitude Part') ylabel('Magnitude') subplot(2,2,3) plot(w/pi,angX) grid

xlabel('frequency in pi units') title('Angle Part') ylabel('Radians') subplot(2,2,2) plot(w/pi,realX) grid

xlabel('frequency in pi units') title('Real Part') ylabel('Real') subplot(2,2,4) plot(w/pi,imagX) grid

xlabel('frequency in pi units') title('Imaginary Part') ylabel('Imaginary')

程序执行结果：

标准文案

实用文档

Magnitude Part221.5Real PartMagnitude1.510.5Real10.500.5frequency in pi unitsAngle Part100.5frequency in pi unitsImaginary Part10-0.20-0.2-0.4-0.6-0.800.5frequency in pi units1-0.4-0.6-0.800.5frequency in pi units1ImaginaryRadians

%8 令xa(t)e1000|t|，绘制其傅立叶变换Xa(j)。用不同频率对其进行采样，分别画出

X(ej)。

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; %k*Wmax 求模拟角频率

KXa=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('anolog signal'); subplot(2,2,2);

plot(W/(2*pi*1000),Xa*1000); %画出连续时间傅立叶变换 xlabel('frequency(kHZ)'); %横坐标为频率（kHz） ylabel('xa(jw)'); %纵坐标为幅度 title('FT');

%下面为采样频率5kHz时的程序

T=0.0002; %采样间隔为1f0.0002s

sn=-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('discrete signal'); subplot(2,2,4);

plot(w/pi,X); %画出离散时间傅立叶变换 xlabel('frequency(radian)'); %横坐标为弧度 ylabel('x1(jw)');title('DTFT');

3、序列的傅里叶变换性质分析

%9 已知序列x(n)(0.9ej/3)n，0n10，求其傅里叶变换，并讨论其傅里叶变换的周期性和对称性。 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('frequency in pi units') ylabel('/X/')

title('Magnitude Part') subplot(2,1,2) plot(w/pi,angX/pi) grid

xlabel('frequency in pi units') ylabel('Radians/pi') title('Angle Part')

Magnitude Part86/X/420-2.5-2-1.5-1-0.500.5frequency in pi unitsAngle Part11.520.5Radians/pi0-0.5-2.5-2-1.5-1-0.500.5frequency in pi units11.52

标准文案

实用文档

由图可见，序列x(n)(0.9ej/3)n的傅里叶变换对w是周期的，但不是共轭对称的。 %10、已知序列x(n)(0.9)n，5n5，求其傅里叶变换，并讨论其傅里叶变换的周期性和对称性。

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('frequency in pi units') ylabel('/X/')

title('Magnitude Part') subplot(2,1,2) plot(w/pi,angX/pi) grid

xlabel('frequency in pi units') ylabel('Radians/pi') title('Angle Part')

Magnitude Part1510/X/50-2.5-2-1.5-1-0.500.5frequency in pi unitsAngle Part11.5210.50-0.5-1-2.5-2-1.5-1-0.500.5frequency in pi units11.52Radians/pi

由图可见，序列x(n)(0.9)n的傅里叶变换对w是周期的，是共轭对称的。

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=[1 zeros(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,'.')

结果图形：上图为零状态响应、下图为全响应。

标准文案

实用文档

3210-1-2-3024681012141618203210-1-202468101214161820

5、离散系统的Z域分析

%11 利用系统函数H(z)分析系统的稳定性。假设系统函数如下式：

H(z)(z9)(z3)，试判断系统是否稳定。

3z43.98z31.17z22.3418z1.5147 解：

%调用roots函数求极点， 并判断系统的稳定性

A=［3， －3.98， 1.17， 2.3418， －1.5147］； %H(z)的分母多项式系数

p=roots(A) %求H(z)的极点

pm=abs(p)； %求H(z)的极点的模

if max(pm)<1 disp(′系统因果稳定′)， else， disp(′系统不因果稳定′)， end 程序运行结果如下：

极点： －0.7486 0.6996-0.7129i 0.6996+0.7129i 0.6760

pm = 0.7486 0.9988 0.9988 0.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=[1 zeros(1,20)];

aa=input(‘aa=’);% 1对应求零状态响应，其他值对应求全响应 if aa==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])

程序流程图：

开始 设置差分方程参数 设置初始条件、输入 读取响应类型aa aa==1？ Y 求零状态响应 N 求全响应 画图 结束

标准文案