第3章Fourier变换

第3章Fourier变换在第一章，我们学习了几种不同频率的信号可以合成一个信号，如“拍”的复杂信号。反过来，能否从复杂的合成信号中分离出原来的信号呢？能，这就要用到我们本章要讲的Fourier变换。本章首先介绍Fourier级数及Fourier谱，然后介绍Fourier变换的性质及快速Fourier变换，最后介绍Fourier变换在数字滤波方面的应用。3.1Fourier级数与Fourier变换在我们生活的世界里充满了各种各样的周期现象，如日常生活中常见的昼夜更替、四季循环、温度气压等气象因素的反复变化、河水水位的周期性涨落、地面植物的岁月枯荣以及科学技术中诸如太阳活动的十一年周期起伏、地磁场的日变化、重力仪上所反映的固体潮、交流电、电磁波和机械振动等无一不是周期现象。描述周期现象的最简单的周期函数是物理学上所说的谐波函数，它由正弦或余弦函数来表示：tAtycos)((3-1)利用三角公式,上式可以写成：tAtAtAtysinsincoscoscos)(（3-2）由于是常数，令a=Acossin,Ab,则可得：tbtatysincos（3-3）这里22baA,abarctg(3-4)由此可以看出：一个带初相位的余弦函数可以看成一个不带相位的正弦函数与一个不带相位的余弦函数的合成。谐波函数是周期函数中最简单的函数，它描述的也是最简单的周期现象，在实际中所碰到的周期现象往往比它复杂得多。但这些复杂函数均在一定近似程度上可分解为不同频率的正弦函数和余弦函数。下面我们就介绍如何将一种复杂的函数分解为一系列不同频率的正弦函数和余弦函数的方法。3.1.1周期函数的Fourier变换回想一下我们在数学中讲的Fourier级数：如何将一个周期为2l函数分解为Fourier级数呢？我们已经学过的Fourier级数展开式为：10sincos2)(nnnlxnblxnaaxf(3-5)其中dxlxnxflbdxlxnxfladxxflallnllnllsin)(1,cos1,10(3-6)如果f(x)是奇函数,积分上下限相互对称，则这时lxnxfcos亦为奇函数，故an均为零，得到的Fourier级数是正弦级数:lxnbxfnnsin1(3-7)其中，bn的积分可简写为：,.....3,2,1sin20ndxlxnxflbln(3-8)如果f(x)是偶函数,因积分上下限相互对称，并且lxnxfsin为奇函数，故bn均为零,得到的Fourier级数是余弦级数：lxnaaxfnncos210(3-9)其中a0和an可简写为：,......3,2,1cos2,2000ndxlxnxfladxxflalnl(3-10)3．1．2离散Fourier变换现在我们设法把上述公式不加证明地应用于离散Fourier级数中。①我们在数字资料处理中经常遇到的不是一个函数，而是一个离散的数列，比如说，等间隔时间取样的时间序列1210,........,,Nxxxx(这里的数据个数为N,一般取N为偶数,并且若取一特殊的偶数还可以使计算速度加快)。因此用上述公式时需要加以改造。首先，我们得到的数字信号只能在正的时间段取值,在负的时间段不能取值。但由于我们认为所取的时间是无限长的周期序列，周期为2l，因此，我们可以把（-l，l）修改为（0,2l）这样就避免了在负时间段取值。x就对应于时间序列的时间，即tix，i=0,1,2,…,N-1。②由于我们要处理的是离散的数据序列，因此不能再用积分，而应用积分的离散形式----求和来表示，即：Nkl120③我们在（0，2l）里等间隔地取N个取值点，取样时间间隔为tdx，其中2tNl。有了上述改正，我们有12,10,,.........,Nxxxxxf（3-11）NkitNtiklxn22（3-12）所以（3-5）公式的离散形式为：NkibNkiaaxkmkki2sin2cos210（3-13）式中NkixNtNkixtNaxNtxtNaNiiNiikNiiNii2cos22cos21,221101010100（3-14）mkNkixNtNkixtNbNiiNiik,.....,2,1,2sin22sin211010（3-15）现在考虑（3-13）式，若我们共有N个观测数据，即{},.....,,,1210Nxxxx,则上述方程可以列出N个方程式，式（3-14）和式（3-15）为方程组的解。那么，我们考虑一下m最大可以取多大呢？在（3-13）式中kkbaa,,0均是未知数，m越大，未知数的个数越多。（3-13）式所列出的N个方程最多只能确定N个未知数，可见m最大只能取2N。这里需要注意的是，求和号里面的未知数的个数已经有N个了，再加上a0这个未知数，未知数的个数是否多于方程的个数呢？如果我们分析一下求和号里面的kb,当k取2N时,最后一项为:022sin2NiNbN，所以2Nb根本没必要求。但这里需要注意的是，与a0一样，这里的2Na也只取一半。这样，取m最大为2N就正好满足了N个方程确定N个未知数的唯一解的条件。对于N为奇数的情况，m的取值应向下取整得到m=floor(N/2),此时的)2/(Nfloorb不为零，)2/(Nfloora按正常计算。xi的第k项为：NkibNkiakk2sin2cos，我们将其写成以下形式：titNkbtitNkakk2sin2cos（3-16）由上式可以看出，第k项为两个周期函数之和，一个为正弦函数，一个为余弦函数，它们的频率均为：TktNkfk（3-17）这里的ti为时间,圆频率为tNkk2。T为所取数据的总的时间长度，随着k的增大，三角函数的频率逐渐增加，周期逐渐减小。它们的周期为：kTktNfTkk1（3-18）即为整个记录长度的k1。因此，对于一个有限长观测数据序列来说，使用Fourier级数法求得的各种周期或频率的波是有限制的。因为k只能取0~2N的整数。①k=0时，a0为该数据序列的直流分量。②k=1时，波的周期为包含全部观测数据记录的时间T，这是波的最大周期。其频率为Tf1，是用Fourier谐波分析所能检测到的最小频率。通常我们称k=1的波为基波，基波的周期由取样数据的长度来决定。③k取最大值2N时，波的周期为tNTNTNTTN22222。此时波的周期最小，频率最大，为tTfNN2112/2/。在数字信号处理中称这个最大频率为Nyquist频率，这与第二章所讲的取样定理（采样定理）是一致的。Nyquist频率由取样间隔t来决定，为取样间隔2倍的倒数。在中国数字地震台网系统中，每秒取50个样，取样间隔为0.02秒，因此能够得到地震波的最大频率为25赫兹。如果想得到高频率的地震波，必须提高取样频率，即缩短取样间隔。④我们称利用（3-13~15）式分离的各种频率（周期）的波，即（3-16）式为k次谐波。k次谐波的频率可由(3-17)式给出，周期可由（3-18）式给出。若我们将（3-13）式中的求和号里的各次谐波写成（3-2）式的形式，则有：2102cos2NkkkiNkicax（3-19）式中22kkkbac，kkkabarctg（3-20）ck表示了k次谐波的振幅大小。在分解地震波时，哪些分量的振幅大，哪些分量的振幅小，对地震波的性质有十分重要的影响。在研究地震波对工程结构的影响时，因为不同的工程构件有不同的固有频率，地震波中与工程构件的固有频率相等或相近的谐波振幅很大时，会使该构件发生共振而造成破坏。我们通常用频率作为横坐标，用振幅作为纵坐标绘出图形来研究频率和振幅的关系，这种图形就是一种重要的谱，通常称为Fourier振幅谱，简称振幅谱(振幅spectra)。k为k次谐波的初相，将（3-20）式中波的各分量的初相按频率画出，与振幅谱对应，我们称之为相位谱（Phasespectra）。如果将k转换为我们平常看到的角度，只需将k乘以0180即可。(3-19)式和（3-20）式之所以重要是因为它将一个数字序列转化为一系列不同频率的波的合成。这样作的好处是通过这种分解，我们可以看出所给波列中含有哪种频率成分。如果我们想滤掉某种频率的波，直接去除该频率的谐波即可。下面我们根据上式给出一个合成波分解为各种不同频率波的例子。3．1．3信号的Fourier分解与合成举例【例3-1】将振幅为1的1Hz正弦波和振幅为0.5的5Hz正弦波相加后进行Fourier分析，研究能否从中分析出含有这两种频率的信号。%Samp3_1clearall%将工作空间中的所有变量清除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两个为零的序列fork=0:m-1forii=0:N-1a(k+1)=a(k+1)+2/N*x(ii+1)*cos(2*pi*k*ii/N);%NkixNaNiik2cos210b(k+1)=b(k+1)+2/N*x(ii+1)*sin(2*pi*k*ii/N);%NkixNbNiik2sin210%MATLAB中的数组序号只能从1开始。endc(k+1)=sqrt(a(k+1).^2+b(k+1).^2);%22kkkbac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('振幅')0123456-1.5-1-0.500.511.5原始信号时间/s051015202500.20.40.60.81Fourier变换频率/Hz振幅图3-1将1Hz和5Hz的正弦振动合成后进行Fourier分析的结果由图3-1可以清楚地看出，从初始信号中明确地识别出了1Hz和5Hz的波。这里1Hz和5Hz的振幅与原来信号振幅并不完全一致，这是由于数据采样点较少造成的。自己作实验可以看出，数据采样点越多，Fourier分析得到的结果与原始振幅越接近。如果将m修改为N,其频谱只是增加了镜像对称的部分，如果将m改为2*N，则频谱增加一个周期（见图3-2），改为3*N，频谱就增加了三个周期，可以明确地看出频谱的周期性。自己修改上面的程序可以验证这一点(参看程序Samp3_1_1)。由此可以看出，离散有限信号的频谱为周期谱。0510152025303540455000.20.40.60.81Fourier变换频率/Hz振幅010203040506070809010000.20.40.60.81Fourier变换频率/Hz振幅图3-2取大于Nyquist频率时振幅谱的例子图示。上图，Nyquist频率的2倍；下图，Nyquist的4倍上面的例子表明，k值取不小于N/2时，可以完全分辨出小于采样频率一半（Nyquist频率）的频率，反过来，根据分解成的ak、bk值能否准确地恢复原来的信号呢？下面我们采用（3-10）式或（3-15）式求得上例合成信号，在前面程序运行后，运行下