全⾯解析傅⽴叶变换(⾮常详细)
前⾔
第⼀部分、  DFT
第⼀章、傅⽴叶变换的由来
第⼆章、实数形式离散傅⽴叶变换(Real DFT)
从头到尾彻底理解傅⾥叶变换算法、下
第三章、复数
第四章、复数形式离散傅⽴叶变换
前⾔: “关于傅⽴叶变换,⽆论是书本还是在⽹上可以很容易到关于傅⽴叶变换的描述,但是⼤都是些故弄⽞虚的⽂章,太过抽象,尽是⼀些让⼈看了就望⽽⽣畏的公式的罗列,让⼈很难能够从感性上得到理解”---dznlong,
那么,到底什么是傅⾥叶变换算法列?傅⾥叶变换所涉及到的公式具体有多复杂列?
傅⾥叶变换(Fourier transform)是⼀种线性的积分变换。因其基本思想⾸先由法国学者傅⾥叶系统地提出,所以以其名字来命名以⽰纪念。
哦,傅⾥叶变换原来就是⼀种变换⽽已,只是这种变换是从时间转换为频率的变化。这下,你就知道了,傅⾥叶就是⼀种变换,⼀种什么变换列?就是⼀种从时间到频率的变化或其相互转化。
ok,咱们再来总体了解下傅⾥叶变换,让各位对其有个总体⼤概的印象,也顺便看看傅⾥叶变换所涉及到的公式,究竟有多复杂:
以下就是傅⾥叶变换的4种变体(摘⾃,)
连续傅⾥叶变换
⼀般情况下,若“傅⾥叶变换”⼀词不加任何限定语,则指的是“连续傅⾥叶变换”。连续傅⾥叶变换将平⽅可积的函数f(t)表⽰成复指数函数的积分或级数形式。
这是将频率域的函数F(ω)表⽰为时间域的函数f(t)的积分形式。
连续傅⾥叶变换的逆变换 (inverse Fourier transform)为:
即将时间域的函数f(t)表⽰为频率域的函数F(ω)的积分。
⼀般可称函数f(t)为原函数,⽽称函数F(ω)为傅⾥叶变换的像函数,原函数和像函数构成⼀个傅⾥叶变换对(transform pair)。
除此之外,还有其它型式的变换对,以下两种型式亦常被使⽤。在通信或是信号处理⽅⾯,常以来代换,⽽形成新的变换对:
或者是因系数重分配⽽得到新的变换对:
⼀种对连续傅⾥叶变换的推⼴称为分数傅⾥叶变换(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)的推⼴,因为积分其实是⼀种极限形式的求和算⼦⽽已。对于周期函数,其傅⾥叶级数是存在的:
其中Fn为复幅度。对于实值函数,函数的傅⾥叶级数可以写成:
其中an和bn是实频率分量的幅度。
离散时域傅⾥叶变换
离散傅⾥叶变换是离散时间傅⾥叶变换(DTFT)的特例(有时作为后者的近似)。DTFT在时域上离散,在频域上则是周期的。DTFT可以被看作是傅⾥叶级数的逆变换。
离散傅⾥叶变换
傅里叶变换公式表信号与系统离散傅⾥叶变换(DFT),是连续傅⾥叶变换在时域和频域上都离散的形式,将时域信号的采样变换为在离散时间傅⾥叶变换(DTFT)频域的采样。在形式上,变换两端(时域和频域上)的序列是有限长的,⽽实际上这两组序列都应当被认为是离散周期信号的主值序列。即使对有限长的离散信号作DFT,也应当将其看作经过周期延拓成为周期信号再作变换。在实际应⽤中通常采⽤快速傅⾥叶变换以⾼效计算DFT。
为了在科学计算和数字信号处理等领域使⽤计算机进⾏傅⾥叶变换,必须将函数xn定义在离散点⽽⾮连续域内,且须满⾜有限性或周期性条件。这种情况下,使⽤离散傅⾥叶变换(DFT),将函数xn表⽰为下⾯的求和形式:
其中Xk是傅⾥叶幅度。直接使⽤这个公式计算的计算复杂度为O(n*n),⽽快速傅⾥叶变换(FFT)可以将复杂度改进为O(n*lgn)。(后⾯会具体阐述FFT是如何将复杂度降为
O(n*lgn)的。)计算复杂度的降低以及数字电路计算能⼒的发展使得DFT成为在信号处理领域⼗分实⽤且重要的⽅法。
下⾯,⽐较下上述傅⽴叶变换的4种变体,
如上,容易发现:函数在时(频)域的离散对应于其像函数在频(时)域的周期性。反之连续则意味着在对应域的信号的⾮周期性。也就是说,时间上的离散性对应着频率上的周期性。同时,注意,离散时间傅⾥叶变换,时间离散,频率不离散,它在频域依然是连续的。
如果,读到此,你不甚明⽩,⼤没关系,不必纠结于以上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):
这四种傅⽴叶变换都是针对正⽆穷⼤和负⽆穷⼤的信号,即信号的的长度是⽆穷⼤的,我们知道这对于计算机处理来说是不可能的,那么有没有针对长度有限的傅⽴叶变换呢?没有。因为正余弦波被定义成从负⽆穷⼩到正⽆穷⼤,我们⽆法把⼀个长度⽆限的信号组合成长度有限的信号。
⾯对这种困难,⽅法是:把长度有限的信号表⽰成长度⽆限的信号。如,可以把信号⽆限地从左右进⾏延伸,延伸的部分⽤零来表⽰,这样,这个信号就可以被看成是⾮周期性离散信号,我们可以⽤到离散时域傅⽴叶变换(DTFT)的⽅法。也可以把信号⽤复制的⽅法进⾏延伸,这样信号就变成了周期性离散信号,这时我们就可以⽤离散傅⽴叶变换⽅法(DFT)进⾏变换。本章我们要讲的是离散信号,对于连续信号我们不作讨论,因为计算机只能处理离散的数值信号,我们的最终⽬的是运⽤计算机来处理信号的。
但是对于⾮周期性的信号,我们需要⽤⽆穷多不同频率的正弦曲线来表⽰,这对于计算机来说是不可能实现的。所以对于离散信号的变换只有离散傅⽴叶变换(DFT)才能被适⽤,对于计算机来说只有离散的和有限长度的数据才能被处理,对于其它的变换类型只有在数学演算中才能⽤到,在计算机⾯前我们只能⽤DFT⽅法,后⾯我们要理解的也正是DFT⽅法。
这⾥要理解的是我们使⽤周期性的信号⽬的是为了能够⽤数学⽅法来解决问题,⾄于考虑周期性信号是从哪⾥得到或怎样得到是⽆意义的。
每种傅⽴叶变换都分成实数和复数两种⽅法,对于实数⽅法是最好理解的,但是复数⽅法就相对复杂许多了,需要懂得有关复数的理论知识,不过,如果理解了实数离散傅⽴叶变换(real DFT),再去理解复数傅⽴叶变换就更容易了,所以我们先把复数的傅⽴叶变换放到⼀边去,先来理解实数傅⽴叶变换,在后⾯我们会先讲讲关于复数的基本理论,然后在理解了实数傅⽴叶变换的基础上再来理解复数傅⽴叶变换。
还有,这⾥我们所要说的变换(transform)虽然是数学意义上的变换,但跟函数变换是不同的,函数变换是符合⼀⼀映射准则的,对于离散数字信号处理(DSP),有许多的变换:傅⽴叶变换、拉普拉斯变换、Z变换、希尔伯特变换、离散余弦变换等,这些都扩展了函数变换的定义,允许输⼊和输出有多种的值,简单地说变换就是把⼀堆的数据变成另⼀堆的数据的⽅法。
三、⼀个关于实数离散傅⽴叶变换(Real DFT)的例⼦
先来看⼀个变换实例,下图是⼀个原始信号图像:
这个信号的长度是16,于是可以把这个信号分解9个余弦波和9个正弦波(⼀个长度为N的信号可以分解成N/2+1个正余弦信号,这是为什么呢?结合下⾯的18个正余弦图,我想从计算机处理精度上就不难理解,⼀个长度为N的信号,最多只能有N/2+1个不同频率,再多的频率就超过了计算机所能所处理的精度范围),如下图:
9个余弦信号:
9个正弦信号:
把以上所有信号相加即可得到原始信号,⾄于是怎么分别变换出9种不同频率信号的,我们先不急,先看看对于以上的变换结果,在程序中⼜是该怎么表⽰的,我们可以看看下⾯这个⽰例图:
上图中左边表⽰时域中的信号,右边是频域信号表⽰⽅法,
从左向右,-->,表⽰正向转换(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个周期,如下图:
上图中⾄于每个波的振幅(amplitude)值(Re X[k],Im X[k])是怎么算出来的,这个是DFT的核⼼,也是最难理解的部分,我们先来看看如何把分解出来的正余弦波合成原始信号(Inverse DFT)。
三、合成运算⽅法(Real Inverse DFT)
DFT合成等式(合成原始时间信号,频率-->时间,逆向变换):
如果有学过傅⽴叶级数,对这个等式就会有似曾相识的感觉,不错!这个等式跟傅⽴叶级数是⾮常相似的:
当然,差别是肯定是存在的,因为这两个等式是在两个不同条件下运⽤的,⾄于怎么证明DFT合成公式,
这个我想需要⾮常强的⾼等数学理论知识了,这是研究数学的⼈的⼯作,对于普通应⽤者就不需要如此的追根究底了,但是傅⽴叶级数是好理解的,我们起码可以从傅⽴叶级数公式中看出DFT合成公式的合理性。
_            _
DFT合成等式中的Im X[k]和Re X[k]跟之前提到的Im X[k]和Re X[k]是不⼀样的,下⾯是转换⽅法(关于此公式的解释,见下⽂):
但k等于0和N/2时,实数部分的计算要⽤下⾯的等式:
上⾯四个式中的N是时域中点的总数,k是从0到N/2的序号。
为什么要这样进⾏转换呢?这个可以从频谱密度(spectral density)得到理解,如下图就是个频谱图:
这是⼀个频谱图,横坐标表⽰频率⼤⼩,纵坐标表⽰振幅⼤⼩,原始信号长度为N(这⾥是
32),经DFT转换后得到的17个频率的频谱,频谱密度表⽰每单位带宽中为多⼤的振幅,那么带宽是怎么计算出来的呢?看上图,除了头尾两个,其余点的所占的宽度是2/N,这个宽度便是每个点的带宽,头尾两个点的带宽是1/N,⽽Im X[k]和Re X[k]表⽰的是频谱密度,即每⼀个单位带宽的振幅⼤⼩,但表⽰2/N(或1/N)带宽的振幅⼤⼩,所以分别应当是Im X[k]和Re X[k]的
2/N(或1/N)。
频谱密度就象物理中物质密度,原始信号中的每⼀个点就象是⼀个混合物,这个混合物是由不同密度的物质组成的,混合物中含有的每种物质的质量是⼀样的,除了最⼤和最⼩两个密度的物质外,这样我们只要把每种物质的密度加起来就可以得到该混合物的密度了,⼜该混合物的质量是单位质量,所以得到的密度值跟该混合物的质量值是⼀样的。
⾄于为什么虚数部分是负数,这是为了跟复数DFT保持⼀致,这个我们将在后⾯会知道这是数学计算上的需要(Im X[k]在计算时就已经加上了⼀个负号(稍后,由下⽂,便可知),再加上负号,结果便是正的,等于没有变化)。
如果已经得到了DFT结果,这时要进⾏逆转换,即合成原始信号,则可按如下步骤进⾏转换:
1、先根据上⾯四个式⼦计算得出的值;
2、再根据DFT合成等式得到原始信号数据。
下⾯是⽤BASIC语⾔来实现的转换源代码:
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

版权声明:本站内容均来自互联网,仅供演示用,请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198,我们将在24小时内删除。