试验二：连续时间信号的频域分析

1、、实验目的1、掌握连续时间周期信号的傅里叶级数的物理意义和分析方法；2、观察截短傅里叶级数而产生的“Gibb觇象”，了解其特点以及产生的原因；3、掌握连续时间傅里叶变换的分析方法及其物理意义；4、学习利用 MATLAB语言编写计算 CTFS和CTFT的仿真程序。基本要求：掌握并深刻理傅里叶变换的物理意义，掌握信号的傅里叶变换的计算方法，掌握利用MATLAB编程完成相关的傅里叶变换的计算。二、实验原理及方法1、连续时间周期信号的傅里叶级数CTFS分析任何一个周期为Ti的正弦周期信号，只要满足 狄利克利条件，就可以展开成傅里叶级数。 其中三角傅里叶级数为：x(t) a0ak cos(k d) bk

2、 sin(k d)9.1k 1或： x(t) a0 Ck cos(k 0tk)9.2k 1，一 2其中0，称为信号的 基本频率，a0, ak,和bk分别是信号x(t)的直流分量、余弦Ti分量幅度和正弦分量幅度，Ck、k为合并同频率项之后各正弦谐波分量的幅度和初相位，它们都是频率k 0的函数，绘制出它们与k 0之间的图像，称为信号的 频谱图(简称 频谱”),Ck k 0图像为幅度谱，k - k 0图像为相位谱。指数形式的傅里叶级数为：x(t)akejk 0t9.3k其中，ak为指数形式的傅里叶级数的系数，按如下公式计算：akTi/2x(t)eTi/2jk 0t dt9.4假设谐波项数为 N,则上

3、面的和成式为: TOC o 1-5 h z N k ot x(t) ake9.5k N显然，N越大，所选项数越多，有限项级数合成的结果越逼近原信号x(t)。2、连续时间信号傅里叶变换-CTFT傅里叶变换在信号分析中具有非常重要的意义，它主要是用来进行信号的频谱分析的。傅里叶变换和其逆变换定义如下：X(j ) x(t)e j tdt9.61-9.7x(t) X(j )ejtd连续时间傅里叶变换主要用来描述连续时间非周期信号的频谱。任意非周期信号，如果满足狄里克利条件，那么，它可以被看作是由无穷多个不同频率(这些频率都是非常的接近)的周期复指数信号 ej t的线性组合构成的，每个频率所对应的周期复

4、指数信号ej 1称为频率分量其相对幅度为对应频率的|X(j )|之值，其相位为对应频率的X(j )的相位。给定一个连续时间非周期信号x(t),它的频谱也是连续且非周期的。对于连续时间周期信号，也可以用傅里变换来表示其频谱，其特点是，连续时间周期信号的傅里叶变换时有冲激序列构成的，是离散的 一一这是连续时间周期信号的傅里叶变换的基本特征。3、连续周期信号的傅里叶级数CTFS的MATLAB 实现3.1傅里叶级数的 MATLAB 计算设周期信号x(t)的基本周期为Ti,且满足狄里克利条件，则其傅里叶级数的系数可由式9.4计算得到。式9.4重写如下：T1/2 ak x(t)e jk otdt T1 T

5、1/22基本频率为：0对周期信号进行分析时，我们往往只需对其在一个周期内进行分析即可，通常选择主周 期。假定x1(t)是x(t)中的主周期，则Ti/21jk 0ta kXi (t)edtT1 T1/2计算机不能计算无穷多个系数，所以我们假设需要计算的谐波次数为N,则总的系数个数为2N+1个。在确定了时间范围和时间变化的步长即Ti和dt之后，对某一个系数，上述系数的积分公式可以近似为：T1 /2ak x1(t)e jk 0tdtx(tn)e jk 0tdt/T)11Tl /2nx(t1),x(t2), x(tM) e jk 0t1,e jk 0t2, e jk 0tMl dt/T对于全部需要的2

6、N+1个系数，上面的计算可以按照矩阵运算实现。MATLAB实现系数计算的程序如下：dt = 0.01;T = 2; t = -T/2:dt:T/2; w0 = 2*pi/T;x1 = in put( 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;例1:给定一个周期为 T1 = 2s的连续时间周期方波信号，如图所示，其一个周期内的数学 表达式为：x1(t)% P

7、rogram2_1 clear, close all T = 2; dt = 0.00001;1,0t10,1t2t = -2:dt:2;图9.1周期方波信号x1 = u(t) - u(t-1-dt); x = 0;for m = -1:1x = x + u(t-m*T) - u(t-1-m*T-dt);endw0 = 2*pi/T;N = 10;L = 2*N+1;for k = -N: N;ak(N+1+k) = (1/T)*x1*exp(-j*k*w0*t)*dt;phi = anglel(ak);end% Evaluate the phase of ak执行程序Program2_1后，

8、就完成了信号的傅里叶级数的系数的计算，在命令窗口键入 ak命令窗口就可以显示傅里叶级数的21个系数:在程序Program2_1中添加相应的计算| ak |和绘图语句，就可以绘制出信号的幅度谱和 相位谱的谱线图。3.2周期信号的合成以及Gibbs现象从傅里叶级数的 合成式(Synthesis equation )jk 0tx(t)akek可以看出，用无穷多个不同频率和不同振幅的周期复指数信号可以合成一个周期信号。然而，我们无法用计算机实现对无穷多个周期复指数信号的合成。但是，用有限项来合成却是可行的，在实际应用中，多半也就是这么做的。然而，这样做的一个必然结果，就是引入了误差。如果一个周期信号在

9、一个周期有内断点存在，那么，引入的误差将除了产生纹波之外，还将在断点处产生幅度大约为9%的过冲(Overshot),这种现象被称为 吉伯斯现象(Gibbsphenomenon )。为了能够观察到合成信号与原信号的不同以及Gibbs现象，我们可以利用前面已经计算出的傅里叶级数的系数，计算出截短的傅里叶级数：Nx(t)akejk 0tkN这个计算可用L = 2N+1次循环来完成：c c , 、 j(r 1 N) 0tx2 x2 ak(r) e其中r作为循环次数，x2在循环之前应先清零。完成这一计算的MATLAB程序为：x2 = 0; L = 2*N+1;for r = 1:L;x2 = x2+ak

10、(r)*exp(j*(r-1-N)*w0*t);end;完成了所有的计算之后，就可以用绘图函数：plot()和stem()将计算结果包括x1, x2,abs(ak)和angle(ak)以图形的形式给出，便于我们观察。例2：修改程序Program2_1,使之能够用有限项级数合成例1所给的周期方波信号，并绘制出原始周期信号、合成的周期信号、信号的幅度谱和相位谱。% Program2_2% This program is used to compute the Fourier series coefficients ak of a periodic square waveclear,close al

11、lT = 2; dt = 0.00001; t = -2:dt:2;x1 = u(t)-u(t-1-dt); x = 0;for m = -1:1x = x + u(t-m*T) - u(t-1-m*T-dt);% Periodically extend x1(t) to form a periodicsignal end w0 = 2*pi/T;N = input(Type in the number of the harmonic components N =:);L=2*N+1;for k = -N:1:N;ak(N+1+k) = (1/T)*x1*exp(-j*k*w0*t)*dt;e

12、nd phi = angle(ak);y=0;for q = 1:L;% Synthesiz the periodic signal y(t) from the finite Fourier seriesy = y+ak(q)*exp(j*(-(L-1)/2+q-1)*2*pi*t/T);end;subplot(221),plot(t,x),title(The original signal x(t),subplot(223), plot(t,y),title(The synthesis signal y(t),subplot(222) k=-N:N;stem(k,abs(ak),k.),ax

13、is(-N,N,-0.1,0.6) subplot(224) stem(k,phi,r.),axis(-2,2,-0.2,1.2),axis(-2,2,-0.2,1.2), xlabel(Time t),title(The amplitude |ak| of x(t),title(The phase phi(k) of x(t),axis(-N,N,-2,2), xlabel(Index k)4用MATLAB 实现CTFT 计算MATLAB 进行傅里叶变换有两种方法，一种利用符号运算的方法计算，另一种是数值 计算，本实验要求采用数值计算的方法来进行傅里叶变换的计算。严格来说，用数值计算的方法计

14、算连续时间信号的傅里叶变换需要有个限定条件，即信号是时限信号(Time limitedsignal),也就是当时间|t|大于某个给定时间时其值衰减为零或接近于零，这个条件与前面提 到的为什么不能用无限多个谐波分量来合成周期信号的道理是一样的。计算机只能处理有限大小和有限数量的数。采用数值计算算法的理论依据是：X(j )xejtdtTm0 x(kT)e jk TTk若信号为时限信号，当时间间隔T取得足够小时，上式可演变为:X(j ) T x(kT)ek Njk Tx(t1),x(t2), ,x(t2N 1) e j tl,e jt2,e j t2N 1T上式用 MATLAB 表示为:X=x*ex

15、p(j*t *w)*T其中X为信号x(t)的傅里叶变换，w为频率Q, 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 )分别以图形的方式表现出来，并对图形加以适当的标注。三实验内容和要求1 ,n 、-sin()

16、cos( n t) n 1 n 21编写程序，绘制下面的信号的波形图:11x(t) cos( 0t)cos(3 0t) cos(5 t) L35其中，0 = 0.5兀，要求将一个图形窗口分割成四个子图，分别绘制cos( 0t)、cos(3 0t)、cos(5 0t)式中的项数。clear,close all, dt=0.00001;t=-2:dt:4;w0=0.5*pi;x1=cos(w0.*t);x2=cos(3*w0.*t);x3=cos(5*w0.*t);和x(t)的波形图，给图形加title,网格线和x坐标标签，并且程序能够接受从键盘输入的和N = input(Type in the

17、number of the harmonic components N =:);x=0for q=1:Nx=x+(sin(q*(pi/2).*cos(q*w0*t)/q;endsubplot(221)plot(t,x1) axis(-2,4,-2,2);grid on, title(signal cos(w0.*t) xlabel(the time t(sec),) subplot(222) plot(t,x2) axis(-2,4,-2,2);subplot(223)plot(t,x3) axis(-2,4,-2,2);grid on, title(signal cos(5*w0.*t) s

18、ubplot(224)plot(t,x)axis(-2,4,-2,2);grid on, title(signal x(t),) xlabel(the time t(sec),)grid on, title(signal cos(3*w0.*t) xlabel(the time t(sec),)the time t(sec)the time t(sec)the time t(sec)2给程序Program2_1增加适当的语句，并存盘，使之能够计算例1中的周期方波信号的傅里叶级数的系数，并绘制出信号的幅度谱和相位谱的谱线图。clear, close allendT = 2; dt = 0.000

19、01; t = -2:dt:2;w0 = 2*pi/T;x1 = u(t) - u(t-1-dt); x = 0;N = 10; % The number of the harmonic componentsfor m = -1:1% PeriodicallyL = 2*N+1;extend x1(t) to form a periodic signalfor k = -N: N; % Evaluate the Fourier seriesx = x + u(t-m*T) - u(t-1-m*T-dt);coefficients akaxis(-10,10,-2,2);endgrid on;p

20、hi = angle(ak);% Evaluate the phase of aktitle(xiangweipu)subplot(211);k=-10:10;stem(k,abs(ak),k);axis(-10,10,0,0.6);grid on;title(fudupu)subplot(212);k=-10:10;stem(k,angle(ak),k);3反复执行程序 Program2_2,每次执行该程序时，输入不同的N值，并观察所合成的周期方波信号。通过观察，描述总结吉伯斯现象的特点。N=10:N=40 :Tl u 产方 4勺门；*承4小C El-Ktau 4N=20:N=80:kak(

21、N+1+k) = (1/T)*x1*exp(-j*k*w0*t)*dt;4给定如下周期信号Xi(t)4.1仿照程序Program2_1 ,编写程序，以计算xi(t)的傅里叶级数的系数。clear,close allT = 2; dt = 0.00001; t = -3:dt:3;x=(t+1).*(u(t+1)-u(t)-(t-1).*(u(t)-u(t-1);x1=0;for m=-1:1x1=x1+(t+1-m*T).*(u(t+1-m*T)-u(t-m*T)-(t-1-m*T).*(u(t-m*T)-u(t-1-m*T);end w0 = 2*pi/T;N = input(Type in

22、 the number of the harmonic components N =:);L = 2*N+1;for k = -N: Nak(N+1+k) = (1/T)*x*exp(-j*k*w0*t)*dt; end y=0;for q = 1:Ly = y+ak(q)*exp(j*(-(L-1)/2+q-1)*2*pi*t/T);end;phi=angle(ak);akak =1 至6列0.0000+ 0.0000i 0.0025+ 0.0000i 0.0081 - 0.0000i7 至12列0.0000- 0.0000i 0.0225 -+ 0.0000i 0.2026 - 0.000

23、0i13 至18列0.0000- 0.0000i 0.0225-0.0000i 0.0041 + 0.0000i19 至21列-0.0000i 0.0000 -0.0000i 0.0000 + 0.0000i 0.0000 +0.0000i 0.0041 -0.0000i 0.2026 +0.0000i 0.0081 +0.0000i 0.00000.0000i 0.50000.0000i 0.00000.0000 + 0.0000i 0.0025 + 0.0000i 0.0000 - 0.0000ix1(t)的波形图，用有限项级4.2仿照程序Program2_2,编写程序，计算并绘制出原始信

24、号数合成的x1(t)的波形图，以及x1(t)的幅度频谱和相位频谱的谱线图。clear,close allT = 2; dt = 0.00001; t=-3:dt:3;title(fudupu)subplot(223);k=-N:N;stem(k,angle(ak),k);axis(-10,10,-2,2);grid on;title(xiangweipu);xlabel(index x);subplot(224);plot(t,y);axis(-4,4,-0.5,4);grid on;title(The signal y(t);xlabel(Time t(sec);ylabel(signal

25、y(t);x=(t+1).*(u(t+1)-u)-(t-1) .*(u(t)-u(t-1);x1=0;for m=-1:1x1=x1+(t+1-m*T).*(u(t+1-m*T)-u(t-m*T)-(t-1-m*T ).*(u(t-m*T)-u(t-1-m*T); endw0 = 2*pi/T;N = input(Type in the number of the harmonic components N =:);L = 2*N+1;for k = -N: Nak(N+1+k)=(1/T)*x*exp(-j*k*w0*t)* dt;y=0;for q = 1:Lyy+ak(q)*exp(j*

26、(-(L-1)/2+q-1)*2*pi*t/T);end;%phi=angle(ak);subplot(221);plot(t,x1);axis(-4,4,-0.5,1.2);grid on;title(The sinal x1(t);xlabel(Time (sec);ylabel(signal x1(t);subplot(222);k=-N:N;stem(k,abs(ak),k);axis(-10,10,0,2.0);endgrid on;10.5D-a.5The signal x1(t)-4-2024Time t(sec)frequency index kThe signal y(t)0.50事.5-4-2024Time t(sec)以下内容，可视个人实验进度决定是否选做5连续时间非周期信号的傅里叶变换给定一个时限信号： TOC o