从头到尾彻底理解傅里叶变换算法、上经典算法研究系列:十、从头到尾彻底理解傅里叶变换算法、上
推荐阅读:The Scientist and Engineer's Guide to Digital Signal Processing,By Steven W. Smith, Ph.D。此书地址:http://www./pdfbook.htm。 博主说明:I、本文中阐述离散傅里叶变换方法,是根据此书:The Scientist and Engineer's Guide to Digital Signal Processing,By Steven W. Smith, Ph.D.而翻译而成的,此书地址:http://www./pdfbook.htm。II、同时,有相当一部分内容编辑整理自dznlong的博客,也贴出其博客地址,向原创的作者表示致敬:http://blog.csdn.net/dznlong 。这年头,真正静下心写来原创文章的人,很少了。 第三章、复数
前言: 那么,到底什么是傅里叶变换算法列?傅里叶变换所涉及到的公式具体有多复杂列? 哦,傅里叶变换原来就是一种变换而已,只是这种变换是从时间转换为频率的变化。这下,你就知道了,傅里叶就是一种变换,一种什么变换列?就是一种从时间到频率的变化或其相互转化。 ok,咱们再来总体了解下傅里叶变换,让各位对其有个总体大概的印象,也顺便看看傅里叶变换所涉及到的公式,究竟有多复杂: 这是将频率域的函数F(ω)表示为时间域的函数f(t)的积分形式。 连续傅里叶变换的逆变换 (inverse Fourier transform)为: 即将时间域的函数f(t)表示为频率域的函数F(ω)的积分。 一般可称函数f(t)为原函数,而称函数F(ω)为傅里叶变换的像函数,原函数和像函数构成一个傅里叶变换对(transform pair)。 除此之外,还有其它型式的变换对,以下两种型式亦常被使用。在通信或是信号处理方面,常以来代换,而形成新的变换对: 或者是因系数重分配而得到新的变换对: 一种对连续傅里叶变换的推广称为分数傅里叶变换(Fractional Fourier Transform)。分数傅里叶变换(fractional Fourier transform,FRFT)指的就是傅里叶变换(Fourier transform,FT)的广义化。 当f(t)为偶函数(或奇函数)时,其正弦(或余弦)分量将消亡,而可以称这时的变换为余弦变换(cosine transform)或正弦变换(sine transform). 另一个值得注意的性质是,当f(t)为纯实函数时,F(?ω) = F*(ω)成立. 傅里叶级数 其中Fn为复幅度。对于实值函数,函数的傅里叶级数可以写成:
离散时域傅里叶变换 离散傅里叶变换 为了在科学计算和数字信号处理等领域使用计算机进行傅里叶变换,必须将函数xn定义在离散点而非连续域内,且须满足有限性或周期性条件。这种情况下,使用离散傅里叶变换(DFT),将函数xn表示为下面的求和形式: 其中Xk是傅里叶幅度。直接使用这个公式计算的计算复杂度为O(n*n),而快速傅里叶变换(FFT)可以将复杂度改进为O(n*lgn)。(后面会具体阐述FFT是如何将复杂度降为O(n*lgn)的。)计算复杂度的降低以及数字电路计算能力的发展使得DFT成为在信号处理领域十分实用且重要的方法。 下面,比较下上述傅立叶变换的4种变体, 如上,容易发现:函数在时(频)域的离散对应于其像函数在频(时)域的周期性。反之连续则意味着在对应域的信号的非周期性。也就是说,时间上的离散性对应着频率上的周期性。同时,注意,离散时间傅里叶变换,时间离散,频率不离散,它在频域依然是连续的。 ok, 本文,接下来,由傅里叶变换入手,后重点阐述离散傅里叶变换、快速傅里叶算法,到最后彻底实现FFT算法,全篇力求通俗易懂、阅读顺畅,教你从头到尾彻底理解傅里叶变换算法。由于傅里叶变换,也称傅立叶变换,下文所称为傅立叶变换,同一个变换,不同叫法,读者不必感到奇怪。 第一部分、DFT 傅立叶是一位法国数学家和物理学家,原名是Jean Baptiste Joseph Fourier(1768-1830), Fourier于1807年在法国科学学会上发表了一篇论文,论文里描述运用正弦曲线来描述温度分布,论文里有个在当时具有争议性的决断:任何连续周期信号都可以由一组适当的正弦曲线组合而成。 当时审查这个论文拉格朗日坚决反对此论文的发表,而后在近50年的时间里,拉格朗日坚持认为傅立叶的方法无法表示带有棱角的信号,如在方波中出现非连续变化斜率。直到拉格朗日死后15年这个论文才被发表出来。 为什么我们要用正弦曲线来代替原来的曲线呢?如我们也还可以用方波或三角波来代替呀,分解信号的方法是无穷多的,但分解信号的目的是为了更加简单地处理原来的信号。 二、傅立叶变换分类
这四种傅立叶变换都是针对正无穷大和负无穷大的信号,即信号的的长度是无穷大的,我们知道这对于计算机处理来说是不可能的,那么有没有针对长度有限的傅立叶变换呢?没有。因为正余弦波被定义成从负无穷小到正无穷大,我们无法把一个长度无限的信号组合成长度有限的信号。 先来看一个变换实例,下图是一个原始信号图像:
9个余弦信号: 9个正弦信号: 把以上所有信号相加即可得到原始信号,至于是怎么分别变换出9种不同频率信号的,我们先不急,先看看对于以上的变换结果,在程序中又是该怎么表示的,我们可以看看下面这个示例图:
第二章、实数形式离散傅立叶变换(Real DFT) 上图中至于每个波的振幅(amplitude)值(Re X[k],Im X[k])是怎么算出来的,这个是DFT的核心,也是最难理解的部分,我们先来看看如何把分解出来的正余弦波合成原始信号(Inverse DFT)。 如果有学过傅立叶级数,对这个等式就会有似曾相识的感觉,不错!这个等式跟傅立叶级数是非常相似的: 当然,差别是肯定是存在的,因为这两个等式是在两个不同条件下运用的,至于怎么证明DFT合成公式,这个我想需要非常强的高等数学理论知识了,这是研究数学的人的工作,对于普通应用者就不需要如此的追根究底了,但是傅立叶级数是好理解的,我们起码可以从傅立叶级数公式中看出DFT合成公式的合理性。 上面a和 b两个图是待检测信号波,图a很明显可以看出是个3个周期的正弦信号波,图b的信号波则看不出是否含有正弦或余弦信号,图c和d都是个3个周期的正弦信号波,图e和f分别是a、b两图跟c、d两图相乘后的结果,图e所有点的平均值是0.5,说明信号a含有振幅为1的正弦信号c,但图f所有点的平均值是0,则说明信号b不含有信号d。这个就是通过信号相关性来检测是否含有某个信号的方法。 第二个式子中加了个负号,是为了保持复数形式的一致,前面我们知道在计算时又加了个负号,所以这只是个形式的问题,并没有实际意义,你也可以把负号去掉,并在计算时也不加负号。 这里有一点必须明白一个正交的概念:两个函数相乘,如果结果中的每个点的总和为0,则可认为这两个函数为正交函数。要确保关联性算法是正确的,则必须使得跟原始信号相乘的信号的函数形式是正交的,我们知道所有的正弦或余弦函数是正交的,这一点我们可以通过简单的高数知识就可以证明它,所以我们可以通过关联的方法把原始信号分离出正余弦信号。当然,其它的正交函数也是存在的,如:方波、三角波等形式的脉冲信号,所以原始信号也可被分解成这些信号,但这只是说可以这样做,却是没有用的。
本人July对本博客所有任何文章、内容和资料享有版权。 从头到尾彻底理解傅里叶变换算法、下经典算法研究系列:十、从头到尾彻底理解傅里叶变换算法、下
推荐阅读:The Scientist and Engineer's Guide to Digital Signal Processing,By Steven W. Smith, Ph.D。此书地址:http://www./pdfbook.htm。 从头到尾彻底理解傅里叶变换算法、下
第三章、复数 复数扩展了我们一般所能理解的数的概念,复数包含了实数和虚数两部分,利用复数的形式可以把由两个变量表示的表达式变成由一个变量(复变量)来表达,使得处理起来更加自然和方便。 其中h表示高度,g表示重力加速度(9.8m/s2),v表示初速度,t表示时间。现在反过来,假如知道了高度,要求计算到这个高度所需要的时间,这时我们又可以通过下式来计算:
上面中右边的两个式子分别是cos(x)和sin(x)的泰勒(Taylor)级数。 这样子我们又可以把复数的表达式表示成指数的形式了:
第四章、复数形式离散傅立叶变换 复数形式的离散傅立叶变换非常巧妙地运用了复数的方法,使得傅立叶变换变换更加自然和简洁,它并不是只是简单地运用替换的方法来运用复数,而是完全从复数的角度来分析问题,这一点跟实数DFT是完全不一样的。 通过欧拉等式可以把正余弦函数表示成复数的形式: 从这个等式可以看出,如果把正余弦函数表示成复数后,它们变成了由正负频率组成的正余弦波,相反地,一个由正负频率组成的正余弦波,可以通过复数的形式来表示。 现在我们的原始信号变成了复数,我们要得到的当然是复数的信号分量,我们是不是可以把它乘以一个复数形式的正交函数呢?答案是肯定的,正余弦函数都是正交函数,变成如下形式的复数后,仍旧还是正交函数(这个从正交函数的定义可以很容易得到证明): cos x + j sin x, cos x – j sin x,…… N/2 ~ N-1(π~ 2π)是负频部分,由于正余弦函数的对称性,所以我们把 –π~ 0表示成π~ 2π,这是出于计算上方便的考虑。 现中我们一般是把它移到正的频谱后面的。 从上图可以看出,时域中的正余弦波(用来组成原始信号的正余弦波)在复数DFT的频谱中被分成了正、负频率的两个组成部分,基于此等式中前面的比例系数是1/N(或1/2π),而不是2/N,这是因为现在把频谱延伸到了2π,但把正负两个频率相加即又得到了2/N,又还原到了实数DFT的形式,这个在后面的描述中可以更清楚地看到。 由于复数DFT生成的是一个完整的频谱,原始信号中的每一个点都是由正、负两个频率组合而成的,所以频谱中每一个点的带宽是一样的,都是1/N,相对实数DFT,两端带宽比其它点的带宽少了一半;复数DFT的频谱特征具有周期性:-N/2 ~ 0与N/2 ~ N-1是一样的,实域频谱呈偶对称性(表示余弦波频谱),虚域频谱呈奇对称性(表示正弦波频谱)。 四、 逆向傅立叶变换 cos(2πkn/N) – j si(2πkn/N), 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] = (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] 但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)) 根据偶对称性和三角函数的性质,把上式化简得到: 这个式子最后的结果是: 2 Re X[ k] cos(2πkn/N)。 -2 Im X[k] sin(2πkn/N) 注意上式前面多了个负符号,这是由于虚数变换的特殊性造成的,当然我们肯定不能把负符号的正弦函数跟余弦来相加,还好,我们前面是用cos(2πkn/N) – j sin(2πkn/N)进行相关性计算,得到的Im X[k]中有个负的符号,这样最后的结果中正弦函数就没有负的符号了,这就是为什么在进行相关性计算时虚数部分要用到负符号的原因(我觉得这也许是复数形式DFT美中不足的地方,让人有一种拼凑的感觉)。
本人July对本博客所有任何文章、内容和资料享有版权。 |
|