从头到尾彻底理解傅里叶变换算法

r囧r小猫 2022-07-15 13:55 240阅读 0赞

**http://www.xuebuyuan.com/2052774.html**

**经典[算法][Link 1]研究系列：十、从头到尾彻底理解傅里叶变换算法、上**

作者：July、dznlong 二零一一年二月二十日

博主说明：I、本文中阐述离散傅里叶变换方法，是根据此书：The Scientist and Engineer's Guide to Digital Signal Processing，By Steven W. Smith, Ph.D.而翻译而成的，此书地址：[http://www.dspguide.com/pdfbook.htm][http_www.dspguide.com_pdfbook.htm]。II、同时，有相当一部分内容编辑整理自dznlong的博客，也贴出其博客地址，向原创的作者表示致敬：[http://blog.csdn.net/dznlong][http_blog.csdn.net_dznlong] 。这年头，真正静下心写来原创[文章][Link 2]的人，很少了。  
\------------------------------------  
**从头到尾彻底理解傅里叶变换算法、上  
前言  
第一部分、 DFT  
第一章、傅立叶变换的由来  
第二章、实数形式离散傅立叶变换（Real DFT）**

[**从头到尾彻底理解傅里叶变换算法、下**][Link 3]

第三章、复数  
第四章、复数形式离散傅立叶变换

**前言：**  
“关于傅立叶变换，无论是书本还是在网上可以很容易找到关于傅立叶变换的描述，但是大都是些故弄玄虚的文章，太过抽象，尽是一些让人看了就望而生畏的公式的罗列，让人很难能够从感性上得到理解”---dznlong，

那么，到底什么是傅里叶变换算法列?傅里叶变换所涉及到的公式具体有多复杂列?  
**傅里叶变换**（Fourier transform）是一种线性的积分变换。因其基本思想首先由法国学者傅里叶系统地提出，所以以其名字来命名以示纪念。

哦，傅里叶变换原来就是一种变换而已，只是这种变换是从时间转换为频率的变化。这下，你就知道了，傅里叶就是一种变换，一种什么变换列?就是一种从时间到频率的变化或其相互转化。

ok，咱们再来总体了解下傅里叶变换，让各位对其有个总体大概的印象，也顺便看看傅里叶变换所涉及到的公式，究竟有多复杂：  
以下就是傅里叶变换的4种变体（摘自，维基百科）  
**连续傅里叶变换**  
 一般情况下，若“傅里叶变换”一词不加任何限定语，则指的是“连续傅里叶变换”。连续傅里叶变换将平方可积的函数f（t）表示成复指数函数的积分或级数形式。

![8394323_1298253739kqo6.jpg][]

这是将频率域的函数F(ω)表示为时间域的函数f（t）的积分形式。

连续傅里叶变换的逆变换 (inverse Fourier transform)为：

![8394323_1298253739WHO6.jpg][]

即将时间域的函数f（t）表示为频率域的函数F(ω)的积分。

一般可称函数**f（t）**为原函数，而称函数**F(ω)**为傅里叶变换的像函数，原函数和像函数构成一个傅里叶变换对（transform pair）。

除此之外，还有其它型式的变换对，以下两种型式亦常被使用。在通信或是信号处理方面，常以![8394323_1298255200A7rr.jpg][]来代换，而形成新的变换对:

![8394323_1298253736OmpV.jpg][]

或者是因系数重分配而得到新的变换对：

![8394323_1298253736ntYg.jpg][]

一种对连续傅里叶变换的推广称为分数傅里叶变换（Fractional Fourier Transform）。分数傅里叶变换(fractional Fourier transform,FRFT)指的就是傅里叶变换(Fourier transform,FT)的广义化。  
分数傅里叶变换的物理意义即做傅里叶变换 a 次，其中 a 不一定要为整数；而做了分数傅里叶变换之后，信号或输入函数便会出现在介于时域(time domain)与频域(frequency domain)之间的分数域(fractional domain)。

当f（t）为偶函数（或奇函数）时，其正弦（或余弦）分量将消亡，而可以称这时的变换为余弦变换（cosine transform）或正弦变换（sine transform）.

另一个值得注意的性质是，当f（t）为纯实函数时，F(−ω) = F\*(ω)成立.

**傅里叶级数**  
 连续形式的傅里叶变换其实是傅里叶级数 (Fourier series)的推广，因为积分其实是一种极限形式的求和算子而已。对于周期函数，其傅里叶级数是存在的：

![8394323_12982537341B19.jpg][]

其中Fn为复幅度。对于实值函数，函数的傅里叶级数可以写成：

![8394323_1298253734CxmE.jpg][]  
其中an和bn是实频率分量的幅度。

**离散时域傅里叶变换**  
 离散傅里叶变换是离散时间傅里叶变换（DTFT）的特例（有时作为后者的近似）。DTFT在时域上离散，在频域上则是周期的。DTFT可以被看作是傅里叶级数的逆变换。

**离散傅里叶变换**  
 离散傅里叶变换（DFT），是连续傅里叶变换在时域和频域上都离散的形式，将时域信号的采样变换为在离散时间傅里叶变换（DTFT）频域的采样。在形式上，变换两端（时域和频域上）的序列是有限长的，而实际上这两组序列都应当被认为是离散周期信号的主值序列。即使对有限长的离散信号作DFT，也应当将其看作经过周期延拓成为周期信号再作变换。在实际应用中通常采用快速傅里叶变换以高效计算DFT。

为了在科学计算和数字信号处理等领域使用计算机进行傅里叶变换，必须将函数xn定义在离散点而非连续域内，且须满足有限性或周期性条件。这种情况下，使用离散傅里叶变换（DFT），将函数xn表示为下面的求和形式：

![8394323_1298253734oZvG.jpg][]

其中Xk是傅里叶幅度。直接使用这个公式计算的计算复杂度为O（n\*n），而快速傅里叶变换（FFT）可以将复杂度改进为O（n\*lgn）。（后面会具体阐述FFT是如何将复杂度降为O（n\*lgn）的。）计算复杂度的降低以及数字电路计算能力的发展使得DFT成为在信号处理领域十分实用且重要的方法。

下面，比较下上述傅立叶变换的4种变体，

![8394323_1298266276GaFF.jpg][]

如上，容易发现：函数在时（频）域的离散对应于其像函数在频（时）域的周期性。反之连续则意味着在对应域的信号的非周期性。也就是说，时间上的离散性对应着频率上的周期性。同时，注意，离散时间傅里叶变换，时间离散，频率不离散，它在频域依然是连续的。  
如果，读到此，你不甚明白，大没关系，不必纠结于以上4种变体，继续往下看，你自会豁然开朗。（有什么问题，也恳请提出，或者批评指正）

**ok， 本文，接下来，由傅里叶变换入手，后重点阐述离散傅里叶变换、快速傅里叶算法，到最后彻底实现FFT算法，全篇力求通俗易懂、阅读顺畅，教你从头到尾彻底理解傅里叶变换算法。**由于傅里叶变换，也称傅立叶变换，下文所称为傅立叶变换，同一个变换，不同叫法，读者不必感到奇怪。

**第一部分、DFT  
第一章、傅立叶变换的由来**  
 要理解傅立叶变换，先得知道傅立叶变换是怎么变换的，当然，也需要一定的高等数学基础，最基本的是级数变换，其中傅立叶级数变换是傅立叶变换的基础公式。  
  
**一、傅立叶变换的提出**

傅立叶是一位法国数学家和物理学家，原名是Jean Baptiste Joseph Fourier(1768-1830), Fourier于1807年在法国科学学会上发表了一篇论文，论文里描述运用正弦曲线来描述温度分布，论文里有个在当时具有争议性的决断：任何连续周期信号都可以由一组适当的正弦曲线组合而成。

当时审查这个论文拉格朗日坚决反对此论文的发表，而后在近50年的时间里，拉格朗日坚持认为傅立叶的方法无法表示带有棱角的信号，如在方波中出现非连续变化斜率。直到拉格朗日死后15年这个论文才被发表出来。  
谁是对的呢？拉格朗日是对的：正弦曲线无法组合成一个带有棱角的信号。但是，我们可以用正弦曲线来非常逼近地表示它，逼近到两种表示方法不存在能量差别，基于此，傅立叶是对的。

为什么我们要用正弦曲线来代替原来的曲线呢？如我们也还可以用方波或三角波来代替呀，分解信号的方法是无穷多的，但分解信号的目的是为了更加简单地处理原来的信号。  
用正余弦来表示原信号会更加简单，因为正余弦拥有原信号所不具有的性质：正弦曲线保真度。一个正余弦曲线信号输入后，输出的仍是正余弦曲线，只有幅度和相位可能发生变化，但是频率和波的形状仍是一样的。且只有正余弦曲线才拥有这样的性质，正因如此我们才不用方波或三角波来表示。

**二、傅立叶变换分类**  
 根据原信号的不同类型，我们可以把傅立叶变换分为四种类别：  
1、非周期性连续信号 傅立叶变换（Fourier Transform）  
2、周期性连续信号 傅立叶级数(Fourier Series)  
3、非周期性离散信号 离散**时域**傅立叶变换（Discrete Time Fourier Transform）  
4、周期性离散信号 离散傅立叶变换(Discrete Fourier Transform)  
下图是四种原信号图例（从上到下，依次是FT，FS，DTFT，DFT）：

![8394323_1298213518qfT8.jpg][]

这四种傅立叶变换都是针对正无穷大和负无穷大的信号，即信号的的长度是无穷大的，我们知道这对于计算机处理来说是不可能的，**那么有没有针对长度有限的傅立叶变换呢？**没有。因为正余弦波被定义成从负无穷小到正无穷大，我们无法把一个长度无限的信号组合成长度有限的信号。  
面对这种困难，方法是：把长度有限的信号表示成长度无限的信号。如，可以把信号无限地从左右进行延伸，延伸的部分用零来表示，这样，这个信号就可以被看成是非周期性离散信号**，**我们可以用到离散时域傅立叶变换（DTFT）的方法。也可以把信号用复制的方法进行延伸，这样信号就变成了周期性离散信号，这时我们就可以用**离散傅立叶变换方法（DFT）**进行变换。本章我们要讲的是**离散信号**，对于连续信号我们不作讨论，因为计算机只能处理离散的数值信号，我们的最终目的是运用计算机来处理信号的。  
  
但是对于非周期性的信号，我们需要用无穷多不同频率的正弦曲线来表示，这对于计算机来说是不可能实现的。所以对于离散信号的变换只有**离散傅立叶变换（DFT）**才能被适用，对于计算机来说只有离散的和有限长度的数据才能被处理，对于其它的变换类型只有在数学演算中才能用到，在计算机面前我们只能用DFT方法，后面我们要理解的也正是DFT方法。  
这里要理解的是我们使用周期性的信号目的是为了能够用数学方法来解决问题，至于考虑周期性信号是从哪里得到或怎样得到是无意义的。  
  
每种傅立叶变换都分成实数和复数两种方法，对于实数方法是最好理解的，但是复数方法就相对复杂许多了，需要懂得有关复数的理论知识，不过，如果理解了实数离散傅立叶变换(real DFT)，再去理解复数傅立叶变换就更容易了，所以我们先把复数的傅立叶变换放到一边去，先来理解实数傅立叶变换，在后面我们会先讲讲关于复数的基本理论，然后在理解了实数傅立叶变换的基础上再来理解复数傅立叶变换。  
  
还有，这里我们所要说的变换(transform)虽然是数学意义上的变换，但跟函数变换是不同的，函数变换是符合一一映射准则的，对于离散数字信号处理（DSP），有许多的变换：傅立叶变换、拉普拉斯变换、Z变换、希尔伯特变换、离散余弦变换等，这些都扩展了函数变换的定义，允许输入和输出有多种的值，简单地说变换就是把一堆的数据变成另一堆的数据的方法。  
  
**三、一个关于实数离散傅立叶变换(Real DFT)的例子**

先来看一个变换实例，下图是一个原始信号图像：

![8394323_1298213518zwwp.jpg][]  
这个信号的长度是16，于是可以把这个信号分解9个余弦波和9个正弦波（一个长度为N的信号可以分解成N/2+1个正余弦信号，这是为什么呢？结合下面的18个正余弦图,我想从计算机处理精度上就不难理解，一个长度为N的信号，最多只能有N/2+1个不同频率，再多的频率就超过了计算机所能所处理的精度范围），如下图：

9个余弦信号：

![8394323_12982135185EPi.jpg][]

9个正弦信号：

![8394323_12982135172NjF.jpg][]

把以上所有信号相加即可得到原始信号，至于是怎么分别变换出9种不同频率信号的，我们先不急，先看看对于以上的变换结果，在[程序][Link 4]中又是该怎么表示的，我们可以看看下面这个示例图：

![8394323_1298213517sJBS.jpg][]  
  
上图中左边表示时域中的信号，右边是频域信号表示方法，  
从左向右，**-->**，表示**正向转换**(Forward DFT)，从右向左，**<--**，表示**逆向转换**(Inverse DFT)，  
用小写x\[\]表示信号在每个时间点上的幅度值数组, 用大写X\[\]表示每种频率的副度值数组（即时间x-->频率**X**）,  
因为有N/2+1种频率，所以该数组长度为N/2+1，  
X\[\]数组又分两种，一种是表示余弦波的不同频率幅度值：Re X\[\]，  
另一种是表示正弦波的不同频率幅度值：Im X\[\]，  
Re是实数(Real)的意思，Im是虚数(Imagine)的意思，采用复数的表示方法把正余弦波组合起来进行表示，但这里我们不考虑复数的其它作用，只记住是一种组合方法而已，目的是为了便于表达（在后面我们会知道，复数形式的傅立叶变换长度是N，而不是N/2+1）。如此，再回过头去，看上面的正余弦各9种频率的变化，相信，问题不大了。

**第二章、实数形式离散傅立叶变换（Real DFT）**  
 上一章，我们看到了一个实数形式离散傅立叶变换的例子，通过这个例子能够让我们先对傅立叶变换有一个较为形象的感性认识，现在就让我们来看看实数形式离散傅立叶变换的正向和逆向是怎么进行变换的。在此，我们先来看一下频率的多种表示方法。  
  
**一、 频域中关于频率的四种表示方法**  
  
1、序号表示方法，根据时域中信号的样本数取0 ~ N/2，用这种方法在程序中使用起来可以更直接地取得每种频率的幅度值，因为频率值跟数组的序号是一一对应的: X\[k\]，取值范围是0 ~ N/2；  
2、分数表示方法，根据时域中信号的样本数的比例值取0 ~ 0.5: X\[ƒ\]，ƒ = k/N，取值范围是0 ~ 1/2；  
3、用弧度值来表示，把ƒ乘以一个2π得到一个弧度值，这种表示方法叫做自然频率(natural frequency)：X\[ω\]，ω = 2πƒ = 2πk/N，取值范围是0 ~ π；  
4、以赫兹(Hz)为单位来表示，这个一般是应用于一些特殊应用，如取样率为10 kHz表示每秒有10,000个样本数：取值范围是0到取样率的一半。  
  
**二、 DFT基本函数**  
  
ck\[i\] = cos(2πki/N)  
sk\[i\] = sin(2πki/N)  
其中k表示每个正余弦波的频率，如为2表示在0到N长度中存在两个完整的周期，10即有10个周期，如下图：

![8394323_1298213516D7V1.jpg][]

上图中至于每个波的振幅(amplitude)值(Re X\[k\],Im X\[k\])是怎么算出来的,这个是DFT的核心，也是最难理解的部分，我们先来看看如何把分解出来的正余弦波合成原始信号(Inverse DFT)。  
  
**三、 合成运算方法(Real Inverse DFT)**  
  
DFT合成等式（合成原始**时间**信号，频率-->时间，逆向变换）：

![8394323_1298213516k4UT.jpg][]

如果有学过傅立叶级数，对这个等式就会有似曾相识的感觉，不错！这个等式跟傅立叶级数是非常相似的：

![8394323_1298213516KUWW.jpg][]

当然，差别是肯定是存在的，因为这两个等式是在两个不同条件下运用的，至于怎么证明DFT合成公式，这个我想需要非常强的高等数学理论知识了，这是研究数学的人的工作，对于普通应用者就不需要如此的追根究底了，但是傅立叶级数是好理解的，我们起码可以从傅立叶级数公式中看出DFT合成公式的合理性。  
\_ \_  
DFT合成等式中的Im X\[k\]和Re X\[k\]跟之前提到的Im X\[k\]和Re X\[k\]是不一样的，下面是转换方法（关于此公式的解释，见下文）:

![8394323_1298213515r5dm.jpg][]

但k等于0和N/2时,实数部分的计算要用下面的等式:

![8394323_12982135155Kvg.jpg][]

上面四个式中的N是时域中点的总数，k是从0到N/2的序号。  
为什么要这样进行转换呢？这个可以从频谱密度(spectral density)得到理解，如下图就是个频谱图：

![8394323_12982135156Bee.jpg][]

这是一个频谱图，横坐标表示频率大小，纵坐标表示振幅大小，原始信号长度为N（这里是**32**），经DFT转换后得到的**17**个频率的频谱，频谱密度表示每单位带宽中为多大的振幅，那么带宽是怎么计算出来的呢？看上图，除了头尾两个，其余点的所占的宽度是2/N，这个宽度便是每个点的带宽，头尾两个点的带宽是1/N,而Im X\[k\]和Re X\[k\]表示的是频谱密度，即每一个单位带宽的振幅大小，但![8394323_1298270500e82R.jpg][]表示2/N（或1/N）带宽的振幅大小，**所以![8394323_1298270500e82R.jpg][]分别应当是Im X\[k\]和Re X\[k\]的2/N（或1/N）。**  
  
频谱密度就象物理中物质密度，原始信号中的每一个点就象是一个混合物，这个混合物是由不同密度的物质组成的，混合物中含有的每种物质的质量是一样的，除了最大和最小两个密度的物质外，这样我们只要把每种物质的密度加起来就可以得到该混合物的密度了，又该混合物的质量是单位质量，所以得到的密度值跟该混合物的质量值是一样的。  
  
至于为什么虚数部分是负数，这是为了跟复数DFT保持一致，这个我们将在后面会知道这是数学计算上的需要（Im X\[k\]在计算时就已经加上了一个负号（稍后，由下文，便可知），![8394323_1298272202VxKp.jpg][]再加上负号，结果便是正的，等于没有变化）。  
  
如果已经得到了DFT结果，这时要进行**逆转换，即合成原始信号**，则可按如下步骤进行转换：  
1、先根据上面四个式子计算得出![8394323_1298270500e82R.jpg][]的值；  
2、再根据DFT合成等式得到原始信号数据。  
下面是用BASIC语言来实现的转换源[代码][Link 4]：  
100 ‘DFT逆转换方法  
110 ‘/XX\[\]数组存储计算结果（时域中的原始信号）  
120 ‘/REX\[\]数组存储频域中的实数分量，IMX\[\]为虚分量  
130 ‘  
140 DIM XX\[511\]  
150 DIM REX\[256\]  
160 DIM IMX\[256\]  
170 ‘  
180 PI = 3.14159265  
190 N% = 512  
200 ‘  
210 GOSUB XXXX ‘转到子函数去获取REX\[\]和IMX\[\]数据  
220 ‘  
230 ‘  
240 ‘  
250 FOR K% = 0 TO 256  
260 REX\[K%\] = REX\[K%\] / (N%/2)  
270 IMX\[K%\] = -IMX\[K%\] / (N%/2)  
280 NEXT k%  
290 ‘  
300 REX\[0\] = REX\[0\] / N  
310 REX\[256\] = REX\[256\] / N  
320 ‘  
330 ‘ 初始化XX\[\]数组  
340 FOR I% = 0 TO 511  
350 XX\[I%\] = 0  
360 NEXT I%  
370 ‘  
380 ‘  
390 ‘  
400 ‘  
410 ‘  
420 FOR K% =0 TO 256  
430 FOR I%=0 TO 511  
440 ‘  
450 XX\[I%\] = XX\[I%\] + REX\[K%\] \* COS(2 \* PI \* K% \* I% / N%)  
460 XX\[I%\] = XX\[I%\] + IMX\[K%\] \* SIN(2 \* PI \* K% \* I% / N%)  
470 ‘  
480 NEXT I%  
490 NEXT K%  
500 ‘  
510 END  
  
上面代码中420至490换成如下形式也许更好理解，但结果都是一样的：  
420 FOR I% =0 TO 511  
430 FOR K%=0 TO 256  
440 ‘  
450 XX\[I%\] = XX\[I%\] + REX\[K%\] \* COS(2 \* PI \* K% \* I% / N%)  
460 XX\[I%\] = XX\[I%\] + IMX\[K%\] \* SIN(2 \* PI \* K% \* I% / N%)  
470 ‘  
480 NEXT I%  
490 NEXT K%  
  
**四、 分解运算方法（DFT）**  
  
有三种完全不同的方法进行DFT：一种方法是通过联[立方][Link 5]程进行求解, 从代数的角度看，要从N个已知值求N个未知值，需要N个联立方程，且N个联立方程必须是线性独立的，但这是这种方法计算量非常的大且极其复杂，所以很少被采用；第二种方法是利用信号的相关性（correlation）进行计算，这个是我们后面将要介绍的方法；第三种方法是快速傅立叶变换（FFT），这是一个非常具有创造性和革命性的的方法，因为它大大提高了运算速度，使得傅立叶变换能够在计算机中被广泛应用，但这种算法是根据复数形式的傅立叶变换来实现的，它把N个点的信号分解成长度为N的频域，这个跟我们现在所进行的实域DFT变换不一样，而且这种方法也较难理解，这里我们先不去理解，等先理解了复数DFT后，再来看一下FFT。有一点很重要，那就是这三种方法所得的变换结果是一样的，经过实践证明，当频域长度为32时，利用相关性方法进行计算效率最好，否则FFT算法效率较高。现在就让我们来看一下相关性算法。  
  
利用**第一种方法、信号的相关性**(correlation)可以从噪声背景中检测出已知的信号，我们也可以利用这个方法检测信号波中是否含有某个频率的信号波：把一个待检测信号波乘以另一个信号波，得到一个新的信号波，再把这个新的信号波所有的点进行相加，从相加的结果就可以判断出这两个信号的相似程度。如下图：

![8394323_1298213514wILB.jpg][]

上面a和 b两个图是待检测信号波，图a很明显可以看出是个3个周期的正弦信号波，图b的信号波则看不出是否含有正弦或余弦信号，图c和d都是个3个周期的正弦信号波，图e和f分别是a、b两图跟c、d两图相乘后的结果，图e所有点的平均值是0.5，说明信号a含有振幅为1的正弦信号c，但图f所有点的平均值是0，则说明信号b不含有信号d。这个就是通过信号相关性来检测是否含有某个信号的方法。  
  
 **第二种方法：**相应地，我也可以通过把输入信号和每一种频率的正余弦信号进行相乘（**关联操作**），从而得到原始信号与每种频率的关联程度（即总和大小），这个结果便是我们所要的傅立叶变换结果，下面两个等式便是我们所要的计算方法：

![8394323_12982135144LLa.jpg][]

第二个式子中加了个负号，是为了保持复数形式的一致，前面我们知道在计算![8394323_1298272202VxKp.jpg][]时又加了个负号，所以这只是个形式的问题，并没有实际意义，你也可以把负号去掉，并在计算![8394323_1298272202VxKp.jpg][]时也不加负号。

这里有一点必须明白一个正交的概念：两个函数相乘，如果结果中的每个点的总和为0，则可认为这两个函数为正交函数。要确保关联性算法是正确的，则必须使得跟原始信号相乘的信号的函数形式是正交的，我们知道所有的正弦或余弦函数是正交的，这一点我们可以通过简单的高数知识就可以证明它，所以我们可以通过关联的方法把原始信号分离出正余弦信号。当然，其它的正交函数也是存在的，如：方波、三角波等形式的脉冲信号，所以原始信号也可被分解成这些信号，但这只是说可以这样做，却是没有用的。  
下面是实域傅立叶变换的BASIC语言代码：

![8394323_1298213514xzg1.jpg][]  
  
到此为止，我们对傅立叶变换便有了感性的认识了吧。但要记住，这只是在实域上的离散傅立叶变换，其中虽然也用到了复数的形式，但那只是个替代的形式，并无实际意义，现实中一般使用的是复数形式的离散傅立叶变换，且**快速傅立叶变换**是根据复数离散傅立叶变换来设计算法的，在后面我们先来复习一下有关复数的内容，然后再在理解实域离散傅立叶变换的基础上来理解复数形式的离散傅立叶变换。**更多见下文**：[**十、从头到尾彻底理解傅里叶变换算法、下**][Link 3]（July、dznlong）

**经典[算法][Link 6]研究系列：十、从头到尾彻底理解傅里叶变换算法、下**

作者：July、dznlong 二零一一年二月二十二日

推荐阅读：***The Scientist and Engineer's Guide to Digital Signal Processing***，By Steven W. Smith, Ph.D。**此书地址**：[http://www.dspguide.com/pdfbook.htm][http_www.dspguide.com_pdfbook.htm]。  
\------------  
从头到尾彻底理解傅里叶变换算法、上  
前言  
第一部分、 DFT  
第一章、傅立叶变换的由来  
第二章、实数形式离散傅立叶变换（Real DFT）

**从头到尾彻底理解傅里叶变换算法、下  
第三章、复数  
第四章、复数形式离散傅立叶变换**

前期回顾，在上一篇：[**十、从头到尾彻底理解傅里叶变换算法、上**][Link 7]里，我们讲了傅立叶变换的由来、和实数形式离散傅立叶变换（Real DFT）俩个问题，  
本文接上文，着重讲下复数、和复数形式离散傅立叶变换等俩个问题。

**第三章、复数**

复数扩展了我们一般所能理解的数的概念，复数包含了实数和虚数两部分，利用复数的形式可以把由两个变量表示的表达式变成由一个变量(复变量)来表达，使得处理起来更加自然和方便。  
我们知道傅立叶变换的结果是由两部分组成的，使用复数形式可以缩短变换表达式，使得我们可以单独处理一个变量（这个在后面的描述中我们就可以更加确切地知道），而且快速傅立叶变换正是基于复数形式的，所以几乎所有描述的傅立叶变换形式都是复数的形式。  
但是复数的概念超过了我们日常生活中所能理解的概念，要理解复数是较难的，所以我们在理解复数傅立叶变换之前，先来专门复习一下有关复数的知识，这对后面的理解非常重要。  
  
**一、 复数的提出**  
  
在此，先让我们看一个物理实验：把一个球从某点向上抛出，然后根据初速度和时间来计算球所在高度，这个方法可以根据下面的式子计算得出：

![8394323_1298213513lll8.jpg][]

其中h表示高度，g表示重力加速度(9.8m/s2)，v表示初速度，t表示时间。现在反过来，假如知道了高度，要求计算到这个高度所需要的时间，这时我们又可以通过下式来计算：

![8394323_1298213513oAs8.jpg][]

（多谢JERRY\_PRI提出：

1、根据公式h=-(gt2/2)+Vt（gt后面的2表示t的平方），我们可以讨论最终情况，也就是说小球运动到最高点时，v=gt，所以，可以得到t=sqt(2h/g)  
且在您给的公式中，根号下为1-(2h)/g，化成分数形式为(g-2h)/g，g和h不能直接做加减运算。

2、g是重力加速度，单位是m/s2，h的单位是m，他们两个相减的话在物理上没有意义，而且使用您给的那个公式反向回去的话推出的是h=-(gt2/2)+gt啊（gt后面的2表示t的平方）。

3、直接推到可以得出t=v/g±sqt((v2-2hg)/g2)（v和g后面的2都表示平方），那么也就是说当v2<2hg时会产生复数，但是如果从实际的v2是不可能小于2hg的，所以我感觉复数不能从实际出发去推到，只能从抽象的角度说明一下。

）  
经过计算我们可以知道，当高度是3米时，有两个时间点到达该高度：球向上运动时的时间是0.38秒，球向下运动时的时间是1.62秒。但是如果高度等于10时，结果又是什么呢？根据上面的式子可以发现存在对负数进行开平方运算，我们知道这肯定是不现实的。  
第一次使用这个不一般的式子的人是意大利数学家Girolamo Cardano（1501-1576），两个世纪后，德国伟大数学家Carl Friedrich Gause（1777-1855）提出了复数的概念，为后来的应用铺平了道路，他对复数进行这样表示：复数由实数（real）和虚数(imaginary)两部分组成，虚数中的根号负1用i来表示（在这里我们用j来表示，因为i在电力学中表示电流的意思）。  
  
我们可以把横坐标表示成实数，纵坐标表示成虚数，则坐标中的每个点的向量就可以用复数来表示，如下图：  
![8394323_1298213513D3wn.jpg][]  
上图中的ABC三个向量可以表示成如下的式子：  
  
A = 2 + 6j  
B = -4 – 1.5j  
C = 3 – 7j  
  
这样子来表达方便之处在于运用一个符号就能把两个原来难以联系起来的数组合起来了，不方便的是我们要分辨哪个是实数和哪个是虚数，我们一般是用Re( )和Im( )来表示实数和虚数两部分，如：  
  
Re A = 2 Im A = 6  
Re B = -4 Im B = -1.5  
Re C = 3 Im C = -7  
  
复数之间也可以进行加减乘除运算：  
![8394323_12982135129xuU.jpg][]  
  
这里有个特殊的地方是j2等于-1，上面第四个式子的计算方法是把分子和分母同时乘以c – dj，这样就可消去分母中的j了。  
  
复数也符合代数运算中的交换律、结合律、分配律：  
  
A B = B A  
(A + B) + C = A + (B + C)  
A(B + C) = AB + AC  
  
  
**二、 复数的极坐标表示形式**  
  
前面提到的是运用直角坐标来表示复数，其实更为普遍应用的是极坐标的表示方法，如下图：

![8394323_1298213512YD5U.jpg][]  
  
上图中的M即是数量积(magnitude)，表示从原点到坐标点的距离，θ是相位角(phase angle)，表示从X轴正方向到某个向量的夹角，下面四个式子是计算方法：  
![8394323_1298213511YE8y.jpg][]  
我们还可以通过下面的式子进行极坐标到直角坐标的转换：  
  
**a + jb = M (cosθ + j sinθ)**

上面这个等式中左边是直角坐标表达式，右边是极坐标表达式。  
  
还有一个更为重要的等式——欧拉等式（欧拉，瑞士的著名数学家，Leonhard Euler，1707-1783）：  
 **ejx = cos x + j sin x**  
  
这个等式可以从下面的级数变换中得到证明：

![8394323_12982135117dZz.jpg][]

上面中右边的两个式子分别是cos(x)和sin(x)的泰勒(Taylor)级数。

这样子我们又可以把复数的表达式表示成指数的形式了：  
  
 **a + jb = M ejθ （这便是复数的两个表达式）**  
  
指数形式是数字信号处理中数学方法的支柱，也许是因为用指数形式进行复数的乘除运算极为简单的缘故吧：

![8394323_1298213511Sn1S.jpg][]

**三、复数是数学分析中的一个工具**  
  
为什么要使用复数呢？其实它只是个工具而已，就如钉子和锤子的关系，复数就象那锤子，作为一种使用的工具。我们把要解决的问题表达成复数的形式（因为有些问题用复数的形式进行运算更加方便），然后对复数进行运算，最后再转换回来得到我们所需要的结果。  
  
有两种方法使用复数，一种是用复数进行简单的替换，如前面所说的向量表达式方法和前一节中我们所讨论的实域DFT，另一种是更高级的方法：数学等价(mathematical equivalence)，复数形式的傅立叶变换用的便是数学等价的方法，但在这里我们先不讨论这种方法，这里我们先来看一下用复数进行替换中的问题。  
  
用复数进行替换的基本思想是：把所要分析的物理问题转换成复数的形式，其中只是简单地添加一个复数的符号j，当返回到原来的物理问题时，则只是把符号j去掉就可以了。  
  
有一点要明白的是并不是所有问题都可以用复数来表示，必须看用复数进行分析是否适用，有个例子可以看出用复数来替换原来问题的表达方式明显是谬误的：假设一箱的苹果是5美元，一箱的桔子是10美元，于是我们把它表示成 5 + 10j，有一个星期你买了6箱苹果和2箱桔子，我们又把它表示成6 + 2j，最后计算总共花的钱是(5 + 10j)(6 + 2j) = 10 + 70j，结果是买苹果花了10美元的，买桔子花了70美元，这样的结果明显是错了，所以复数的形式不适合运用于对这种问题的解决。  
  
**四、用复数来表示正余弦函数表达式**  
  
对于象M cos (ωt + φ)和A cos(ωt ) + B sin(ωt )表达式，用复数来表示，可以变得非常简洁，对于直角坐标形式可以按如下形式进行转换：  
![8394323_1298213511RH49.jpg][]  
上式中余弦幅值A经变换生成a，正弦幅值B的相反数经变换生成b：A <=> a，B<=> -b，但要注意的是，这不是个等式，只是个替换形式而已。  
  
对于极坐标形式可以按如下形式进行转换：  
![8394323_1298213510p4ct.jpg][]  
  
上式中，M <=> M，θ<=>φ。  
这里虚数部分采用负数的形式主要是为了跟复数傅立叶变换表达式保持一致，对于这种替换的方法来表示正余弦，符号的变换没有什么好处，但替换时总会被改变掉符号以跟更高级的等价变换保持形式上的一致。  
  
在离散信号处理中，运用复数形式来表示正余弦波是个常用的技术，这是因为利用复数进行各种运算得到的结果跟原来的正余弦运算结果是一致的，但是，我们要小心使用复数操作，如加、减、乘、除，有些操作是不能用的，如两个正弦信号相加，采用复数形式进行相加，得到的结果跟替换前的直接相加的结果是一样的，但是如果两个正弦信号相乘，则采用复数形式来相乘结果是不一样的。幸运的是，我们已严格定义了正余弦复数形式的运算操作条件：  
1、参加运算的所有正余弦的频率必须是一样的；  
2、运算操作必须是线性的，如两个正弦信号可以进行相加减，但不能进行乘除，象信号的放大、衰减、高低通滤波等系统都是线性的，象平方、缩短、取限等则不是线性的。要记住的是卷积和傅立叶分析也只有线性操作才可以进行。  
  
下图是一个相量变换(我们把正弦或余弦波变成复数的形式称为相量变换，Phasor transform)的例子，一个连续信号波经过一个线性处理系统生成另一个信号波，从计算过程我们可以看出采用复数的形式使得计算变化十分的简洁：  
![8394323_1298213510caW4.jpg][]  
在第二章中我们描述的实数形式傅立叶变换也是一种替换形式的复数变换，但要注意的是那还不是复数傅立叶变换，只是一种代替方式而已。下一章、即，第四章，我们就会知道复数傅立叶变换是一种更高级的变换，而不是这种简单的替换形式。

**第四章、复数形式离散傅立叶变换**

复数形式的离散傅立叶变换非常巧妙地运用了复数的方法，使得傅立叶变换变换更加自然和简洁，它并不是只是简单地运用替换的方法来运用复数，而是完全从复数的角度来分析问题，这一点跟实数DFT是完全不一样的。  
  
**一、 把正余弦函数表示成复数的形式**

通过欧拉等式可以把正余弦函数表示成复数的形式：  
  
cos( x ) = 1/2 e j(-x) + 1/2 ejx  
sin( x ) = j (1/2 e j(-x) - 1/2 ejx)

从这个等式可以看出，如果把正余弦函数表示成复数后，它们变成了由正负频率组成的正余弦波，相反地，一个由正负频率组成的正余弦波，可以通过复数的形式来表示。  
  
我们知道，在实数傅立叶变换中，它的频谱是0 ~ π(0 ~ N/2),但无法表示-π~ 0的频谱，可以预见，如果把正余弦表示成复数形式，则能够把负频率包含进来。  
  
**二、 把变换前后的变量都看成复数的形式**  
  
复数形式傅立叶变换把原始信号x\[n\]当成是一个用复数来表示的信号，其中实数部分表示原始信号值，虚数部分为0，变换结果X\[k\]也是个复数的形式，但这里的虚数部分是有值的。  
在这里要用复数的观点来看原始信号，是理解复数形式傅立叶变换的关键（如果有学过复变函数则可能更好理解，即把x\[n\]看成是一个复数变量，然后象对待实数那样对这个复数变量进行相同的变换）。  
  
**三、 对复数进行相关性算法（正向傅立叶变换）**  
  
从实数傅立叶变换中可以知道，我们可以通过原始信号乘以一个正交函数形式的信号，然后进行求总和，最后就能得到这个原始信号所包含的正交函数信号的分量。

现在我们的原始信号变成了复数，我们要得到的当然是复数的信号分量，我们是不是可以把它乘以一个复数形式的正交函数呢？答案是肯定的，正余弦函数都是正交函数，变成如下形式的复数后，仍旧还是正交函数（这个从正交函数的定义可以很容易得到证明）：

cos x + j sin x, cos x – j sin x，……  
  
这里我们采用上面的第二个式子进行相关性求和，为什么用第二个式子呢?，我们在后面会知道，正弦函数在虚数中变换后得到的是负的正弦函数，这里我们再加上一个负号，使得最后的得到的是正的正弦波，根据这个于是我们很容易就可以得到了复数形式的**DFT正向变换等式**：

![8394323_12983812909Ndf.jpg][]  
这个式子很容易可以得到欧拉变换式子：  
![8394323_1298381290l3Fp.jpg][]  
  
其实我们是为了表达上的方便才用到欧拉变换式，在解决问题时我们还是较多地用到正余弦表达式。  
  
对于上面的等式，我们要清楚如下几个方面（也是区别于实数DFT的地方）：  
1、X\[k\]、x\[n\]都是复数，但x\[n\]的虚数部分都是由0组成的，实数部分表示原始信号；  
2、k的取值范围是0 ~ N-1 (也可以表达成0 ~ 2π)，其中0 ~ N/2（或0 ~ π）是正频部分，

N/2 ~ N-1（π~ 2π）是负频部分，由于正余弦函数的对称性，所以我们把 –π~ 0表示成π~ 2π，这是出于计算上方便的考虑。  
3、其中的j是一个不可分离的组成部分，就象一个等式中的变量一样，不能随便去掉，去掉之后意义就完全不一样了，但我们知道在实数DFT中，j只是个符号而已，把j去掉，整个等式的意义不变；  
4、下图是个连续信号的频谱，但离散频谱也是与此类似的，所以不影响我们对问题的分析：  
![8394323_1298381289adoZ.jpg][]

上面的频谱图把负频率放到了左边，是为了迎合我们的思维习惯，但在实际实

现中我们一般是把它移到正的频谱后面的。

从上图可以看出，时域中的正余弦波（用来组成原始信号的正余弦波）在复数DFT的频谱中被分成了正、负频率的两个组成部分，基于此等式中前面的比例系数是1/N（或1/2π），而不是2/N，这是因为现在把频谱延伸到了2π,但把正负两个频率相加即又得到了2/N,又还原到了实数DFT的形式，这个在后面的描述中可以更清楚地看到。

由于复数DFT生成的是一个完整的频谱，原始信号中的每一个点都是由正、负两个频率组合而成的，所以频谱中每一个点的带宽是一样的，都是1/N，相对实数DFT，两端带宽比其它点的带宽少了一半；复数DFT的频谱特征具有周期性：-N/2 ~ 0与N/2 ~ N-1是一样的，实域频谱呈偶对称性（表示余弦波频谱），虚域频谱呈奇对称性（表示正弦波频谱）。

**四、 逆向傅立叶变换**  
  
假设我们已经得到了复数形式的频谱X\[k\]，现在要把它还原到复数形式的原始信号x\[n\]，当然应该是把X\[k\]乘以一个复数，然后再进行求和，最后得到原始信号x\[n\]，这个跟X\[k\]相乘的复数首先让我们想到的应该是上面进行相关性计算的复数：

cos(2πkn/N) – j si(2πkn/N)，  
  
但其中的负号其实是为了使得进行逆向傅立叶变换时把正弦函数变为正的符号，因为虚数j的运算特殊性，使得原来应该是正的正弦函数变为了负的正弦函数（我们从后面的推导会看到这一点），所以这里的负号只是为了纠正符号的作用，在进行逆向DFT时，我们可以把负号去掉，于是我们便得到了这样的逆向DFT变换等式：

**x\[n\] = X\[k\] (cos(2πkn/N) + j sin(2πkn/N))**

我们现在来分析这个式子，会发现这个式其实跟实数傅立叶变换是可以得到一样结果的。我们先把X\[k\]变换一下：

X\[k\] = Re X\[k\] + j Im X\[k\]

这样我们就可以对x\[n\]再次进行变换，如：

x\[n\] = (Re X\[k\] + j Im X\[k\]) (cos(2πkn/N) + j sin(2πkn/N))

= ( Re X\[k\] cos(2πkn/N) + j Im X\[k\] cos(2πkn/N) +j Re X\[k\] sin(2πkn/N) - Im X\[k\] sin(2πkn/N) )

= ( Re X\[k\] （cos(2πkn/N) + j sin(2πkn/N)） + ---------------------(1)

Im X\[k\] （ - sin(2πkn/N) + j cos(2πkn/N))） ---------------(2)

这时我们就把原来的等式分成了两个部分，第一个部分是跟实域中的频谱相乘，第二个部分是跟虚域中的频谱相乘，根据频谱图我们可以知道，Re X\[k\]是个偶对称的变量，Im X\[k\]是个奇对称的变量，即

Re X\[k\] = Re X\[- k\]  
Im X\[k\] = - Im X\[-k\]

但k的范围是0 ~ N-1，0~N/2表示正频率，N/2~N-1表示负频率，为了表达方便我们把N/2~N-1用-k来表示，这样在从0到N-1的求和过程中对于(1)和(2)式分别有N/2对的k和-k的和，对于（1）式有：

Re X\[k\] (cos(2πkn/N) + j sin(2πkn/N)) + Re X\[- k\] (cos( - 2πkn/N) + j sin( -2πkn/N))

根据偶对称性和三角函数的性质，把上式化简得到：  
  
Re X\[k\] (cos(2πkn/N) + j sin(2πkn/N)) + Re X\[ k\] (cos( 2πkn/N) - j sin( 2πkn/N))

这个式子最后的结果是：

2 Re X\[ k\] cos(2πkn/N)。  
  
再考虑到求Re X\[ k\]等式中有个比例系数1/N，把1/N乘以2，这样的结果不就是跟实数DFT中的式子一样了吗？  
  
对于(2)式，用同样的方法，我们也可以得到这样的结果：

\-2 Im X\[k\] sin(2πkn/N)

注意上式前面多了个负符号，这是由于虚数变换的特殊性造成的，当然我们肯定不能把负符号的正弦函数跟余弦来相加，还好，我们前面是用cos(2πkn/N) – j sin(2πkn/N)进行相关性计算，得到的Im X\[k\]中有个负的符号，这样最后的结果中正弦函数就没有负的符号了，这就是为什么在进行相关性计算时虚数部分要用到负符号的原因（我觉得这也许是复数形式DFT美中不足的地方，让人有一种拼凑的感觉）。  
  
从上面的分析中可以看出，实数傅立叶变换跟复数傅立叶变换，在进行逆变换时得到的结果是一样的，只不过是殊途同归吧。本文完。（July、dznlong）

**本人July对本博客所有任何文章、内容和资料享有版权。  
转载务必注明作者本人及出处，并通知本人。二零一一年二月二十二日。**

[Link 1]: http://www.xuebuyuan.com/category/%E7%AE%97%E6%B3%95
[http_www.dspguide.com_pdfbook.htm]: http://www.dspguide.com/pdfbook.htm
[http_blog.csdn.net_dznlong]: http://blog.csdn.net/dznlong
[Link 2]: http://www.it-crazy.com/
[Link 3]: http://blog.csdn.net/v_JULY_v/archive/2011/02/22/6200945.aspx
[8394323_1298253739kqo6.jpg]: /images/20220715/72f4ef28d14242bcbb8d6436b6867c06.png
[8394323_1298253739WHO6.jpg]: /images/20220715/a3666c9342704ef1a9d2d26df802668b.png
[8394323_1298255200A7rr.jpg]: /images/20220715/7a39ca047ec047c0a40f2d88850abcc9.png
[8394323_1298253736OmpV.jpg]: /images/20220715/3523c558c8104eca9082b28acd1ff0e7.png
[8394323_1298253736ntYg.jpg]: /images/20220715/ccddb01ca1d140488aec744bd895c573.png
[8394323_12982537341B19.jpg]: /images/20220715/5d2422da03ae463db6d4e2478e713e8f.png
[8394323_1298253734CxmE.jpg]: /images/20220715/c7f4d997f06d42d6835077c48c91f6fa.png
[8394323_1298253734oZvG.jpg]: /images/20220715/9e30c55c9eef4f2396daf7259beceea7.png
[8394323_1298266276GaFF.jpg]: /images/20220715/c6531f44eb7946f1b7247bf0f06ee685.png
[8394323_1298213518qfT8.jpg]: /images/20220715/f8344998673d433d8e2644022c428cc7.png
[8394323_1298213518zwwp.jpg]: /images/20220715/da1aa53626264fabb054f8ac422936b7.png
[8394323_12982135185EPi.jpg]: /images/20220715/e91a0f38174e4c369c2acbbf35939fa5.png
[8394323_12982135172NjF.jpg]: /images/20220715/ef8ea495e2c54923905afc8655104d20.png
[Link 4]: http://www.xuebuyuan.com/
[8394323_1298213517sJBS.jpg]: /images/20220715/12a82847526d49feb93097ae19e1aba0.png
[8394323_1298213516D7V1.jpg]: /images/20220715/791010c068ed4e3eab98c2967171904f.png
[8394323_1298213516k4UT.jpg]: /images/20220715/cae198ad9ff4414aa57a25c6c889a6e2.png
[8394323_1298213516KUWW.jpg]: /images/20220715/e54d3d6b669c44a3b2bee5644ebed2f6.png
[8394323_1298213515r5dm.jpg]: /images/20220715/15dd22140b984951ade84b9ac8783175.png
[8394323_12982135155Kvg.jpg]: /images/20220715/2ff0d1e2a4fd412d9669e7625249e2bc.png
[8394323_12982135156Bee.jpg]: /images/20220715/704f10b36b204a1c9a69e94cf8d7d9ed.png
[8394323_1298270500e82R.jpg]: /images/20220715/c1a566a346894254b593a0b866dd2b16.png
[8394323_1298272202VxKp.jpg]: /images/20220715/17b6bccfd0dc42ea995f2f2a6a55fd2d.png
[Link 5]: http://www.carcube.com.cn/
[8394323_1298213514wILB.jpg]: /images/20220715/ee0e619d234b4db0be74a84b2ed552b0.png
[8394323_12982135144LLa.jpg]: /images/20220715/0359a5f9024045f0815a0ea76f26380e.png
[8394323_1298213514xzg1.jpg]: /images/20220715/8a519f1fc95f4907902649c408d1b54a.png
[Link 6]: http://lib.csdn.net/base/datastructure
[Link 7]: http://blog.csdn.net/v_JULY_v/archive/2011/02/20/6196862.aspx
[8394323_1298213513lll8.jpg]: /images/20220715/71a29a3b77ff40059b09f06ce34cf05d.png
[8394323_1298213513oAs8.jpg]: /images/20220715/3f4e63bd56004a76913b784030461e99.png
[8394323_1298213513D3wn.jpg]: /images/20220715/702f3565b81a4de6b5aa0d7f6a14242c.png
[8394323_12982135129xuU.jpg]: /images/20220715/442021baa6784df280767531bb43b3a8.png
[8394323_1298213512YD5U.jpg]: /images/20220715/b4f11059f81b4e54b5e3607cce10a713.png
[8394323_1298213511YE8y.jpg]: /images/20220715/986ebf3fb9f74d9082842c5338091640.png
[8394323_12982135117dZz.jpg]: /images/20220715/834c44e40a8048ff9d35a04b76d78db4.png
[8394323_1298213511Sn1S.jpg]: /images/20220715/f59e0f748f094705a0c614bb2263afca.png
[8394323_1298213511RH49.jpg]: /images/20220715/f49ad2d53d614610804251392a36c2f6.png
[8394323_1298213510p4ct.jpg]: /images/20220715/db552aaee0f94082b92d72be352bc347.png
[8394323_1298213510caW4.jpg]: /images/20220715/7f29afe7cb2a434b915a1684ddc0510a.png
[8394323_12983812909Ndf.jpg]: /images/20220715/041f0b8f27324a9e81663e4ddaed5e3f.png
[8394323_1298381290l3Fp.jpg]: /images/20220715/4a9aedf7a1ee4c8d808d93f9c61d351d.png
[8394323_1298381289adoZ.jpg]: /images/20220715/246d094d87404935b75930dee59adc86.png