1、指数形式的傅里叶级数告诉我们，如果一个周期信号x（t），满足狄里克利条件，那么，它就可以被看作是由很多不同频率的互为谐波关系（harmonically related）的周期复指数信号所组成，其中每一个不同频率的周期复指数信号称为基本频率分量，其复幅度（complex amplitude）为这里“复幅度（complex amplitude）”指的是通常是复数。上面的傅里叶级数的合成式说明，我们可以用无穷多个不同频率的周期复指数信号来合成任意一个周期信号。然而，用计算机（或任何其它设备）合成一个周期信号，显然不可能做到用无限多个谐波来合成，只能取这些有限个谐波分量来近似合成。假设谐波项数为N，则
2、上面的和成式为：显然，N越大，所选项数越多，有限项级数合成的结果越逼近原信号x（t）。本实验可以比较直观地了解傅里叶级数的物理意义，并观察到级数中各频率分量对波形的影响包括“Gibbs”现象：即信号在不连续点附近存在一个幅度大约为9%的过冲，且所选谐波次数越多，过冲点越向不连续点靠近。这一现象在观察周期矩形波信号和周期锯齿波信号时可以看得很清楚。2、连续时间信号傅里叶变换-CTFT 傅里叶变换在信号分析中具有非常重要的意义，它主要是用来进行信号的频谱分析的。傅里叶变换和其逆变换定义如下： 连续时间傅里叶变换主要用来描述连续时间非周期信号的频谱。按照教材中的说法，任意非周期信号，如果满足狄里克利
3、条件，那么，它可以被看作是由无穷多个不同频率（这些频率都是非常的接近）的周期复指数信号ejt的线性组合构成的，每个频率所对应的周期复指数信号ejt称为频率分量（frequency component），其相对幅度为对应频率的|X（j）|之值，其相位为对应频率的X（j）的相位。X（j）通常为关于的复函数，可以按照复数的极坐标表示方法表示为：X（j）=| X（j）|ej X（j）其中，| X（j）|称为x（t）的幅度谱，而X（j）则称为x（t）的相位谱。 给定一个连续时间非周期信号x（t），它的频谱也是连续且非周期的。对于连续时间周期信号，也可以用傅里变换来表示其频谱，其特点是，连续时间周期信号的
4、傅里叶变换时有冲激序列构成的，是离散的这是连续时间周期信号的傅里叶变换的基本特征。3、离散时间序列的傅里叶变换-DTFT给定一个非周期离散时间序列xn，它的傅里叶变换定义为其反变换定义为 式称为离散时间傅里叶正变换，式称为离散时间傅里叶反变换。由式可以看出，一个非周期离散时间序列，总是可以被看作是由无穷多个不同频率的加权的基本频率分量ejn组合而成的。对序列中的频率为的频率分量来说，其权为X（ej），通常是复数，也可以将它表示为 X（ej） = | X（ej）|ejX（ej）|X（ej）|称为序列的幅度谱，而X（ej）称为序列的相位谱，它们都是频率的周期函数。4、连续周期信号的傅里叶级数CTF
5、S的MATLAB实现 傅里叶级数的MATLAB计算设周期信号x（t）的基本周期为T1，且满足狄里克利条件，则其傅里叶级数的系数可由式计算得到。式重写如下：基本频率为：对周期信号进行分析时，我们往往只需对其在一个周期内进行分析即可，通常选择主周期（Principle period）。假定x1（t）是x（t）中的主周期，则计算机不能计算无穷多个系数，所以我们假设需要计算的谐波次数为N，则总的系数个数为2N+1个。在确定了时间范围和时间变化的步长即T1和dt之后，对某一个系数，上述系数的积分公式可以近似为：对于全部需要的2N+1个系数，上面的计算可以按照矩阵运算实现。MATLAB实现系数计算的程序如
6、下：dt = ;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越小，精度越高，但是，计算机计算所花的时间越长。例题2-1：给定一个周期为T1 = 2s的连续时间周期方波信号，如图所示，其一个
7、周期内的数学表达式为：解：首先，我们根据前面所给出的公式，计算该信号的傅里叶级数的系数。因为：0 = 2/T1 = ，代入上式得到：在MATLAB命令窗口，依次键入： k = -10:10; ak = （-j）.k）.* （sin（k+eps）*pi/2）./（k+eps）*pi） % The expression of akak = Columns 1 through 4 0 + 0 + Columns 5 through 8 Columns 9 through 12 0 + 0 - Columns 13 through 16 0 - 0 - Columns 17 through 20 Co
8、lumn 21 从MATLAB命令窗口，我们得到了该周期信号从到共21个系数。紧接着再键入以下命令： subplot（221） stem（k,abs（ak）,k.） title（The Fourier series coefficients xlabel（Frequency index k就得到一幅如右图所示的描述与k之间的关系的图形。以上是我们通过手工计算得到的这个周期信号的傅里叶级数表达式及其频谱图，下面给出完成傅里叶级数系数计算的相应MATLAB范例程序。% Program2_1% This program is used to evaluate the Fourier series c
9、oefficients ak of a periodic square waveclear, close all dt = ; t = -2:2;x1 = u（t） - u（t-1-dt）; x = 0;for m = -1:1 % Periodically extend x1（t） to form a periodic signal x = x + u（t-m*T） - u（t-1-m*T-dt）;endw0 = 2*pi/T;N = 10; % The number of the harmonic componentsL = 2*N+1;for k = -N: N; % Evaluate
10、the Fourier series coefficients ak ak（N+1+k） = （1/T）*x1*exp（-j*k*w0*t）*dt;phi = anglel（ak）; % Evaluate the phase of ak执行程序Program2_1后，就完成了信号的傅里叶级数的系数的计算，在命令窗口键入 ak命令窗口就可以显示傅里叶级数的21个系数： + + - + - + - + - + - + - + - - 将这里的ak之值同前面手工计算得到的ak比较，可见两者是完全相同的。再次特别提示：程序中，时间变量的变化步长dt的大小对傅里叶级数系数的计算精度的影响非常大，dt越小
11、，精度越高，本程序中的dt之所以选择就是为了提高计算精度。但是，计算机所花的计算时间越长。在程序Program2_1中添加相应的计算| ak |和绘图语句，就可以绘制出信号的幅度谱和相位谱的谱线图。 周期信号的合成以及Gibbs现象从傅里叶级数的合成式（Synthesis equation）可以看出，用无穷多个不同频率和不同振幅的周期复指数信号可以合成一个周期信号。然而，我们无法用计算机实现对无穷多个周期复指数信号的合成。但是，用有限项来合成却是可行的，在实际应用中，多半也就是这么做的。然而，这样做的一个必然结果，就是引入了误差。如果一个周期信号在一个周期有内断点存在，那么，引入的误差将除了产
12、生纹波之外，还将在断点处产生幅度大约为9%的过冲（Overshot），这种现象被称为吉伯斯现象（Gibbs phenomenon）。为了能够观察到合成信号与原信号的不同以及Gibbs现象，我们可以利用前面已经计算出的傅里叶级数的系数，计算出截短的傅里叶级数：这个计算可用L = 2N+1次循环来完成：其中r作为循环次数，x2在循环之前应先清零。完成这一计算的MATLAB程序为：x2 = 0;for r = 1:L; x2 = x2+ak（r）*exp（j*（r-1-N）*w0*t）;end;完成了所有的计算之后，就可以用绘图函数：plot（）和stem（）将计算结果包括x1, x2, abs（a
13、k）和angle（ak）以图形的形式给出，便于我们观察。观察吉伯斯现象的最好的周期信号就是图2-1所示的周期方波信号，这种信号在一个周期内有两个断点，用有限项级数合成这个信号时，吉伯斯现象的特征非常明显，便于观察。例题2-2：修改程序Program2_1，使之能够用有限项级数合成例题2-1所给的周期方波信号，并绘制出原始周期信号、合成的周期信号、信号的幅度谱和相位谱。为此，只要将前述的for循环程序段和绘图程序段添加到程序Program2_1中即可，范例程序如下：% Program2_2% This program is used to compute the Fourier series c
14、oefficients ak of a periodic square waveclear,close allx1 = u（t）-u（t-1-dt）;1 % Periodically extend x1（t） to form a periodic signalN = input（Type in the number of the harmonic components N = :）;1:phi = angle（ak）;y=0;for q = 1: % Synthesiz the periodic signal y（t） from the finite Fourier series y = y+
15、ak（q）*exp（j*（-（L-1）/2+q-1）*2*pi*t/T）;subplot（221）,plot（t,x）, title（The original signal x（t）, axis（-2,2,）,subplot（223）, plot（t,y）, title（The synthesis signal y（t）, axis（-2,2,）, xlabel（Time t）,subplot（222）k=-N:）, title（The amplitude |ak| of x（t）, axis（-N,N,）subplot（224）stem（k,phi,r.The phase phi（k） of
16、 x（t）, axis（-N,N,-2,2）, xlabel（Index k在用这个程序观察吉伯斯现象时，可以反复执行该程序，每次执行时，输入不同之N值，比较所的图形的区别，由此可以观察到吉伯斯现象的特征。5 用MATLAB实现CTFT及其逆变换的计算 用MATLAB实现CTFT的计算MATLAB进行傅里叶变换有两种方法，一种利用符号运算的方法计算，另一种是数值计算，本实验要求采用数值计算的方法来进行傅里叶变换的计算。严格来说，用数值计算的方法计算连续时间信号的傅里叶变换需要有个限定条件，即信号是时限信号（Time limited signal），也就是当时间|t|大于某个给定时间时其值衰减为
17、零或接近于零，这个条件与前面提到的为什么不能用无限多个谐波分量来合成周期信号的道理是一样的。计算机只能处理有限大小和有限数量的数。采用数值计算算法的理论依据是：若信号为时限信号，当时间间隔T取得足够小时，上式可演变为：上式用MATLAB表示为： X=x*exp（j*t*w）*T其中X为信号x（t）的傅里叶变换，w为频率，T为时间步长。相应的MATLAB程序：T = ; dw = ; %时间和频率变化的步长t = -10:T:w = -4*pi:dw:4*pi;X（j）可以按照下面的矩阵运算来进行：X=x*exp（-j*t*）*T; %傅里叶变换X1=abs（X）; %计算幅度谱phai=ang
18、le（X）; %计算相位谱为了使计算结果能够直观地表现出来，还需要用绘图函数将时间信号x（t），信号的幅度谱|X（j）|和相位谱 X（j）分别以图形的方式表现出来，并对图形加以适当的标注。在Matlab中有FFT函数可以直接计算傅里叶变换，该函数使用了快速傅里叶变换算法（参考数字信号处理课程）。但该函数计算结果将负频率部分平移到了正频率结果的后面，因此需要使用FFTSHIFT函数将负频率谱数据移到它原来该处的位置。请看下面的演示示例。例：sin（2*pi*1*t）函数的频谱。（这里为了清楚展示频谱平移的情况，采用了离散谱的画法）程序：f=1;fs=4*f;t=0:1/fs:1;x=sin（2*
19、pi*t*f）;%注意x的长度是5y=fft（x）figure（1）subplot（311）stem（abs（y）;z=fftshift（y）;subplot（312）stem（abs（z）;faxis=linspace（-fs/2,fs/2,length（z）;subplot（313）stem（faxis,abs（z）;结果图形：第三幅图加上了频率坐标，请注意频率坐标向量的产生方法。 用MATLAB实现傅里叶逆变换 连续时间傅里叶逆变换可用式进行计算。从定义式可看出，其计算方法与傅里叶变换是一样的，因此可以采用同样的矩阵运算的方法来计算，即 x（t）=X（j）*exp（j*t）*d具体的MA
20、TLAB函数如下：t = -5:;5; % 指定信号的时间范围，此范围应根据信号的持续时间确定。dw = ; w = -4*pi:d:X = input（Type in the expression of X（jw）;x = X* exp（jw*t）*dw;然后用绘图函数就可以绘制出逆变换得到的时域信号波形图。 同样，在Matlab中也有IFFT函数用来求逆傅里叶变换，此函数一般用来对FFT函数计算出的傅里叶变换向量做逆变换，无需使用fftshift。计算sin（2*pi*10t）的频谱，并做逆变换。f=10;fs=100*f;%x5plot（t,x）;iy=real（ifft（y）;axis
21、（-40 40 0 max（abs（z）plot（t,iy）;6 连续时间LTI系统的频率响应所谓频率特性，也称为频率响应特性，简称频率响应（Frequency response），是指系统在正弦信号激励下的稳态响应随频率变化的情况，包括响应的幅度随频率的变化情况和响应的相位随频率的变化情况两个方面。连续时间LTI系统的时域及频域分析图 上图中x（t）、y（t）分别为系统的时域激励信号和响应信号，h（t）是系统的单位冲激响应，它们三者之间的关系为：，由傅里叶变换的时域卷积定理可得到： 或者：为系统的频域数学模型，它实际上就是系统的单位冲激响应h（t）的傅里叶变换。即 由于H（j）实际上是系统单
22、位冲激响应h（t）的傅里叶变换，如果h（t）是收敛的，或者说是绝对可积（Absolutly integrabel）的话，那么H（j）一定存在，而且H（j）通常是复数，因此，也可以表示成复数的不同表达形式。在研究系统的频率响应时，更多的是把它表示成极坐标形式：上式中，称为幅度频率相应（Magnitude response），反映信号经过系统之后，信号各频率分量的幅度发生变化的情况，称为相位特性（Phase response），反映信号经过系统后，信号各频率分量在相位上发生变换的情况。和都是频率的函数。对于一个系统，其频率响应为H（j），其幅度响应和相位响应分别为，如果作用于系统的信号为，则其响应
23、信号为若输入信号为正弦信号，即x（t） = sin（0t），则系统响应为可见，系统对某一频率分量的影响表现为两个方面，一是信号的幅度要被加权，二是信号的相位要被移相。由于都是频率的函数，所以，系统对不同频率的频率分量造成的幅度和相位上的影响是不同的。7 用MATLAB计算系统频率响应在本实验中，表示系统的方法仍然是用系统函数分子和分母多项式系数行向量来表示。实验中用到的MATLAB函数如下：H,w = freqs（b,a）：b,a分别为连续时间LTI系统的微分方程右边的和左边的系数向量（Coefficients vector），返回的频率响应在各频率点的样点值（复数）存放在H中，系统默认的样点
24、数目为200点；Hm = abs（H）：求模数，即进行运算，求得系统的幅度频率响应，返回值存于Hm之中。real（H）：求H的实部；imag（H）：求H的虚部；phi = atan（-imag（H）./（real（H）+eps）：求相位频率相应特性，atan（）用来计算反正切值；或者phi = angle（H）：求相位频率相应特性；tao = grpdelay（num,den,w）：计算系统的相位频率响应所对应的群延时。计算频率响应的函数freqs（）的另一种形式是：H = freqs（b,a,w）：在指定的频率范围内计算系统的频率响应特性。在使用这种形式的freqs/freqz函数时，要在前
25、面先指定频率变量w的范围。例如在语句H = freqs（b,a,w）之前加上语句：w = 0:2*pi/256:2*pi。下面举例说明如何利用上述函数计算并绘制系统频率响应特性曲线的编程方法。假设给定一个连续时间LTI系统，下面的微分方程描述其输入输出之间的关系编写的MATLAB范例程序，绘制系统的幅度响应特性、相位响应特性、频率响应的实部和频率响应的虚部。程序如下：% Program2_3% This Program is used to compute and draw the plots of the frequency response % of a continuous-time systemb = 1; % The coefficient vector of the right side of the differential equationa = 1 3 2; % The c
