实验二连续时间信的频域分析

实验二连续时间信号的频域分析一、实验目的1、掌握连续时间周期信号的傅里叶级数的物理意义和分析方法；2、观察截短傅里叶级数而产生的“Gibbs现象”，了解其特点以及产生的原因；3、掌握连续时间傅里叶变换的分析方法及其物理意义；4、掌握各种典型的连续时间非周期信号的频谱特征以及傅里叶变换的主要性质；5、学习掌握利用Matlab语言编写计算CTFS、CTFT和DTFT的仿真程序，并能利用这些程序对一些典型信号进行频谱分析，验证CTFT、DTFT的若干重要性质。

基本要求：掌握并深刻理傅里叶变换的物理意义，掌握信号的傅里叶变换的计算方法，掌握利用Matlab编程完成相关的傅里叶变换的计算。

二、原理说明1、连续时间周期信号的傅里叶级数CTFS分析任何一个周期为T1的正弦周期信号，只要满足狄利克利条件，就可以展开成傅里叶级数。

三角傅里叶级数为：∑∞=++=10 0)]sin()cos([)(kkktkbtkaatxωω 2.1或： ∑∞=++=100)cos()(k k k t k c a t x ϕω 2.2 其中102T πω=，称为信号的基本频率（Fundamental frequency ），k k b a a ，和，0分别是信号)(t x 的直流分量、余弦分量幅度和正弦分量幅度，k k c ϕ、为合并同频率项之后各正弦谐波分量的幅度和初相位，它们都是频率0ωk 的函数，绘制出它们与0ωk 之间的图像，称为信号的频谱图（简称“频谱”），k c －0ωk 图像为幅度谱，k ϕ－0ωk 图像为相位谱。

三角形式傅里叶级数表明，如果一个周期信号x(t)，满足狄里克利条件，就可以被看作是由很多不同频率的互为谐波关系（harmonically related ）的正弦信号所组成，其中每一个不同频率的正弦信号称为正弦谐波分量 (Sinusoid component)，其幅度（amplitude ）为k c 。

也可以反过来理解三角傅里叶级数：用无限多个正弦谐波分量可以合成一个任意的非正弦周期信号。

指数形式的傅里叶级数为：∑∞-∞==k t jk k e a t x 0)(ω 2.3其中，k a 为指数形式的傅里叶级数的系数，按如下公式计算： ⎰--=2/2/1110)(1T T t jk k dt e t x T a ω 2.4 指数形式的傅里叶级数告诉我们，如果一个周期信号x(t)，满足狄里克利条件，那么，它就可以被看作是由很多不同频率的互为谐波关系（harmonically related ）的周期复指数信号所组成，其中每一个不同频率的周期复指数信号称为基本频率分量，其复幅度（complex amplitude ）为k a 。

这里“复幅度（complex amplitude ）”指的是k a 通常是复数。

上面的傅里叶级数的合成式说明，我们可以用无穷多个不同频率的周期复指数信号来合成任意一个周期信号。

然而，用计算机（或任何其它设备）合成一个周期信号，显然不可能做到用无限多个谐波来合成，只能取这些有限个谐波分量来近似合成。

假设谐波项数为N ，则上面的和成式为：∑-==N N k t jk k e a t x 0)(ω 2.5显然，N 越大，所选项数越多，有限项级数合成的结果越逼近原信号x(t)。

本实验可以比较直观地了解傅里叶级数的物理意义，并观察到级数中各频率分量对波形的影响包括“Gibbs ”现象：即信号在不连续点附近存在一个幅度大约为9%的过冲，且所选谐波次数越多，过冲点越向不连续点靠近。

这一现象在观察周期矩形波信号和周期锯齿波信号时可以看得很清楚。

2、连续时间信号傅里叶变换----CTFT傅里叶变换在信号分析中具有非常重要的意义，它主要是用来进行信号的频谱分析的。

傅里叶变换和其逆变换定义如下：⎰∞∞--=dt e t x j X t j ωω)()( 2.6 ⎰∞∞-=ωωπωd e j X t x t j )(21)( 2.7连续时间傅里叶变换主要用来描述连续时间非周期信号的频谱。

按照教材中的说法，任意非周期信号，如果满足狄里克利条件，那么，它可以被看作是由无穷多个不同频率（这些频率都是非常的接近）的周期复指数信号e j?t 的线性组合构成的，每个频率所对应的周期复指数信号e j?t 称为频率分量（frequency component ），每个频率的幅度为对应频率的|X(j?)|之值，其相位为对应频率的X(j?)的相位。

X(j?)通常为复函数，可以按照复数的极坐标表示方法表示为：X(j?)=| X(j?)|e j? X(j?)其中，| X(j?)|称为x(t)的幅度谱，而?X(j?)则称为x(t)的相位谱。

给定一个连续时间非周期信号x(t)，它的频谱也是连续且非周期的。

对于连续时间周期信号，也可以用傅里叶变换来表示其频谱，其特点是，连续时间周期信号的傅里叶变换由冲激序列构成的，是离散的——这是连续时间周期信号的傅里叶变换的基本特征。

4.1 傅里叶级数的Matlab 计算设周期信号x(t)的基本周期为T 1，且满足狄里克利条件，则其傅里叶级数的系数可由式2.4计算得到。

式2.4重写如下：基本频率为： 102T πω=对周期信号进行分析时，我们往往只需对其在一个周期内进行分析即可，通常选择主周期（Principle period ）。

假定x 1(t)是x(t)中的主周期，则计算机不能计算无穷多个系数，所以我们假设需要计算的谐波次数为N ，则总的系数个数为2N+1个。

在确定了时间范围和时间变化的步长即T 1和dt 之后，对某一个系数，上述系数的积分公式可以近似为： ∑⎰---==n t jk n T T t jk k T dt e t x dt e t x T a 12/2/11/)()(10110ωω 对于全部需要的2N+1个系数，上面的计算可以按照矩阵运算实现。

Matlab 实现系数计算的程序如下：dt = 0.01;T = 2; t = -T/2:dt:T/2; w0 = 2*pi/T;x1 = input(‘Type in the periodic signal x(t) over one period x1(t)=’); N = input(‘Type in the number N=’);k = -N:N; L = 2*N+1;ak = x1*exp(-j*k*w0*t’)*dt/T;需要强调的是，时间变量的变化步长dt 的大小对傅里叶级数系数的计算精度的影响非常大，dt 越小，精度越高，但是，计算机计算所花的时间越长。

5、用Matlab 实现CTFT 及其逆变换的计算5.1 用Matlab 实现CTFT 的计算Matlab 进行傅里叶变换有两种方法，一种利用符号运算的方法计算，另一种是数值计算，本实验要求采用数值计算的方法来进行傅里叶变换的计算。

严格来说，用数值计算的方法计算连续时间信号的傅里叶变换需要有个限定条件，即信号是时限信号（Time limitedsignal），也就是当时间|t|大于某个给定时间时其值衰减为零或接近于零，这个条件与前面提到的为什么不能用无限多个谐波分量来合成周期信号的道理是一样的。

计算机只能处理有限大小和有限数量的数。

采用数值计算算法的理论依据是：若信号为时限信号，当时间间隔T取得足够小时，上式可演变为：上式用Matlab表示为：X=x*exp(j*t’*w)*T其中X为信号x(t)的傅里叶变换，w为频率Ω，T为时间步长。

相应的Matlab程序：T = 0.01; dw = 0.1; %时间和频率变化的步长t = -10:T:10;w = -4*pi:dw:4*pi;X(j?)可以按照下面的矩阵运算来进行：X=x*exp(-j*t’*?)*T; %傅里叶变换X1=abs(X); %计算幅度谱phai=angle(X); %计算相位谱为了使计算结果能够直观地表现出来，还需要用绘图函数将时间信号x(t)，信号的幅度谱|X(j?)|和相位谱? X(j?)分别以图形的方式表现出来，并对图形加以适当的标注。

四、实验步骤及内容6、编写程序，绘制下面的信号的波形图：其中，?0= 0.5π，要求将一个图形窗口分割成四个子图，分别绘制cos(?t)、cos(3?t)、cos(5?t) 和x(t) 的波形图，给图形加title，网格线和x坐标标签。

dt=0.002; %定义时间间隔t=-2:dt:2; %设置时间起始和截止w0=0.5*pi; %题目要求，定义wox1=cos(w0*t); %输入信号x1subplot(221); %将信号1的图形放在两行两列中的第一个plot(t,x1),grid on; %输出信号1的图形，并加网格x2=cos(3*w0*t); %输入信号x2subplot(222); %将信号2的图形放在两行两列中的第二个plot(t,x2),grid on; %输出信号2的图形，并加网格x3=cos(5*w0*t); %输入信号3subplot(223); %将信号3的图形放在两行两列中的第三个plot(t,x3),grid on; %输出信号3的图形，并加网格ak=0; %设傅里叶系数的初始值为0for n=1:9 %循环计算9次傅里叶系数ak=ak+(-1).^(n+1).*(1/(2*n-1)).*cos((2*n-1).*w0*t);%计算傅里叶系数endsubplot(224); %将该图形放在两行两列中的第四个plot(t,ak),grid on; %输入该图形，以t为横坐标，ak为纵坐标，并加网格7、给定如下两个周期信号：编写程序，计算x1(t)和x2(t)共21个系数；仿照程序Program2_212合成的y1(t)和y2(t)的波形图，以及x1(t) 和x2(t)的幅度频谱和相位频谱的谱线图。

反复执行该程序，输入不同的N值，观察合成的信号波形中，是否会产生Gibbs现象？为什么？clear,close allT=2;dt=0.00001; %定义周期T=2，时间间隔dt为0.00001，t=-T/2:dt:T/2; %时间轴的起始和截止x1=u(t+0.2)-u(t-0.2); %输入信号1x=0; %给x附初值为0t1=-2*T:dt:2*T; %定义t1时间轴的起始和截止for m=-2:2x=x+u(t1+0.2-m*T)-u(t1-0.2-m*T); %定期延长x1(t)形成一个周期性的信号endsubplot(221) %将该图形放在两行两列的第一个plot(t1,x) %以t1为横坐标，x为纵坐标画图axis([-2*T,2*T,-0.2,1.2]) %图形的横坐标从-2T到2T，纵坐标从-0.2到1.2title('x(t)') %命名x（t）xlabel('Time t') %轴命名Time tw0=2*pi/T; %定义频率w0N=10; %谐波组件的数量为10L=2*N+1;for k=-N:N;ak(N+1+k)=x1*exp(-j*k*w0*t')*dt/T; %计算傅里叶系数的值endphi=angle(ak); %计算ak的阶段值y=0;for q=1:L;y=y+ak(q)*exp(j*(q-1-N)*2*pi*t1/T); %分析周期性信号y从有限傅里叶级数(t) end;subp lot(222) %将图形放在两行两列中的第二位plot(t1,y) %以t1为横坐标，y为纵坐标画图axis([-2*T,2*T,-0.2,1.2]) %图形的横坐标从-2T到2T，纵坐标从-0.2到1.2title('有限项级数合成后的信号') %命名‘有限项级数合成后的信号’xlabel('Time t') %横轴命名Time tk=-N:N;subplot(223)stem(k,abs(ak))axis([-N,N,0,1]) %限定横纵坐标范围title('abs(ak)') %命名xlabel('Frequency index k') %横轴命名subplot(224)stem(k,phi) %作图axis([-N,N,-pi,pi]) %图形的横坐标从-N到N，纵坐标从-PI到PI title('phi') %命名phixlabel('Frequency index k') %横轴命名'Frequency index k'T = 2; dt = 0.00001; t = -2:dt:2; %定义周期，时间间隔，时间轴长度x1=(1-t).*(u(t)-u(t-1))+(t+1).*(u(t+1)-u(t)); x = 0; %输入信号x1for m = -1:1x = x +(1-t+2*m).*(u(t-2*m)-u(t-1-2*m))+(t+1-2*m).*(u(t+1-2*m)-u(t-2*m));%使信号周期化endw0 = 2*pi/T; %定义w0N = input('Type in the number of the harmonic components N = :'); %输入N值L = 2*N+1; %定义Lfor k = -N:1:N; %定义区间ak(N+1+k) = (1/T)*x1*exp(-j*k*w0*t')*dt; %做傅里叶变换endphi = angle(ak); %计算相位谱y=0;for q = 1:L;y = y+ak(q)*exp(j*(-(L-1)/2+q-1)*2*pi*t/T); %周期化end;subp lot(221) %将图形放在两行两列的第一位plot(t,x), title('The original signal x(t)'), axis([-2,2,-0.2,1.2]), %以t为横坐标x为纵坐标画图，并命名和确定横纵轴的长度subplot(223),plot(t,y), title('The synthesis signal y(t)'), axis([-2,2,-0.2,1.2]),xlabel('Time t'),subplot(222)k=-N:N; stem(k,abs(ak),'k.'), title('The amplitude |ak| of x(t)'), axis([-N,N,-0.1,0.6])subplot(224)stem(k,phi,'r.'), title('The phase phi(k) of x(t)'), axis([-N,N,-2,2]), xlabel('Index k')clear,close allT=2;dt=0.001;t=-T/2:dt:T/2; %定义周期，时间，时间间隔w=0.5*pi; %定义频率wN=10; %输入谐波组件的数量L = 2*N+1; %x1=cos((pi/2)*t).*(u(t+1)-u(t-1)); %输入信号x1subplot(241) %将图形放在两行两列的第一位plot(t,x1) %以t为横坐标，x1为纵坐标画图for k=-N:1:N %ak(N+1+k)=(1/T)*x1*exp(-j*k*w*t')*dt; %计算傅里叶级数endphi = angle(ak); %计算ak的阶段值y=0; %for q = 1:L; %y = y+ak(q)*exp(j*(-(L-1)/2+q-1)*2*pi*t/T); %分析周期性信号y从有限傅里叶级数(t)endsubplot(242) %将图形放在两行两列的第二位plot(t,y) %以t为横坐标，y为纵坐标画图subplot(243), %将图形放在两行两列的第三位k=-N:N; stem(k,abs(ak),'k.') %作图subplot(244) %将图形放在两行两列的第四位stem(k,phi,'r.') %作图t=-2:0.01:2; %定义时间轴长度w=-4*pi:0.1:4*pi; %定义频率w的长度及间隔X=x1*exp(-j*t'*w)*T; %x1的傅里叶变换x1=(t+2).*(u(t+2)-u(t+1))+1.*(u(t+1)-u(t-1))+(-t+2).*(u(t-1)-u(t-2)) %输入信号x1X2=abs(X); %计算幅度谱phai=angle(X); %计算相位谱x3=real(X); %X的实部信号赋给信号3x4=imag(X); % X的虚部信号赋给信号4subplot(231) %plot(w,X2); %title('幅度频谱') %命名为'幅度频谱'subplot(232) %plot(w,phai) %title('相位谱') %命名为'相位谱'subplot(233) %plot(t,x1) %title('x1波形') %命名为(' x1波形')axis([-2,2,-2,2]); %限定横纵坐标subplot(234) %plot(w,X) %title('傅里叶变换') %命名为‘傅里叶变换’subplot(235) %plot(w,x3) %title('实部') %命名为‘实部’subplot(236) %plot(w,x4) %title('虚部') %命名为‘虚部’9、从信号分解的角度，谈谈你对傅里叶级数、傅里叶变换及其物理意义和信号频谱概念的理解。