第3章 Fourier变换.docx

第3章 Fourier变换在第一章，我们学习了几种不同频率的信号可以合成一个信号，如“拍”的复杂信号。反过来，能否从复杂的合成信号中分离出原来的信号呢？能，这就要用到我们本章要讲的Fourier变换。本章首先介绍Fourier级数及Fourier谱，然后介绍Fourier变换的性质及快速Fourier变换，最后介绍Fourier变换在数字滤波方面的应用。3.1 Fourier级数与Fourier变换在我们生活的世界里充满了各种各样的周期现象，如日常生活中常见的昼夜更替、四季循环、温度气压等气象因素的反复变化、河水水位的周期性涨落、地面植物的岁月枯荣以及科学技术中诸如太阳活动的十一年周期起伏、地磁场的日变化、重力仪上所反映的固体潮、交流电、电磁波和机械振动等无一不是周期现象。描述周期现象的最简单的周期函数是物理学上所说的谐波函数，它由正弦或余弦函数来表示： (3-1)利用三角公式,上式可以写成： （3-2）由于是常数，令a=Acos, 则可得： （3-3）这里 , (3-4)由此可以看出：一个带初相位的余弦函数可以看成一个不带相位的正弦函数与一个不带相位的余弦函数的合成。谐波函数是周期函数中最简单的函数，它描述的也是最简单的周期现象，在实际中所碰到的周期现象往往比它复杂得多。但这些复杂函数均在一定近似程度上可分解为不同频率的正弦函数和余弦函数。下面我们就介绍如何将一种复杂的函数分解为一系列不同频率的正弦函数和余弦函数的方法。3.1.1 周期函数的Fourier变换回想一下我们在数学中讲的Fourier级数：如何将一个周期为2l函数分解为Fourier级数呢？我们已经学过的Fourier级数展开式为 ： (3-5)其中 (3-6)如果f(x)是奇函数,积分上下限相互对称，则这时亦为奇函数，故an均为零，得到的Fourier级数是正弦级数: (3-7)其中，bn的积分可简写为： (3-8)如果f (x)是偶函数, 因积分上下限相互对称，并且为奇函数，故bn均为零,得到的Fourier级数是余弦级数： (3-9)其中a0和an可简写为： (3-10)312 离散Fourier 变换 现在我们设法把上述公式不加证明地应用于离散Fourier级数中。我们在数字资料处理中经常遇到的不是一个函数，而是一个离散的数列，比如说，等间隔时间取样的时间序列(这里的数据个数为N,一般取N为偶数,并且若取一特殊的偶数还可以使计算速度加快)。因此用上述公式时需要加以改造。首先，我们得到的数字信号只能在正的时间段取值,在负的时间段不能取值。但由于我们认为所取的时间是无限长的周期序列，周期为2l，因此，我们可以把（-l，l）修改为（0,2l）这样就避免了在负时间段取值。x就对应于时间序列的时间，即，i=0,1,2,N-1。由于我们要处理的是离散的数据序列，因此不能再用积分，而应用积分的离散形式-求和来表示，即：我们在（0，2 l）里等间隔地取N个取值点，取样时间间隔为，其中。 有了上述改正，我们有 （3-11） （3-12）所以（3-5）公式的离散形式为： （3-13）式中 （3-14） （3-15）现在考虑（3-13）式，若我们共有N个观测数据，即,则上述方程可以列出N个方程式，式（3-14）和式（3-15）为方程组的解。那么，我们考虑一下m最大可以取多大呢？在（3-13）式中均是未知数，m越大，未知数的个数越多。（3-13）式所列出的N个方程最多只能确定N个未知数，可见m最大只能取。这里需要注意的是，求和号里面的未知数的个数已经有N个了，再加上a0这个未知数，未知数的个数是否多于方程的个数呢？如果我们分析一下求和号里面的,当k取时,最后一项为:，所以根本没必要求。但这里需要注意的是，与a0一样，这里的也只取一半。这样，取m最大为就正好满足了N个方程确定N个未知数的唯一解的条件。对于N为奇数的情况，m的取值应向下取整得到m=floor(N/2),此时的不为零，按正常计算。xi的第k项为：，我们将其写成以下形式： （3-16）由上式可以看出，第k项为两个周期函数之和，一个为正弦函数，一个为余弦函数，它们的频率均为： （3-17）这里的为时间,圆频率为。T为所取数据的总的时间长度，随着k的增大，三角函数的频率逐渐增加，周期逐渐减小。它们的周期为： （3-18）即为整个记录长度的。因此，对于一个有限长观测数据序列来说，使用Fourier级数法求得的各种周期或频率的波是有限制的。因为k只能取0 的整数。k=0时，a0为该数据序列的直流分量。k=1时，波的周期为包含全部观测数据记录的时间T，这是波的最大周期。其频率为，是用Fourier谐波分析所能检测到的最小频率。通常我们称k=1的波为基波，基波的周期由取样数据的长度来决定。k取最大值时，波的周期为。此时波的周期最小，频率最大，为。在数字信号处理中称这个最大频率为Nyquist频率，这与第二章所讲的取样定理（采样定理）是一致的。Nyquist频率由取样间隔来决定，为取样间隔2倍的倒数。在中国数字地震台网系统中，每秒取50个样，取样间隔为0.02秒，因此能够得到地震波的最大频率为25赫兹。如果想得到高频率的地震波，必须提高取样频率，即缩短取样间隔。我们称利用（3-1315）式分离的各种频率（周期）的波，即（3-16）式为k次谐波。k次谐波的频率可由(3-17)式给出，周期可由（3-18）式给出。若我们将（3-13）式中的求和号里的各次谐波写成（3-2）式的形式，则有： （3-19）式中 ， （3-20）ck表示了k次谐波的振幅大小。在分解地震波时，哪些分量的振幅大，哪些分量的振幅小，对地震波的性质有十分重要的影响。在研究地震波对工程结构的影响时，因为不同的工程构件有不同的固有频率，地震波中与工程构件的固有频率相等或相近的谐波振幅很大时，会使该构件发生共振而造成破坏。我们通常用频率作为横坐标，用振幅作为纵坐标绘出图形来研究频率和振幅的关系，这种图形就是一种重要的谱，通常称为Fourier振幅谱，简称振幅谱(振幅 spectra)。为k次谐波的初相，将（3-20）式中波的各分量的初相按频率画出，与振幅谱对应，我们称之为相位谱（Phase spectra）。如果将转换为我们平常看到的角度，只需将乘以即可。(3-19) 式和（3-20）式之所以重要是因为它将一个数字序列转化为一系列不同频率的波的合成。这样作的好处是通过这种分解，我们可以看出所给波列中含有哪种频率成分。如果我们想滤掉某种频率的波，直接去除该频率的谐波即可。下面我们根据上式给出一个合成波分解为各种不同频率波的例子。313信号的Fourier分解与合成举例【例3-1】将振幅为1的1Hz正弦波和振幅为0.5的5Hz正弦波相加后进行Fourier分析，研究能否从中分析出含有这两种频率的信号。%Samp3_1clear all %将工作空间中的所有变量清除N=256;dt=0.02; %数据的个数和采样间隔n=0:N-1;t=n*dt; %序号序列和时间序列x=sin(2*pi*t)+0.5*sin(2*pi*5*t); %信号加得到的合成信号m=floor(N/2)+1; %分解a,b的最大序号值，为分解的N/2个参数再加参数a0%floor函数为向下取整a=zeros(1,m);b=zeros(1,m); %产生a,b两个为零的序列for k=0:m-1 for ii=0:N-1 a(k+1)=a(k+1)+2/N*x(ii+1)*cos(2*pi*k*ii/N); % b(k+1)=b(k+1)+2/N*x(ii+1)*sin(2*pi*k*ii/N); % %MATLAB中的数组序号只能从1开始。 end c(k+1)=sqrt(a(k+1).2+b(k+1).2); %endsubplot(2,1,1),plot(t,x);title(原始信号),xlabel(时间/s) %绘出时间域信号subplot(2,1,2),plot(0:m-1)/(N*dt),c) %绘出频率域信号，对应频率点用（3-17）给出title(Fourier变换),xlabel(频率/Hz),ylabel(振幅)图3-1 将1Hz和5Hz的正弦振动合成后进行Fourier分析的结果由图3-1可以清楚地看出，从初始信号中明确地识别出了1Hz和5Hz的波。这里1Hz和5Hz的振幅与原来信号振幅并不完全一致，这是由于数据采样点较少造成的。自己作实验可以看出，数据采样点越多，Fourier分析得到的结果与原始振幅越接近。如果将m修改为N,其频谱只是增加了镜像对称的部分，如果将m改为2*N，则频谱增加一个周期（见图3-2），改为3*N，频谱就增加了三个周期，可以明确地看出频谱的周期性。自己修改上面的程序可以验证这一点(参看程序Samp3_1_1)。由此可以看出，离散有限信号的频谱为周期谱。图3-2 取大于Nyquist频率时振幅谱的例子图示。上图，Nyquist频率的2倍；下图，Nyquist的4倍上面的例子表明，k值取不小于N/2时，可以完全分辨出小于采样频率一半（Nyquist频率）的频率，反过来，根据分解成的ak、bk值能否准确地恢复原来的信号呢？下面我们采用（3-10）式或（3-15）式求得上例合成信号，在前面程序运行后，运行下面的程序：%Samp3_1_2if(mod(N,2)=1)a(m)=a(m)/2; end %此时b(m)为零，a(m)减半for ii=0:N-1 xx(ii+1)=a(1)/2; for k=1:m-1 xx(ii+1)=xx(ii+1)+a(k+1)*cos(2*pi*k*ii/N)+b(k+1)*sin(2*pi*k*ii/N);% end endsubplot(2,1,1)plot(0:N-1)*dt,x) %绘制原信号title(原始信号)subplot(2,1,2)plot(0:N-1)*dt,xx) %绘制合成信号title(合成信号),xlabel(时间/s)图3-3 采用1Hz和5Hz合成振动及其Fourier分解系数合成振动的比较。运行结果为图3-3，可以看到，运用Fourier分解得到的系数合成后的振动和原振动完全一致。验证了上述理论的正确性。我们从另一个角度来考虑问题，如果k取值不到floor(N/2)+1，所合成的波形与原波形有何区别呢？下面以方波为例来说明这个问题。【例3-2 】 用振幅为0.8的方波进行Fourier分析，并用分析得到的系数求解当k取不同值时的合成信号图。程序如下：%Samp3_2clear all %清除内存所有变量close all %关闭所有打开的图形窗口N=200;dt=4/N; %数据点数和采样间隔for n=1:N %得到方波信号 if (n*dt=2) x(n)=0.8; else x(n)=-0.8; endendfigure(1) %打开第一个图形窗口，绘制原始信号及其Fourier分析系数subplot(2,1,1),plot(1:N)*dt,x),hold on;plot(1:N)*dt,zeros(1,N),k),xlabel(时间/s)title(原始信号)a=zeros(1,N);b=zeros(1,N);nn=floor(N/2)+1;for k=0:nn-1 a(k+1)=0; b(k+1)=0; for ii=0:N-1 a(k+1)=a(k+1)+2/N*x(ii+1)*cos(2*pi*k*ii/N); %求解Fourier系数 b(k+1)=b(k+1)+2/N*x(ii+1)*sin(2*pi*k*ii/N); end c(k+1)=sqrt(a(k+1).2+b(k+1).2);endsubplot(2,1,2),plot(0:nn-1)/(N*dt),c);title(Fourier变换)xlabel(频率/Hz),ylabel(振幅) %绘制振幅谱m=input(输入谐波最大阶数?); %输入最大k值if(m(floor(N/2)+1) %如果最大k值大于Nyquist频率对应的点数，则显示出错信息 error(谐波最大阶数必须小于Nyquist频率对应的阶数.)endif(mod(N,2)=1)a(nn)=a(nn)/2; end %此时b(nn)为零，a(nn)减半for ii=0:N-1 %合成信号 xx(ii+1)=a(1)/2; for k=1:m xx(ii+1)=xx(ii+1)+a(k+1)*cos(2*pi*k*ii/N)+b(k+1)*sin(2*pi*k*ii/N); endendfigure(2) %打开第二个图形窗口，绘制合成信号图plot(1:N)*dt,xx,(0:N-1)*dt,x) %绘制合成信号和方波图形便于比较hold onplot(1:N)*dt,zeros(1,N),k),xlabel(时间/s) %绘出横轴title(合成信号)程序运行后得到的第一个图为图3-4。该图表明方波信号进行Fourier分解后，得到的系数几乎分布在整个频率轴上。即说明方波可以分解为含有多种频率的振动。从图中还可以发现，低频成分的振幅较大，随着频率的增高，振幅逐渐衰减。图3-4 方波及其Fourier分解系数图示。当m=1,3,5,9时(运行程序Samp3_2_1)，得到的合成信号组合在一起为图3-5，可以看出，m越大，即k值取得越多，与原始的方波信号越接近。当m=1，只将方波分解为一个正弦振动；当m=3时，在原来m=1振动的基础上又叠加了一次谐波，使得该振动波形与方波更加接近；当m=5时，得到的振动波形与方波更加接近，以此类推。可以想象，当m为floor(N/2)+1时，合成波形与原波形完全一致。请读者自己修改程序观察合成的效果。这里应当注意的是，本例中的方波信号为奇函数，根据（3-6）式，只有bk在起作用。图3-5 当m=1,3,5,7时得到的合成振动与原始方波信号的比较通过前面的Fourier级数分析，我们知道采用Fourier级数分解技术可以得到信号中含有哪种频率成分，振幅为多少。并且运用分解得到的Fourier系数可以合成为原振动图像，这是Fourier分析的基础，我们必须掌握。3.2 复数形式的Fourier级数及其应用3.2.1 基 本 理 论前面我们已经得出了离散Fourier变换的 (3-13) (3-15)式。如果我们将和（即分解成的k次谐波余弦和正弦函数的振幅）表示为一个复数实部和虚部的某种函数的组合，并且将、也表示为一个复数的实部和虚部的某种组合，我们可得到下式（可用欧拉公式证明，在此从略）： , i=0,1,2,N-1 （3-21） , k=0,1,2,N-1 （3-22）注意，这里的ck为复数序列，其模序列对应于Fourier变换所分解的各谐波的振幅，其频率由确定，即，为采样间隔。若我们假定采样间隔为1秒，则频率为。其初相由实部和虚部的相对大小来确定。这样我们就将一个信号分解为多种频率的信号。若无某种频率的信号，则该频率信号的振幅为零。另外，考虑到欧拉公式，当k由0 N-1变化时，后半部分复序列是前半部分复序列的共轭（即实部相同，虚部正负号相反）。这样Fourier变换（3-22）式得到的复序列也具有后半部分复序列是前半部分复序列的共轭的性质。由于后半部分复序列是前半部分复序列的共轭，因此，我们将k取到N-1，并没有增加独立参数的个数。这种形式是Fourier变换的最基本的形式，以后所讲得Fourier变换都是指的这种形式。这样，就建立了复数形式Fourier变换之间的关系。如果已知时间域内的等间隔数据，要求频率域内的振幅，采用（3-22）式，称之为Fourier变换。已知频率域内的情况求取时间域的值，我们采用（3-21）式，称之为Fourier逆变换。注意：在作Fourier变换时，所求的结果为复数形式。我们求其模，就得到振幅谱的相对值。真正的振幅大小与ck的关系为： （3-23）相位的计算方法为： （3-24）3.2.2 Fourier变换程序为了运用上面的Fourier变换公式得到MATLAB中的子程序，我们首先解剖一个例子。【例3-3 】求周期序列x(n)的离散Fourier变换，其中x(n)=0 1 2 3；解： 基本周期为N=4，令W4=，则Fourier变换公式可写为 ，将数据代入（3-22）式可得：X(0)=6 X(1)=-2+2j X(2)=-2 X(3)=-2-2j根据该例的解法，可以将Fourier变换写成下面的程序function Xk = dfs(xn,N)% 计算离散Fourier序列的系数% 调用方式：Xk = dfs(xn,N)% Xk：离散Fourier变换的系数，序号：0 = k = N-1% xn：取周期信号的一个周期，序号：0 = n = N-1% N： xn的数据个数n = 0:1:N-1; % 数据数组序列k = 0:1:N-1; % 频率数组序列WN = exp(-j*2*pi/N); % Wn因子nk = n*k; % 产生N*N的矩阵nkWNnk = WN . nk; %离散Fourier变换矩阵Xk = xn * WNnk; % Fourier变换系数这里巧妙地利用了MATLAB中的矩阵相乘的概念使程序简化，式中，对于例3-3WNnk=WNnk也是一个4*4的矩阵。利用矩阵xn=0 1 2 3与矩阵WNnk相乘中的相加，即，相加的运算蕴含在矩阵相乘的运算中，请大家注意MATLAB中的这种简便方法。这里有一个问题，即得出的Fourier变换与第一节所讲的变换差一个系数，即所得的结果要乘以2/N才为某种频率的真实振幅，因此如想得到真实振幅，必须将结果乘以2/N。Fourier逆变换可以用下列程序表示，其中分析原理与前面相同。function xn = idfs(Xk,N)%计算离散Fourier逆变换% 调用方式：xn = idfs(Xk,N)% xn：输出的逆Fourier逆变换后的序列，序号：0 = n = N-1% Xk：输入数组，序号：0 = k = N-1% N：Xk的元素个数 %n = 0:1:N-1; % 数据数组序列号k = 0:1:N-1; %频率数组序列号WN = exp(-j*2*pi/N); % Wn因子nk = n*k; %产生NN的矩阵nkWNnk = WN . (-nk); %离散Fourier逆变换矩阵xn = (Xk * WNnk)/N; % Fourier逆变换后的数组值【例3-4】 重新用程序求解上例的Fourier变换xn=0 1 2 3;N=4;Xk=dfs(xn,N)可以发现得到的结果与上面结果一致，也进一步证明了程序的正确性。3.2.3 应 用 实 例前面我们对Fourier变换进行了理论分析，并给出了程序，这里我们给出Fourier变换的实例。在第一章中，我们了解到振幅相同、频率相差不大的两个振动可以合成为一种称为“拍”的、振幅周期性变化的振动。那么从拍的数据中通过Fourier变换能否分析出其中的频率成分呢？请看下面的例子。【例3-5】已知序列x(n)=cos(2*pi*0.24*n)+cos(2*pi*0.26*n), n=0:99,试绘制x(n)及其傅立叶变换的幅值图，并用变换后的数值求解逆变换。其中采样频率为1s。我们知道，假定采样频率为1Hz（即以1秒的采样间隔进行采样）。则该信号包含0.24Hz和0.26Hz两种频率的波，其振幅均为1。%Samp3_5clf %清除图形框的内容N=100;dt=1; %设置最大点数n=0:N-1; t=n*dt; %给出时间序列xn=cos(2*pi*0.24*t)+cos(2*pi*0.26*t); %给出原始信号的值序列Xk=dfs(xn,N); %对原始信号进行Fourier变换magXk=abs(Xk);phaXk=angle(Xk); %求出Fourier变换的振幅和相位subplot(2,2,1),plot(t,xn); xlabel(时间/s) %绘出原始信号title(原始信号(N=100);xx=idfs(Xk,N); %Fourier逆变换x=real(xx); %取变换后的实部，如做实验可以验证其虚部为零subplot(2,2,2),plot(t,x),xlabel(时间/s),title(运用Fourier逆变换得到的合成信号)k=0:length(magXk)-1;subplot(2,1,2),plot(k/(N*dt),magXk*2/N); %绘出Fourier变换的振幅谱，频率采用（3-17）%采用真实振幅（3-23）式绘图xlabel(频率/Hz);ylabel(振幅);title( X(k)振幅(N=100);程序的运行结果见图3-6。程序中函数abs(x)用于计算复向量x的幅值，angle(x) 用于计算复向量x的相角，介于-，用弧度表示。可以看出，Fourier变换后确实分析出频率为0.24Hz和0.26Hz的振动。而其他频率成分的幅值为零，表明信号中不存在其他频率成分。运用逆变换完全恢复了原始信号，表明了该程序的正确性。由于前半部分复序列是后半部分复序列的共轭，因此，Fourier变换后的振幅谱以0.5Hz(k=50)为对称轴将振幅序列分为对称的两部分。本例中我们没有给出信号的相位谱。要想得到相位谱，只要在程序中加以修改即可，请大家自己修改程序完成。图3-6 对0.24Hz和0.26Hz组成的拍信号进行Fourier变换，及与逆变换得到结果的比较【例3-6】对于上例中的x(n)通过补零增长至120个数,绘制x(n)和它的离散Fourier幅值图，将原始信号与逆变换恢复的结果进行比较。该例只要将例3-5程序中调用dfs函数之前的语句置换为下面的语句，并改变后面的标记即可进行分析：%Samp3_6clf %清除图形框的内容N=120;dt=1; %设置最大点数n=0:N-1; t=n*dt; %给出时间序列xn=cos(2*pi*0.24*0:99)+cos(2*pi*0.26*0:99); xn=xn,zeros(1,N-100);%给出原始信号的值序列图3-7对0.24Hz和0.26Hz组成的拍信号加20个零后进行Fourier变换并与逆变换得到结果的比较运行结果为图3-7。可以看到，上例中的“拍”信号加上20个零后，其Fourier变换分解的频率成分除了含有原有的0.24Hz和0.26Hz的频率信号外，还含有很多频率的小振幅振动，也就是说，多种频率的小振幅振动叠加才能使后面的时间域信号为零。3.3 Fourier变换的性质前面介绍了Fourier变换的基本理论，为了进一步了解Fourier变换，本节对Fourier变换的性质进行讲解。3.3.1 线 性 性若x(n)和y(n)分别具有离散Fourier变换X(k)和Y(k)，则x(n)+y(n)的Fourier变换为X(k)+Y(k)。此定理通过公式可立即得出。下面的例子简单说明了Fourier变换的线性性。【例3-7】设有两个振动cos(2*pi*0.24*t)和cos(2*pi*0.12*t)。求两个振动及其合成振动的Fourier变换，采样间隔为1s，数据长度为100。%Samp3_7clfN=100; dt=1; %数据