搜索
理解离散傅立叶变换（一）------傅立叶变换的由来一、傅立叶变换的提出让我们先看看为什么会有傅立叶变换？傅立叶是一位法国数学家和物理学家的名字，英语原名是JeanBaptisteJosephFourier(1768-1830),Fourier对热传递很感兴趣，于1807年在法国科学学会上发表了一篇论文，运用正弦曲线来描述温度分布，论文里有个在当时具有争议性的决断：任何连续周期信号可以由一组适当的正弦曲线组合而成。当时审查这个论文的人，其中有两位是历史上著名的数学家拉格朗日(JosephLouisLagrange,1736-1813)和拉普拉斯(PierreSimondeLaplace,1749-1827)，当拉普拉斯和其它审查者投票通过并要发表这个论文时，拉格朗日坚决反对，在近50年的时间里，拉格朗日坚持认为傅立叶的方法无法表示带有棱角的信号，如在方波中出现非连续变化斜率。法国科学学会屈服于拉格朗日的威望，拒绝了傅立叶的工作，幸运的是，傅立叶还有其它事情可忙，他参加了政治运动，随拿破仑远征埃及，法国大革命后因会被推上断头台而一直在逃避。直到拉格朗日死后15年这个论文才被发表出来。拉格朗日是对的：正弦曲线无法组合成一个带有棱角的信号。但是，我们可以用正弦曲线来非常逼近地表示它，逼近到两种表示方法不存在能量差别，基于此，傅立叶是对的。为什么我们要用正弦曲线来代替原来的曲线呢？如我们也还可以用方波或三角波来代替呀，分解信号的方法是无穷的，但分解信号的目的是为了更加简单地处理原来的信号。用正余弦来表示原信号会更加简单，因为正余弦拥有原信号所不具有的性质：正弦曲线保真度。一个正弦曲线信号输入后，输出的仍是正弦曲线，只有幅度和相位可能发生变化，但是频率和波的形状仍是一样的。且只有正弦曲线才拥有这样的性质，正因如此我们才不用方波或三角波来表示。二、傅立叶变换分类根据原信号的不同类型，我们可以把傅立叶变换分为四种类别：1非周期性连续信号傅立叶变换（FourierTransform）2周期性连续信号傅立叶级数(FourierSeries)3非周期性离散信号离散时域傅立叶变换（DiscreteTimeFourierTransform）4周期性离散信号离散傅立叶变换(DiscreteFourierTransform)下图是四种原信号图例：这四种傅立叶变换都是针对正无穷大和负无穷大的信号，即信号的的长度是无穷大的，我们知道这对于计算机处理来说是不可能的，那么有没有针对长度有限的傅立叶变换呢？没有。因为正余弦波被定义成从负无穷小到正无穷大，我们无法把一个长度无限的信号组合成长度有限的信号。面对这种困难，方法是把长度有限的信号表示成长度无限的信号，可以把信号无限地从左右进行延伸，延伸的部分用零来表示，这样，这个信号就可以被看成是非周期性离散信号，我们就可以用到离散时域傅立叶变换的方法。还有，也可以把信号用复制的方法进行延伸，这样信号就变成了周期性离解信号，这时我们就可以用离散傅立叶变换方法进行变换。这里我们要学的是离散信号，对于连续信号我们不作讨论，因为计算机只能处理离散的数值信号，我们的最终目的是运用计算机来处理信号的。但是对于非周期性的信号，我们需要用无穷多不同频率的正弦曲线来表示，这对于计算机来说是不可能实现的。所以对于离散信号的变换只有离散傅立叶变换（DFT）才能被适用，对于计算机来说只有离散的和有限长度的数据才能被处理，对于其它的变换类型只有在数学演算中才能用到，在计算机面前我们只能用DFT方法，后面我们要理解的也正是DFT方法。这里要理解的是我们使用周期性的信号目的是为了能够用数学方法来解决问题，至于考虑周期性信号是从哪里得到或怎样得到是无意义的。每种傅立叶变换都分成实数和复数两种方法，对于实数方法是最好理解的，但是复数方法就相对复杂许多了，需要懂得有关复数的理论知识，不过，如果理解了实数离散傅立叶变换(realDFT)，再去理解复数傅立叶就更容易了，所以我们先把复数的傅立叶放到一边去，先来理解实数傅立叶变换，在后面我们会先讲讲关于复数的基本理论，然后在理解了实数傅立叶变换的基础上再来理解复数傅立叶变换。还有，这里我们所要说的变换(transform)虽然是数学意义上的变换，但跟函数变换是不同的，函数变换是符合一一映射准则的，对于离散数字信号处理（DSP），有许多的变换：傅立叶变换、拉普拉斯变换、Z变换、希尔伯特变换、离散余弦变换等，这些都扩展了函数变换的定义，允许输入和输出有多种的值，简单地说变换就是把一堆的数据变成另一堆的数据的方法。三、一个关于实数离散傅立叶变换(RealDFT)的例子先来看一个变换实例，下图是一个原始信号图像：这个信号的长度是16，于是可以把这个信号分解9个余弦波和9个正弦波（一个长度为N的信号可以分解成N/2+1个正余弦信号，这是为什么呢？结合下面的18个正余弦图,我想从计算机处理精度上就不难理解，一个长度为N的信号，最多只能有N/2+1个不同频率，再多的频率就超过了计算机所能所处理的精度范围），如下图：9个余弦信号：9个正弦信号：把以上所有信号相加即可得到原始信号，至于是怎么分别变换出9种不同频率信号的，我们先不急，先看看对于以上的变换结果，在程序中又是该怎么表示的，我们可以看看下面这个示例图：上图中左边表示时域中的信号，右边是频域信号表示方法，从左向右表示正向转换(ForwardDFT)，从右向左表示逆向转换(InverseDFT)，用小写x[]表示信号在每个时间点上的幅度值数组,用大写X[]表示每种频率的幅度值数组,因为有N/2+1种频率，所以该数组长度为N/2+1，X[]数组又分两种，一种是表示余弦波的不同频率幅度值：ReX[]，另一种是表示正弦波的不同频率幅度值：ImX[]，Re是实数(Real)的意思，Im是虚数(Imagine)的意思，采用复数的表示方法把正余弦波组合起来进行表示，但这里我们不考虑复数的其它作用，只记住是一种组合方法而已，目的是为了便于表达（在后面我们会知道，复数形式的傅立叶变换长度是N，而不是N/2+1）。下一节我们将来看一下实数傅立叶变换的具体方法。理解离散傅立叶变换（二）------实数形式离散傅立叶变换（RealDFT）上一节我们看到了一个实数形式离散傅立叶变换的例子，通过这个例子能够让我们先对傅立叶变换有一个较为形象的感性认识，现在就让我们来看看实数形式离散傅立叶变换的正向和逆向是怎么进行变换的。在此，我们先来看一下频率的多种表示方法。一、频域中关于频率的四种表示方法1、序号表示方法，根据时域中信号的样本数取0~N/2，用这种方法在程序中使用起来可以更直接地取得每种频率的幅度值，因为频率值跟数组的序号是一一对应的:X[k]，取值范围是0~N/2；2、分数表示方法，根据时域中信号的样本数的比例值取0~0.5:X[ƒ]，ƒ=k/N，取值范围是0~N/2；3、用弧度值来表示，把ƒ乘以一个2π得到一个弧度值，这种表示方法叫做自然频率(naturalfrequency)：X[ω]，ω=2πƒ=2πk/N，取值范围是0~π；4、以赫兹(Hz)为单位来表示，这个一般是应用于一些特殊应用，如取样率为10kHz表示每秒有10,000个样本数：取值范围是0到取样率的一半。二、DFT基本函数ck[i]=cos(2πki/N)sk[i]=sin(2πki/N)其中k表示每个正余弦波的频率，如为2表示在0到N长度中存在两个完整的周期，10即有10个周期，如下图：上图中至于每个波的振幅(amplitude)值(ReX[k],ImX[k])是怎么算出来的,这个是DFT的核心，也是最难理解的部分，我们先来看看如何把分解出来的正余弦波合成原始信号(InverseDFT)。三、合成运算方法(RealInverseDFT)DFT合成等式：如果有学过傅立叶级数，对这个等式就会有似曾相识的感觉，不错！这个等式跟傅立叶级数是非常相似的：当然，差别是肯定是存在的，因为这两个等式是在两个不同条件下运用的，至于怎么证明DFT合成公式，这个我想需要非常强的高等数学理论知识了，这是研究数学的人的工作，对于普通应用者就不需要如此的追根究底了，但是傅立叶级数是好理解的，我们起码可以从傅立叶级数公式中看出DFT合成公式的合理性。DFT合成等式中的ImX[k]和ReX[k]跟ImX[k]和ReX[k]是不一样的，下面是转换方法:但k等于0和N/2时,实数部分的计算要用下面的等式:上面四个式中的N是时域中点的总数，k是从0到N/2的序号。为什么要这样进行转换呢？这个可以从频谱密度(spectraldensity)得到理解，如下图就是个频谱图：这是一个频谱图，横坐标表示频率大小，纵坐标表示振幅大小，原始信号长度为N（这里是32），经DFT转换后得到的17个频率的频谱，频谱密度表示每单位带宽中为多大的振幅，那么带宽是怎么计算出来的呢？看上图，除了头尾两个，其余点的所占的宽度是2/N，这个宽度便是每个点的带宽，头尾两个点的带宽是1/N,而ImX[k]和ReX[k]表示的是频谱密度，即每一个单位带宽的振幅大小，但Im[k]和Re[k]表示2/N（或1/N）带宽的振幅大小，所以Im[k]和Re[k]分别应当是ImX[k]和ReX[k]的2/N（或1/N）。频谱密度就象物理中物质密度，原始信号中的每一个点就象是一个混合物，这个混合物是由不同密度的物质组成的，混合物中含有的每种物质的质量是一样的，除了最大和最小两个密度的物质外，这样我们只要把每种物质的密度加起来就可以得到该混合物的密度了，又该混合物的质量是单位质量，所以得到的密度值跟该混合物的质量值是一样的。至于为什么虚数部分是负数，这是为了跟复数DFT保持一致，这个我们将在后面会知道这是数学计算上的需要（ImX[k]计算时加上了一个负号，Im[k]再加上负号，结果便是正的，等于没有变化）。如果已经得到了DFT结果，这时要进行逆转换，即合成原始信号，则可按如下步骤进行转换：1、先根据上面四个式子计算得出Im[k]和Re[k]的值；2、再根据DFT合成等式得到原始信号数据。下面是用BASIC语言来实现的转换源代码：100‘DFT逆转换方法110‘/XX[]数组存储计算结果（时域中的原始信号）120‘/REX[]数组存储频域中的实数分量，IMX[]为虚分量130‘140DIMXX[511]150DIMREX[256]160DIMIMX[256]170‘180PI=3.14159265190N%=512200‘210GOSUBXXXX‘转到子函数去获取REX[]和IMX[]数据220‘230‘240‘250FORK%=0TO256260REX[K%]=REX[K%]/(N%/2)270IMX[K%]=-IMX[K%]/(N%/2)280NEXTk%290‘300REX[0]=REX[0]/N310REX[256]=REX[256]/N320‘330‘初始化XX[]数组340FORI%=0TO511350XX[I%]=0360NEXTI%370‘380‘390‘400‘410‘420FORK%=0TO256430FORI%=0TO511440‘450XX[I%]=XX[I%]+REX[K%]*COS(2*PI*K%*I%/N%)460XX[I%]=XX[I%]+IMX[K%]*SIN(2*PI*K%*I%/N%)470‘480NEXTI%490NEXTK%500‘510END上面代码中420至490换成如下形式也许更好理解，但结果都是一样的：420FORI%=0TO511430FORK%=0TO256440‘450XX[I%]=XX[I%]+REX[K%]*COS(2*PI*K%*I%/N%)460XX[I%]=XX[I%]+IMX[K%]*SIN(2*PI*K%*I%/N%)470‘480NEXTI%490NEXTK%四、分解运算方法（DFT）有三种完全不同的方法进行DFT：一种方法是通过联立方程进行求解,从代数的角度看，要从N个已知值求N个未知值，需要N个联立方程，且N个联立方程必须是线