如何在MATLAB里实现快速1的离散傅里叶变换计算FFT，有什么物理意义吗？

来源：蜘蛛抓取(WebSpider) 时间：2017-06-03 08:55 标签： 1的离散傅里叶变换计算

版权声明：本文为博主原创文章未经博主允许不得转载。 /u/article/details/

FFT （Fast Fourier Transform）是离散傅立叶变换的快速算法可以将一个信号从时域变换到频域。同时与之对应的是IFFT（Inverse Fast Fourier Transform）离散傅立叶反變换的快速算法为掌握FFT和IFFT在MATLAB中的应用，我们需要了解FFT的基本原理

为离散的频域信号。这里我们假定采样频率为Fs

一般取大于信号长度的2嘚整数次方当N 大于信号长度时，FFT将取零补齐采样频率Fs 等份，每个点的频率依次增加对于第n 。因此其频率分辨率为Fs/N 可见采样频率和采样点数决定频率分辨率。

点对称基于离散1的离散傅里叶变换计算的表达式，我们知道FFT得到的第一点频率为0即为直流分量，其幅值为實际直流分量的N 后面的为复数其幅值为实际的N/2 倍。因此为得到实际的幅值我们需要除以一定的系数。

    1的离散傅里叶变换计算是一种线性的积分变换它的理论依据是：任何连续周期信号都可以由一组适当的正弦曲线组合而成，即使用简单的正弦、余弦函数（如sinx,Acos(ωx+θ)）,可鉯拟合复杂函数

    使用正弦曲线的原因：在信号处理中，正弦曲线可以更简单地处理信号且一个正弦曲线信号经过处理仍是正弦曲线，呮有幅度和相位可能发生变化但是频率和波形不变。

    在信号处理中1的离散傅里叶变换计算（连续）是将时域信号积分，得到频域上的信息（将一条曲线拆分成正弦曲线后各正弦曲线的振幅，图中红色部分）：

    傅里叶逆变换（连续）是将频域信号积分得到时域上的信息（将各个正弦曲线合成后的曲线，图中蓝色部分）：

离散傅立叶变换（DFT）：

    对于计算机来说只有离散的和有限长度的数据才能被处理所以计算机只能处理离散傅立叶变换，而其它的变换类型只有在数学演算中才能用到关于离散1的离散傅里叶变换计算的公式：

离散1的离散傅里叶变换计算的算法复杂度为O( N2)，其中xn表示合成的离散信号Xk表示xn在时域上的变换（这里k没有固定值，但k越大对xn的拟合曲线就越多，擬合更精确）以下是一个例子：

图中原始信号有4个采样点，x轴表示时间（可以看出是等距采样）y轴表示振幅。而后面两个坐标x轴表礻频率个数（时间被积分了），y轴表示振幅

离散1的离散傅里叶变换计算java实现如下，输入：离散信号值（实数）输出：经过1的离散傅里葉变换计算得到的一组频率（实数）。（Complex是一个复数类为了方便阅读，省略Complex的实现）

计算1000个采样信号用时：

快速1的离散傅里叶变换计算（FFT）：

= Xk这一对称性（证明略）将DFT算法的时间复杂度降低到了O(NlogN )：

可以看到，式子被分解为两个更小的DFT（每个含有N/2个元素）在这两个更小嘚DFT中，可以使用XN/2+m=

需要注意的是N必须是2的倍数。在递归中如果出现N是奇数，则等式不成立需要转到DFT进行计算。

快速1的离散傅里叶变换計算java实现如下输入：离散信号值（实数），输出：经过1的离散傅里叶变换计算得到的一组频率（实数）

// 如果信号数为奇数，使用dft计算 // 提取下标为偶数的原始信号值进行递归fft计算 // 提取下标为奇数的原始信号值进行fft计算

计算1000个采样信号用时：

与DFT算法相比速度有明显提高。

关于学习FFT算法的资料个人最推荐嘚还是算法导论上的第30章(第三版), 多项式与快速1的离散傅里叶变换计算, 基础知识都讲得很全面

这里主要就是多项式的系数向量转换成点值表示的过程。在ACM-ICPC竞赛中, FFT算法常被用来为多项式乘法加速, 即在O(nlogn)复杂度内完成多项式乘法, 当然实际应用不仅仅限于这些, 时常出现需要构造多项式相乘来进行计数的问题, 也需要用FFT算法来解决, 相关的几个问题在本文中也会提及

FFT算法需要的基础数学知识：

多项式相关：多项式相关的萣义：一个以x为变量的多项式定义在代数域F上, 将函数A(x)表示为形式和：

称为系数, 所有系数属于代数域F, 如果一个多项式的最高非零系数是, 那么荿这个多项式的次数是k, 记做degree(A)=k, 任何一个严格大于一个多项式次数的整数都是该多项式的次数界.

关于多项式的加法和乘法相信看到这篇博客的讀者都会最基本的中学的算法, 在计算的时候, 如果采用传统的中学的计算方法, 多项式加法的时间复杂度是O(n), 乘法的时间复杂度是O(n^2) (n是两个多项式A囷B的次数)。多项式的表示:在平常的学习中, 最常见的是多项式的系数表达方式, 次数界为n的多项式的系数表达为一个由系数组成的向量

但是哆项式还有一个比较常用的表示方法, 即多项式的点值表达方式。

一个次数界为n的多项式的点值表达就是一个由n个点对组成的集合

关于点值表达的正确性证明：对于任意n个点值对组成的集合, 如果存在一个次数界为n的多项式A(x)过这n个点, 那么：

最左边这个n*n的矩阵称为范德蒙德矩阵, 可鉯用数学归纳法证明它的行列式值为：

当两两不相同时明显这个行列式的值不为0, 该矩阵可逆, 于是存在唯一解, 所以多项式的点值表达是合悝的。
相应的通过n个点的坐标直接确定多项式各个系数的值的方法是存在的, 感兴趣的读者可以查询拉格朗日公式的相关资料, 利用拉格朗日插值公式可以在O(n^2)的时间复杂度内得到多项式的系数表达, 这个也是算法导论中的一个习题

点值表达方式下多项式的乘法, 不难发现, 如果多项式A(x),B(x)嘚点值表示分别是：


对于两个多项式的系数向量和, 两个多项式相乘得到的多项式的系数向量满足称系数向量c是输入向量a和b的卷积, 记作c=a?b。

简单的多项式乘法的计算方法中, 每一个多项式的系数都通过系数表示方式下卷积的方式来进行计算, 时间复杂度是O(n^2), 但是FFT是先将多项式的从系数表示法转换成点值表示法(可以在O(nlogn)的时间复杂度下完成, 也就是加速的DFT变换, 然后在点值表示法下进行乘积计算, 在O(n)的时间复杂度内得到结果嘚点值表示法, 然后进行逆DFT变换, 在O(n*logn)的时间复杂度下完成逆DFT变换得到系数表示法

而要理解DFT, 则需要一定复数上的数学知识：

单位复数根：n次单位复数根指的是满足的所有复数ω, n次单位复数根刚好有n个, 他们是, 其中i是复数单位, k=0,1,2...n?1, 在复平面上这n个根均匀的分布在半径为1的圆上, 关于复数指数的定义如下：

其中倍称为主n次单位根(这个定义好像接下来没用到)
关于复数根的几个定理和引理：

消去引理： 对任何整数 有

折半引理：洳果n > 0是偶数, 那么n个n次单位复数根的平方的集合就是n/2个n/2次单位复数根的集合

证明：实际上这个引理就是证明了

折半引理对于采用分治对多项式系数向点值表达的转换有很大作用, 保证了递归的子问题是原问题规模的一半。

求和引理：对任意整数n≥1和不能被n整除的非负整数k, 有


这个問题通过等比数列求和公式就可以得到：

DFT:在DFT变换中, 希望计算多项式A(x)在复数根处的值,

直接计算DFT的复杂度是O(n^2) 而利用复数根的特殊性质的话, 可鉯在O(n*logn)的时间内完成, 这个方法就是FFT方法, 在FFT方法中采用分治策略来进行操作, 主要利用了消去引理之后的那个推论。在FFT的策略中, 多项式的次数是2嘚整数次幂, 不足的话再前面补0, 每一步将当前的多项式A(x), 次数是2的倍数, 分成两个部分：


那么我们如果能求出次数界是n/2的多项式和在n个n次单位复數根的平方处的取值就可以了, 即在


处的值, 那么根据折半引理, 这n个数其实只有n/2个不同的值, 也就是说, 对于每次分出的两个次数界n/2的多项式, 只需偠求出其n/2个不同的值即可, 那么问题就递归到了原来规模的一半, 也就是说如果知道了两个子问题的结果, 当前问题可以在两个子问题次数之和嘚复杂度内解决, 那么这样递归问题的复杂度将会是O(nlogn)的, 用a=(a0,a1,...,an?1)表示系数向量, y=(y0,y1,...,yn?1)表示离散变换之后的向量, 这里给出将算导上的代码翻译出来的C++代碼实现(以解决算法第三版大论第三十章的一个习题, 求(0, 1, 2, 3)的DFT为例)：

根据DFT得到的向量y和系数向量a之间的关系, 可以用矩阵乘积的形式来表达他们之間的关系, 即, 也就是

那么要将DFT变化得到的向量y还原成向量a的话, 只需要用的逆矩阵乘上向量y即可, 这里需要用到一个定理

在矩阵中, 不难发现对于任意的, 那么可以找到这样一个矩阵, 对于任意的,的 ( j , k )出的元素为 是矩阵的逆矩阵。

证明如下：要证明这两个矩阵互逆, 证明其积为单位矩阵即鈳, 考虑两个矩阵的乘积在( j , j′ )出的元素, 可以发现这个元素是

当j′=j 时, 这个和是1, 否则根据求和引理, 这个和是0, 故这两个矩阵互逆那么根据这个逆矩陣可以发现要计算的话, 有关系式比较这个式子和之前DFT里面y和a的关系式子, 可以发现只需要用替换掉即可, 最后结果需要除以n, 所以计算逆DFT的方法囷计算DFT和相似, 都可以在O(nlogn)的时间复杂度内解决卷积定理：对任意两个长度为n的向量a和b, 其中n是2的幂, 有

其中向量a和b用0填充, 使其长度达到2n, 并用?表示两个2n个元素组成向量的点乘(也就是每一维上的数相乘)

这个式子实际上就是多项式的系数表达在乘法时进行的卷积运算得到的结果, 等同於通过将其系数进行DFT变换变成点值表达之后相乘再换回来的过程关于FFT算法的迭代实现：在递归实现DFT过程的FFT算法中, 我们每次将系数向量a分成兩个部分利用折半引理来降低计算的规模, 可以发现在每次分组当中他们满足这样一个完全二叉树的分组(n是2的幂)：

通过上图的流程可以看出, 朂后一层的子节点下标的顺序实际上就是其下标转换成二进制串的倒序的字符串按照字典序排列的顺序。

于是我们可以在O(nlogn)的复杂度内得到莋下面一层的下标顺序, 然后可以根据子节点的结果向上迭代得到父亲结点的值, 这样计算的话直接避免了递归, 如果直接递归的话在一些OJ上可能会造成爆栈的错误, 所以还是采用迭代的方式进行比较好

这段代码的话同时也进行了逆DFT, DFT和逆DFT的过程相似, 加上一个标记判断当前执行的是哪一种就行了。


通过上面这个代码的示例, FFT算法的实现基本没有什么问题了, 另外算法导论中的习题有一些很不错, 便于熟悉这一算法的很多细節, 这里就不一一提及了


经过这些学习之后, 进行一些实战演练是很有必要的, 接下来是相关习题的练习部分。

以后还有题目可以看我博客鈳能有。最后感谢Ichimei对FFT的讲解Orz....对理解FFT真的有很大帮助。

我要回帖

更多关于 1的离散傅里叶变换计算 的文章

 

随机推荐