信号、系统及系统响应,离散系统的时域分析实验报告

实验报告

实验二 信号、系统及系统响应，离散系统的时域分析

一、实验目的

（1） 熟悉连续信号经理想采样前后的频谱变换关系，加深对时域采样定理的理

解；

（2） 熟悉时域离散系统的时域特性； （3） 利用卷积方法观察分析系统的时域特性；

（4） 掌握序列傅里叶变换的计算机实现方法，利用序列的傅里叶变换对连续信

号、离散信号及系统响应进行频域分析。

（5） 熟悉并掌握离散系统的差分方程表示法； （6） 加深对冲激响应和卷积分析方法的理解。

二、实验原理与方法

1、信号、系统及系统响应

采样是连续信号数字处理的第一个关键环节。对采样过程的研究不仅可以了解采样前后信号时域和频域特性发生的变化以及信号信息不丢失的条件，而且可以加深对傅里叶变换、Z 变换和序列傅里叶变换之间关系式的理解。

我们知道，对一个连续信号xa(t)进行理想采样的过程可用(2-1)表示。

^

()()()

(21)

a a x t x t p t =-

其中^

()a x t 为()a x t 的理想采样，()p t 为周期冲激脉冲，即

()()

(22)

n p t t nT δ∞

=-∞=

--∑

^

()a x t 的傅里叶变换^

()a X j Ω为

^

1()[()]

(23)

a a s m X j X j m T ∞

=-∞

Ω=Ω-Ω-∑

（2-3）式表明^

()a X j Ω为()a X j Ω的周期延拓，其延拓周期为采样角频率

(2/)s T πΩ=。其采样前后信号的频谱只有满足采样定理时，才不会发生频率混叠失真。

将（2-2）带入（2-1）式并进行傅里叶变换：

^

()[()()]j t a a n X j x t t nT e dt

δ∞

∞

-Ω-∞

=-∞

Ω=-∑?

[()()]j t a n x t t nT e dt

δ∞

∞

-Ω-∞

=-∞

=

-∑?

()(24)

j nT

a

n x nT e

∞

-Ω=-∞

=

-∑

式中()a x nT 就是采样后得到的序列()x n ，即

()()a x n x nT =

()x n 的傅里叶变换()j X e ω为

()()(25)

j j n

n X e x n e

ω

ω∞

-=-∞

=

-∑

比较(2-5)和(2-4)可知

在数字计算机上观察分析各种序列的频域特性， 通常对X(ej ω)在［0, 2π］上进行M 点采样来观察分析。 对长度为N 的有限长序列x(n), 有

一个时域离散线性非移变系统的输入/输出关系为

上述卷积运算也可以在频域实现

2、离散系统时域分析

^

()()

(26)

j a T

X j X e ωω=ΩΩ=-1

()()(27)

2,0,1,,1k N j n

j k

n k X e

x m e k k M M

ωωπ

ω--==-=

=???-∑()()()()()

(28)

m y n x n h n x m h n m ∞

=-∞

=*=

--∑()()()

(29)

j j j Y e X e H e ωωω=-式中

离散系统

]

[n x ]

[n y Discrete-time

systme

其输入、输出关系可用以下差分方程描述：

∑∑==-=-M

k k N

k k

k n x p k n y d

][][

输入信号分解为冲激信号，

∑∞

-∞

=-=

m m n m x n x ][][][δ。

记系统单位冲激响应

][][n h n →δ，

则系统响应为如下的卷积计算式：

∑∞

-∞

=-=

*=m m n h m x n h n x n y ][][][][][

当N k d k ,...2,1,0==时，h[n]是有限长度的（n ：[0，M]），称系统为FIR 系统；反之，称系统为IIR 系统。

三、实验内容及步骤

（1） 认真复习采样理论、 离散信号与系统、 线性卷积、 序列的傅里叶变换

及性质等有关内容， 阅读本实验原理与方法。

（2） 编制实验用主程序及相应子程序。

① 信号产生子程序， 用于产生实验中要用到的下列信号序列： a.采样信号序列：对下面连续信号：

x a (t)=Ae -at sin(Ω0t)u(t)

进行采样， 可得到采样序列

xa(n)=xa(nT)=Ae-anT sin(Ω0nT)u(n), 0≤n<50 其中A 为幅度因子， a 为衰减因子， Ω0是模拟角频率， T 为采样间隔。 这些参数都要在实验过程中由键盘输入， 产生不同的xa(t)和xa(n)。 b. 单位脉冲序列： xb(n)=δ(n)

c. 矩形序列： xc(n)=RN(n), N=10

②系统单位脉冲响应序列产生子程序。本实验要用到两种FIR系统。

a. h

a (n)=R

10

(n);

b. h

b

(n)=δ(n)+2.5δ(n-1)+2.5δ(n-2)+δ(n-3)

③有限长序列线性卷积子程序，用于完成两个给定长度的序列的卷积。可以直接调用MATLAB语言中的卷积函数conv。 conv用于两个有限长度序列的卷积，它假定两个序列都从n=0 开始。调用格式如下：

y=conv (x, h)

程序如下：

figure

subplot(121);Tp = 50/1000;Fs = 1000;%频率T = 1 / Fs;%周期

M = Tp * Fs;%采样点n = 0:Tp * Fs;dt = 0.00005;

t = -0.005:dt:0.005;A = 444.128;a = 50*2^0.5*pi;omiga = 50*2^0.5*pi;

xat = A*exp(-a*t).*sin(omiga*t);%连续信号

xnt = A*exp(-a*n*T).*sin(omiga*n*T);%50个采样点

stem(n,xnt,'.');box on;title('1kHz时域抽样曲线');%Fs=1000Hz的时域抽样曲线

Xk=T * fft(xnt,M);k = 0 : M-1;fk = k / Tp;

subplot(122);plot(fk,abs(Xk));title('1kHz时幅频特性曲线');%幅频特性曲线

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

%取采样频率fs=300 Hz

Figure;subplot(121);Fs=300;T=1/Fs;M=ceil(Tp*Fs);n=0:M-1;%取整

xnt=A*exp(-a*n*T).*sin(omiga*n*T); %xnt=xa（nt）

stem(n,xnt,'k.');box on;title('300Hz时域抽样曲线'); %Fs=300Hz的时域抽样

Xk=T*fft(xnt,M); %M点FFT[xnt)] k=0 : M-1;fk = k / Tp;

subplot(122);plot(fk,abs(Xk));title('300Hz时幅频特性曲线'); %xa（t）的幅频特性曲线

xlabel('f(Hz)');ylabel('幅度');axis([0,Fs,0,1.2*max(abs(Xk))]); %设置坐标比例

%取采样频率fs=200 Hz

Figure;subplot(121);Fs=200;T=1/Fs;M=ceil(Tp*Fs);n=0:M-1;%取整

xnt=A*exp(-a*n*T).*sin(omiga*n*T); %xnt=xa（nt）

stem(n,xnt,'k.');box on;title('200Hz时域抽样曲线'); %Fs=200Hz的时域抽样

Xk=T*fft(xnt,M); %M点FFT[xnt)] k=0 : M-1;fk = k / Tp;

subplot(122);plot(fk,abs(Xk));title('200Hz时幅频特性曲线'); %xa（t）的幅频特性曲线xlabel('f(Hz)');ylabel('幅度');axis([0,Fs,0,1.2*max(abs(Xk))]); %设置坐标比例

%绘制hb(n)函数

%hb(n)的时域序列

hb=[1,2.5,2.5,1];i=0:3;figure;subplot(121);stem(i,hb,'k.');axis([0 3 0 2.5]);

xlabel('n');ylabel('hb(n)');title('hb(n)的时域序列');

%hb(n)的傅氏变换|Hb(jw)|采样点为4

N=4;k=-200:200;w=k*pi/100;Hb=DFT(hb,N);subplot(122);plot(w/pi,abs(Hb));

xlabel('w/pi');ylabel('|Hb(jw)|');

title('hb(n)的傅氏变换|Hb(jw)|');

%通过单位冲击信号xb(n)的时域序列

xb=[1,0,0,0,0,0,0,0,0,0];i=0:9;figure;subplot(121);stem(i,xb,'k.');

xlabel('n');ylabel('xb(n)');

title('xb(n)的时域序列');

%xb(n)的傅氏变换(Xb|jw|)采样点为10

N=10;k=-200:200;w=k*pi/100;Xb=DFT(xb,N);magXb=abs(Xb);subplot(122);plot(w/pi,magXb); xlabel('w/pi');ylabel('|Xb(jw)|');

title('xb(n)的傅氏变换(Xb|jw|)');

%卷积yb(n)=xb(n)*hb(n)的时域序列

yb=conv(xb,hb);figure;subplot(121);stem(0:12,yb,'k.');

xlabel('n');ylabel('yb(n)=xb(n)*hb(n)');

title('yb(n)=xb(n)*hb(n)的时域序列');

%卷积yb(n)的傅氏变换(Yb|jw|)

N=13;k=-200:200;w=k*pi/100;Yb=DFT(yb,N);subplot(122);plot(w/pi,abs(Yb));

xlabel('w/pi');ylabel('|Yb(jw)|');

title('yb(n)的傅氏变换|Yb(jw)|');

%绘制ha(n)函数

%ya(n)=xc(n)*ha(n)的时域序列

ha=[1,1,1,1,1,1,1,1,1,1];xc=ha;ya=conv(ha,xc);figure;subplot(121);stem(0:18,ya,'k.'); xlabel('n');ylabel('ya(n)=xc(n)*ha(n)');

title('ya(n)=xc(n)*ha(n)的时域序列');

%ya(n)的傅氏变换(Ya|jw|)

N=19;k=-200:200;w=k*pi/100;Ya=DFT(ya,N);subplot(122);plot(w/pi,abs(Ya));

xlabel('w/pi');ylabel('|Ya(jw)|');

title('ya(n)的傅氏变换|Ya(jw)|');

%ya(n)=xc(n)*ha(n)的时域序列

ha=[1,1,1,1,1,1,1,1,1,1];xc=[1,1,1,1,1];ya=conv(ha,xc);figure;subplot(121);

stem(0:13,ya,'k.');xlabel('n');ylabel('ya(n)=xc(n)*ha(n)');

title('ya(n)=xc(n)*ha(n)的时域序列');

%ya(n)的傅氏变换(Ya|jw|)，采样点变

N=14;k=-200:200;w=k*pi/100;Ya=DFT(ya,N);subplot(122);plot(w/pi,abs(Ya));

xlabel('w/pi');ylabel('|Ya(jw)|');

（3） 调通并运行实验程序， 完成下述实验内容：

① 分析采样序列的特性。产生采样信号序列x a (n)，使

0444.128,

502,502

A a ππ==Ω=。右图给出了连续信号x a (t)的幅频特性曲线，由此图可以确定对x a (t)应采用的采样频率。

a. 取采样频率fs=1 kHz, 即T=1 ms 。 观察所得采样x a (n)的幅频特性|X(e j ω

)|和右图中的|X(e j Ω

)|，在折叠 频率附近有无明显差别。

图：1kHZ 时域抽样曲线与时幅频特性曲线

由上图可得知，当采样频率为1kHz 时，采样序列在折叠频率附近处，无明显频谱混叠。

b. 改变采样频率， fs=300 Hz ， 观察|X(e j ω)|的变化， 并做记录(打印曲线)； 进一步降低采样频率， fs=200 Hz ， 观察频谱混叠是否明显存在， 说明原因， 并记录(打印)这时的|X(e j ω)|曲线。

10.8

0.60.40.200

100

200300

400

500

x a (j f )

f /Hz

图：300HZ 时域抽样曲线与时幅频特性曲线

图：200HZ 时域抽样曲线与时幅频特性曲线

当fz 为200Hz 时，由于fz 太小，存在严重的频率混叠失真。

由以上两张图可得知，随着采样频率的降低，主瓣宽度随之变宽，频率混叠现象也逐渐严重。同时，伴随较为明显的失真现象。

② 时域离散信号、 系统和系统响应分析。

a. 观察信号x b (n)和系统h b (n)的时域和频域特性； 利用线性卷积求信号x b (n)通过系统h b (n)的响应y(n)， 比较所求响应y(n)和h b (n)的时域及频域特性， 注意它们之间有无差别，绘图说明并用所学理论解释所得结果。

图：h b (n)时间序列及其幅频谱

图：x b (n)时间序列及其幅频谱

图：y(n)=x b (n)*h b (n)时间序列及其幅频谱

b. 观察系统h a (n)对信号x c (n)的响应特性。 利用线性卷积求系统响应y(n)，并判断比较y(n)图形及其非零值序列长度是否与理论结果一致。调用序列傅里叶 变换数值计算子程序，求得()k j Y e ω,观察|()|k j Y e ω特性曲线，定性判断结果的正确性。改变x c (n)的长度，取N=5，重复该实验。注意参数变化的影响，说明变化前后的差异，并解释所得结果。

图：y(n)=x c (n)*h a (n)及()k j Y e ω，x c (n)=h a (n)=R 10(n)

图：y(n)=x c (n)*h a (n)及()k j Y e ω，x c (n)=R 5(n)

当N=10时，峰值较高，且峰值很窄，变换之后图形频带主值部分比较集中；N=5时，峰值较低，且峰值较宽，变换之后图形频带主值部分相对分散。

（4） 编制程序求解下列两个系统的单位冲激响应和阶跃响应，并绘出其图形。要求分别用 filter 、conv 、impz 三种函数完成。

]1[][]2[125.0]1[75.0][--=-+-+n x n x n y n y n y ]}4[]3[]2[]1[{25.0][-+-+-+-=n x n x n x n x n y

给出理论计算结果和程序计算结果并讨论。

式子一：]1[][]2[125.0]1[75.0][--=-+-+n x n x n y n y n y 。 在filter 函数下的单位冲激响应和阶跃响应的程序与图形如下： 程序：

%使用filter 函数,第一式子 ya = [1,0.75,0.125];xa = [1,-1];

n = 0 : 20; x1 = [1 zeros(1,20)];h = filter(xa,ya,x1);%调用单位冲激响应函数 figure;subplot(121);stem(n,h,'k.');xlabel('第一式子：n');ylabel('h(n)');

title('脉冲响应-filter函数');x2 = ones(1,21);

g = filter(xa,ya,x2); %调用单位阶跃响应函数

subplot(122);stem(n,g,'k.');xlabel('第一式子：n');ylabel('g(n)');title('阶跃响应-filter函数');

结果：

图：filter函数下的单位冲激响应和阶跃响应

在conv函数下的单位冲激响应和阶跃响应的程序与图形如下：

程序：

%使用conv函数,第一式子

ya = [1,0.75,0.125];xa = [1,-1];x2 = ones(1,21);h = impz(xa,ya,20);

y1 = conv(h,x2);

y1conv = y1(1:21); %为何y1conv要取y1中1：21的值，解释见

n1 = 0:20; %y2?单位阶跃响应?用conv函数中注释

figure;subplot(121);stem(n1,y1conv,'k.');

xlabel('第一式子：n');ylabel('g(n)');title('阶跃响应-conv函数');

x1 = [1 zeros(1,10)];h = impz(xa,ya,10);y1conv = conv(h,x1);n = 0 : 19;

subplot(122);stem(n,y1conv,'k.');xlabel

('第一式子：n');ylabel('h(n)');title('脉冲响应-conv函数');

结果：

图：conv函数下的单位冲激响应和阶跃响应

在impz 函数下的单位冲激响应和阶跃响应的程序与图形如下： 程序：

%使用impz 函数,第一式子

Figure ；subplot(121);ya = [1,0.75,0.125];xa = [1,-1,0];N = 20;n = 0 : N-1; h = impz(xa,ya,n);%调用单位冲激响应函数 g = dstep(xa,ya,n); %调用单位阶跃响应函数

stem(n,h,'k.'); ylabel('h(n)');xlabel('第一式子：n'); title('脉冲响应-impz 函数'); subplot(122);stem(n,g,'k.'); ylabel('g(n)');xlabel('第一式子：n'); title(' 阶跃响应

-impz 函数');

结果：

图：impz 函数下的单位冲激响应和阶跃响应

式子一的理论计算结果：

单位冲激响应:]2[125.0]1[75.0]1[][][------=n h n h n n n H δδ n 0 1 2 3 4 5 ... h[n] 1 -1.75 1.19 -0.67 0.355 -0.18 ... Y[n]

1

-1.75

1.19

-0.67

0.355

-0.18

...

单位阶跃响应:][]2[125.0]1[75.0][n n y n y n y δ+----= n 0 1 2 3 4 5 ... Y[n]

1

-0.75

0.44

-0.234

0.12

-0.06

...

结论：

由此可见，式子一的程序结果和理论计算结果一致。

式子二：]}4[]3[]2[]1[{25.0][-+-+-+-=n x n x n x n x n y

在filter 函数下的单位冲激响应和阶跃响应的程序与图形如下： 程序：

%使用filter 函数,第二式子

ya = [1,0,0,0];xa = [0.25,0.25,0.25,0.25];n = 0 : 20;x1 = [1 zeros(1,20)]; h = filter(xa,ya,x1);%调用单位冲激响应函数

figure;subplot(121);stem(n,h,'k.');xlabel('第二式子：n');ylabel('h(n)');title('脉冲响应-filter函数');n = 0 : 20;x2 = ones(1,21);g = filter(xa,ya,x2); %调用单位阶跃响应函数subplot(122);stem(n,g,'k.');xlabel('第二式子：n');ylabel('g(n)');title('阶跃响应-filter 函数');fprintf('实验结果: %f\n',sum(abs(h)));

结果：

图：filter函数下的单位冲激响应和阶跃响应

在conv函数下的单位冲激响应和阶跃响应的程序与图形如下：

程序：

%使用conv函数,第二式子

ya = [1,0,0,0];xa = [0.25,0.25,0.25,0.25];

x1 = [1 zeros(1,10)];n = 0 : 19;h = impz(xa,ya,10);y2conv = conv(h,x1);

figure;subplot(121);stem(n,y2conv,'k.');xlabel('第二式子：n');ylabel('h(n)');title('脉冲响应-conv函数');

h = [0 0.25*ones(1,4)];n = 0:20;x2 = ones(1,21);y2 = conv(h,x2);y2conv = y2(1:21);

subplot(122);stem(n,y2conv,'k.');xlabel('第二式子：n');ylabel('g(n)');title('阶跃响应-conv函数');

结果：

图：conv函数下的单位冲激响应和阶跃响应

在impz函数下的单位冲激响应和阶跃响应的程序与图形如下：

程序：

%使用impz函数，第二式子

ya = [1,0,0,0];xa = [0.25,0.25,0.25,0.25];N = 20;n = 0 : N-1;

h = impz(xa,ya,n);%调用单位冲激响应函数g = dstep(xa,ya,n); %调用单位阶跃响应函数

figure;subplot(121);stem(n,h,'k.'); ylabel('h(n)');xlabel('第二式子：n'); title('脉冲响应-impz函数');

subplot(122);stem(n,g,'k.'); ylabel('g(n)');xlabel('第二式子：n'); title('阶跃响应-impz函数');

结果：

图：impz函数下的单位冲激响应和阶跃响应

式子二的理论计算结果：

单位冲激响应:

n 0 1 2 3 4 5 ...

Y[n] 0 0.25 0.25 0.25 0.25 0 ...

单位阶跃响应:

n 0 1 2 3 4 5 ...

Y[n] 0 0.25 0.5 0.75 1 1 ...

结论：

由上可知，式子二的程序结果和理论计算结果一致。

2、思考题

在分析理想采样序列特性的实验中，采样频率不同时，相应理想采样序列的傅里叶变换频谱的数字频率度量是否都相同? 它们所对应的模拟频率是否相同? 为什么?

ωΩ?可知，当采样的频率不相同时，则周期T不同，于此同时，相应由=T

的数字频率ω也不相同；由于是同一信号，所以模拟频率Ω将保持不变。