信号与线性系统综合实验

一、实验目的

  1、掌握连续时间信号与系统的时域、频域综合分析方法；
  2、掌握运用Matlab软件分析连续时间信号与系统的时域、频域特性；
  3、通过对连续时间信号与系统的综合分析，加深对信号频谱、系统函数、系统频率特性、冲激响应、阶跃响应等概念的理解，了解系统函数零、极点分布与系统的频率特性、稳定性之间的关系。

二、实验内容及其结果分析

（一）基础部分

  （1）设计一个2阶或3阶连续时间系统，分别用不同的实现结构（直接、级联、并联）实现所设计的系统，给出系统的S域模拟框图和信号流图，并运用Matlab软件对系统的时域、频域特性进行分析。
答：假设3阶连续时间系统为$$H(s)=\frac{s+1}{s^3+5s^2+6s}=\frac{s+1}{s(s+2)(s+3)}$$
S域模拟框图和信号流图如下：

运用Matlab软件对系统的时域特性进行分析：

>>syms s;
Hs=(s+1)/(s^3+5*s^2+6*s);
ht=ilaplace(Hs)
ezplot(ht,[0:0.01:20]);
title('时域波形
');xlabel('t');ylabel('h(t)');
ht =
exp(-2*t)/2 - (2*exp(-3*t))/3 + 1/6

运用Matlab软件对系统的频域特性进行分析。

>> T=0.02;
t=0:T:10;                                  
ht=exp(-2*t)/2 - (2*exp(-3*t))/3 + 1/6;
W1=pi/T;                                 
N=200;                                      
k=-N:N;
W=k*W1/(2*N+1);                             
Hjw=ht*exp(-j*t'*W);                          
subplot(2,1,1)
plot(W/2/pi,abs(Hjw));
xlabel('f(Hz)');ylabel('幅度');
title('幅度频谱');
subplot(2,1,2)
plot(W/2/pi,angle(Hjw)/pi*180)
xlabel('f(Hz)');ylabel('相角(度）');
title('相位频谱');

结果分析：设计的3阶连续时间系统$H(s)=(s+1)/(s^3+5s^2+6s)$的时域波形和幅度频谱如上图所示。其中从时域波形上来看，随着时间的增大，h(t)会逐渐的趋于1/6。而它的幅度频谱只在频率分量0附件有一个类似于冲激函数的幅度。

（2）改变上述系统的参数，观察并分析系统的各种响应特性，观察并分析系统的零、极点的变化情况与系统稳定性的关系。
方案一：上述的三阶系统$$H(s)=\frac{s+1}{s^3+5s^2+6s}=\frac{s+1}{s(s+2)(s+3)}$$

>>figure(1)
a=[1,5,6,0];
b=[1,1];
lxljdt(a,b);
figure(2)
[H,w]=freqs(b,a);
subplot(2,1,1)
plot(w,abs(H));
title('幅频响应');
subplot(2,1,2)
plot(w,angle(H))
title('相频响应');
figure(3)
step(b,a);
title('阶跃响应');
figure(4)
impulse(b,a);
title('冲激响应');

结果分析：从H(s)的零极点图中可以看出来，它的两个极点位于左半平面，一个极点正好位于虚轴上，这个系统是稳定的。

方案二：改变系统参数为$$H(s)=\frac{s+1}{s^3-s^2-6s}=\frac{s+1}{s(s+2)(s-3)}$$

>> figure(1)
a=[1,-1,-6,0];
b=[1,1];
lxljdt(a,b);
figure(2)
[H,w]=freqs(b,a);
subplot(2,1,1)
plot(w,abs(H));
title('幅频响应');
subplot(2,1,2)
plot(w,angle(H))
title('相频响应');
figure(3)
step(b,a);
title('阶跃响应');
figure(4)
impulse(b,a);
title('冲激响应') 

结果分析：从H(s)的零极点图中可以看出来，它的某个极点位于右半平面，因此这个系统是不稳定的。

方案三：改变系统参数为$$H(s)=\frac{s+1}{s^3+4s}=\frac{s+1}{s(s+2i)(s-2i)}$$

>> figure(1)
a=[1,0,4,0];
b=[1,1];
lxljdt(a,b);
figure(2)
[H,w]=freqs(b,a);
subplot(2,1,1)
plot(w,abs(H));
title('幅频响应');
subplot(2,1,2)
plot(w,angle(H))
title('相频响应');
figure(3)
step(b,a);
title('阶跃响应');
figure(4)
impulse(b,a);
title('冲激响应');

结果分析：从H(s)的零极点图中可以看出来，它的三个极点均位于虚轴上，这个系统是不稳定的。
总体结果分析：1.程序设计：我百度了一下，zplane(B,A)是求离散系统函数H(z)的零极点图，而我用的是连续系统函数H(s)，因此我在这里用的是lxljdt(a,b)，这是实验三指导书中的一个M函数，作用就是求连续系统函数H(s) 的零极点图。
2.幅频响应：从图中可以观察得出，当H(s)的极点为a+jb时，其中b表示幅频响应的峰值，如方案三所示，当b=0时，幅频响应只在0处有个峰值，如方案一、二所示。
3. 零、极点与系统稳定性：当零、极点全部位于左半平面或者仅且只有一个极点位于虚轴上，而右边平面没有极点，此时系统是稳定的

（3）给定一个输入激励信号，观察并分析该信号通过上述系统后输出响应的时域的变化情况，并分析其原因。

>> subplot(1,2,1)
syms t s;
f=heaviside(t); 
ezplot(f,[0 100]) 
subplot(1,2,2)
H=(s+1)/(s^3+5*s^2+6*s);
F=1/s;
f1=ilaplace(H*F)
ezplot(f1,[0 100])
f1 =
t/6 - exp(-2*t)/4 + (2*exp(-3*t))/9 + 1/36

结果分析：假定输入激励信号为单位阶跃信号，从时域上来看，阶跃信号通过该系统后变成了斜变信号。变化原因：可以从求解的表达式入手。首先我先求S域的H(s)和F(s)乘积，然后再求其拉普拉斯逆变换，求出的表达式为$\frac{t}{6}-\frac{1}{4}{e}^{-2t}+\frac{2}{9}{e}^{-3t}+\frac{1}{36}$，对应上图用matlab画出来的图。

（二）拓展部分

  （1）构建一个包含若干个不同频率分量的周期连续信号（各分量频率自定）f(t)，截取该信号的不同长度（注意截取长度应不小于最低频率分量的一个周期），分别用Matlab软件分析所截取信号的频谱（画出频谱图，含幅度频谱和相位频谱）。比较所截取的不同长度信号频谱的差异，同时与理论频谱进行比较，并运用所学知识，分析产生这些差异的原因。（复习频域卷积定理）
答：假定$f(t)=cos(t)+cos(10t)+cos(20t)$
方案一：截取长度为10，对应的门函数为：$g_{10}(t)$

>> T=0.02;                                    
t=-40:T:40;
f0=cos(t)+cos(10*t)+cos(20*t);
f1=heaviside(t+5)-heaviside(t-5);
f=f1.*f0;
W1=pi/T;                                   
N=1600;                                      
k=-N:N;
W=k*W1/(2*N+1);                             
F=f*exp(-j*t'*W);                          
subplot(2,1,1)
plot(W,abs(F));
xlabel('频率');ylabel('幅度');
title('幅度频谱');
subplot(2,1,2)
plot(W,angle(F)/pi*180)
xlabel('频率');ylabel('相角(度）');
title('相位频谱');

 
>> w=-40:0.01:40;
y=10*sin(5*(w+1))./(5*(w+1));
y0=10*sin(5*(w-1))./(5*(w-1));
y1=10*sin(5*(w+10))./(5*(w+10));
y2=10*sin(5*(w-10))./(5*(w-10));
y3=10*sin(5*(w+20))./(5*(w+20));
y4=10*sin(5*(w-20))./(5*(w-20));
plot(w,y);
hold on;plot(w,y0);
hold on;plot(w,y1);
hold on;plot(w,y2);
hold on;plot(w,y3);
hold on;plot(w,y4);
title('方案一的理论频谱')

方案二：截取长度为30，对应的门函数为：g_30 (t)

> T=0.02;                                    
t=-40:T:40;
f0=cos(t)+cos(10*t)+cos(20*t);
f1=heaviside(t+15)-heaviside(t-15);
f=f1.*f0;
W1=pi/T;                                   
N=1600;                                      
k=-N:N;
W=k*W1/(2*N+1);                             
F=f*exp(-j*t'*W);                          
subplot(2,1,1)
plot(W,abs(F));
xlabel('频率');ylabel('幅度');
title('幅度频谱');
subplot(2,1,2)
plot(W,angle(F)/pi*180)
xlabel('频率');ylabel('相角(度）');
title('相位频谱');
 
>> w=-40:0.01:40;
y=30*sin(15*(w+1))./(15*(w+1));
y0=30*sin(15*(w-1))./(15*(w-1));
y1=30*sin(15*(w+10))./(15*(w+10));
y2=30*sin(15*(w-10))./(15*(w-10));
y3=30*sin(15*(w+20))./(15*(w+20));
y4=30*sin(15*(w-20))./(15*(w-20));
plot(w,y);
hold on;plot(w,y0);
hold on;plot(w,y1);
hold on;plot(w,y2);
hold on;plot(w,y3);
hold on;plot(w,y4);
title('方案二的理论频谱')

方案三：截取长度为60，对应的门函数为：g_30 (t)

>> T=0.02;                                    
t=-40:T:40;
f0=cos(t)+cos(10*t)+cos(20*t);
f1=heaviside(t+30)-heaviside(t-30);
f=f1.*f0;
W1=pi/T;                                   
N=1600;                                      
k=-N:N;
W=k*W1/(2*N+1);                             
F=f*exp(-j*t'*W);                          
subplot(2,1,1)
plot(W,abs(F));
xlabel('频率');ylabel('幅度');
title('幅度频谱');
subplot(2,1,2)
plot(W,angle(F)/pi*180)
xlabel('频率');ylabel('相角(度）');
title('相位频谱');
 
>> w=-40:0.01:40;
y=60*sin(30*(w+1))./(30*(w+1));
y0=60*sin(30*(w-1))./(30*(w-1));
y1=60*sin(30*(w+10))./(30*(w+10));
y2=60*sin(30*(w-10))./(30*(w-10));
y3=60*sin(30*(w+20))./(30*(w+20));
y4=60*sin(30*(w-20))./(30*(w-20));
plot(w,y);
hold on;plot(w,y0);
hold on;plot(w,y1);
hold on;plot(w,y2);
hold on;plot(w,y3);
hold on;plot(w,y4);
title('方案三的理论频谱')

结果分析：从图中可以看出，当截取信号的长度越长时，幅度频谱在频率分量为±1，±10，±20上的外围包络线越接近为一个冲激响应，同时相同频率分量上对应的幅度变化也越来越大。和它的理论频谱的变化差异不是太大，但是他的波动幅度还是有一定的差异的。
原因：
  1.由$g\tau (t)\leftrightarrow \tau s_a(\omega \tau /2)$可知，截取信号的长度越长，即脉宽τ越大时，对应的取样函数的高度越大，频谱搬移后，相同频率分量上对应的幅度变化同样也是越来越大。另外，长度越长，对应的频带也就越窄，由这些取样函数构成的图形也就越窄，所以上图也就越接近冲激响应的形式。
  2.当f(t)被截取不同的长度，也就是f(t)乘以不同长度的门函数时，理论频谱对应是门函数的傅里叶变换取样函数在频域上的搬移，但是在用matlab操作时，理论频谱上重叠的部分会叠加起来，因此会和理论频谱有所差异。

  （2）设计一个LTI连续系统作为滤波器，对上一步构建的信号进行滤波（滤除部分频率分量，保留另一部分频率分量）。通过在S平面分布适当的零极点来设计该滤波器的系统函数H(s)，用Matlab软件分析所设计滤波器的频率响应（含幅频响应和相频响应），并画出频率响应波形。运用所学知识，分析幅频响应与零极点位置的关系。
答：$f(t)=cos(t)+cos(10t)+cos(20t)$，该信号中角频率为±10的分量滤除掉，另外的分量保留。此时H(s)的零点为±j10，极点为-0.2±j, -0.2±j20,故$$H(s)=\frac{S^2+100}{(S^2+0.4S+1.04)(S^2+0.4S+400.04)}=\frac{s^2+100}{s^4+0.8s^3+401.24s^2+160.432s+416.0416}$$

结果分析：H(s)虚轴上的极点a±jb中，b对应的幅频响应的峰值，虚轴上的零点a±jb中，b对应的幅频响应的谷值。在设计这个滤波器时，我考虑到了系统的稳定性，所以设计极点a±jb时，a不能直接为0，我取一个小的负数-0.2，使其极点位于左半平面，保证系统的稳定性。而零点a±jb不用考虑这些，a=0即可。从幅频响应的图中也可以看出保留了角频率为±1，±20的分量。

  （3）基于构建的f(t)和H(s)，采用数值计算的方法，求出滤波输出信号y(t)，并求出y(t)的频率响应Y(jω)，利用Matlab画出y(t)和Y(jω)的波形。

>> T=0.02;                                    
t=-40:T:40;
f=cos(t)+cos(20*t);
subplot(3,1,1)
plot(t,f),title('y(t)')
xlabel('t');ylabel('y(t)');
subplot(3,1,2)
W1=pi/T;                                   
N=1600;                                      
k=-N:N;
W=k*W1/(2*N+1);                             
F=f*exp(-j*t'*W);                          
plot(W,abs(F));
xlabel('频率');ylabel('幅度');
title('Y(jω)的幅度频谱');
subplot(3,1,3)
plot(W,angle(F)/pi*180)
xlabel('频率');ylabel('相角(度）');
title('Y(jω)的相位频谱');

结果分析：最开始我是利用s域求逆变换得到y(t)，但是不知道为什么画不出对应的波形。所以上面的程序是按照理论上的$f(t)=cos(t)+cos(20t)$来进行绘制的。

（三）应用设计部分

  利用Matlab进行音频信号的获取与输出，并进行相应的分析，绘制信号波形及频谱图等。

[pyr,fs]=audioread('pyr.wav');%声音读取
sound(pyr,fs);   %声音回放
n=length(pyr);
pyr1=fft(pyr,n);  %快速傅里叶变换
figure(1)
subplot(2,1,1);
plot(pyr);  
xlabel('时间');
ylabel('幅度');
title('初始信号波形');   %绘出时域波
subplot(2,1,2);  %绘出频域频谱
plot(abs(fftshift(pyr1)));
title('初始信号频谱');
xlabel('频率');
ylabel('幅度');

三、心得体会

  本次实验是综合性实验，我负责的是基础部分和扩展部分。通过对连续时间信号与系统的综合分析，我加深了对信号频谱、系统函数、系统频率特性、冲激响应、阶跃响应等概念的理解，进一步了解到了系统函数零、极点分布与系统的频率特性、稳定性之间的关系。通过本次实验，我拾起来了很多我遗忘的知识。比如：基础部分中H(s)的S域模拟框图和信号流图；再者，扩展部分中频谱的搬移，当f(t)乘以一个门函数时，相当于频谱中门函数对应的取样函数进行了频谱搬移，还有就是当门函数的脉宽变化时，频谱的特点。同时，我也学到了很多新的知识，例如：用zplane(B,A)求离散系统函数H(z)的零极点图；用freqs直接求H(s)对应的频率响应。我收益颇多。