数字信号处理——Python实现快速傅里叶变换FFT

文章首发于我的个人博客

1、FFT背景

快速傅里叶变换（FFT）是离散傅里叶变换（DFT）的快速算法，它是根据离散傅里叶的奇、偶、虚、实等特性，在DFT的基础上进行改进获得的。它对傅里叶变换的理论没有新的发现，但它的出现让离散傅里叶变换在计算机系统中得到了广泛的应用。

设x(n)为N项的复数序列，对其进行DFT，任一X(m)的计算都需要N次复数乘法和N-1次复数加法，因此要求出N项复序列的X(m)，大约需要 N 2 N^2 N2次运算。当N逐渐增大时，运算量将呈指数式增长，这在实际应用中是一场灾难。而在FFT中，利用 W N W_N WN的周期性和对称性，把一个N项序列（设N=2k）分为两个N/2项的子序列，每个N/2点DFT变换需要 ( N 2 ) 2 (\frac{N}{2})^2 (2N)2次运算，再用N点运算把两个N/2点的DFT变换组合成一个N点的DFT变换。这样的变换总运算次数为 N + 2 ( N 2 ) 2 = N + N 2 2 N+2(\frac{N}{2})^2=N+\frac{N^2}{2} N+2(2N)2=N+2N2。如果将这种一分为二的思想一直继续下去，知道分成两两一组的DFT运算单元，那么N点的DFT变换就只需要Nlog2N次运算，当N为1024点时，运算量仅为10240次，是使用DFT算法的1%，点数越多，运算量的减少就越显著，这是FFT的优越性。

2、FFT推导

详细的推导在这里我不再重复，我推荐去看程佩青版本的《数字信号处理》基2DIT-FFT部分。

在此，我只通过一个8点DFT流图来说明FFT算法的核心原理，在我使用Python实现FFT代码的过程中，也是先实现上图所示的FFT流图，并不断地根据代码输出结果与正确结果之间的差别来DEBUG代码，不断完善，最终拓展到2048点的FFT。

实现基2DIT-FFT主要有三个步骤：

补零
抽取奇偶序列（倒序）
三个循环
- L级蝶形运算
- 2 L − m 2^{L-m} 2L−m个蝶形运算块
- 每个蝶形运算块中有 2 m − 1 2^{m-1} 2m−1次蝶形运算

首先，将序列输入FFT算法之前，首先要判断序列长度是否满足 N = 2 L N=2^L N=2L，如果不满足情况，需要对输入序列进行补零操作，然后再送入算法。其次，要使最终输出FFT结果为正序排列，需要对输入的数据进行奇偶抽取（倒序）。倒序是指将指数n。倒序是指将数n的二进制码颠倒后的数 n ^ \hat{n} n^。例如当N=8时，采用三位二进制码， n = ( 110 ) = 6 n=(110)=6 n=(110)=6，则 n ^ = ( 011 ) = 3 \hat{n}=(011)=3 n^=(011)=3。因而要将输入序列x(6)和x(3)互换，即做变址运算。N=8时变址结果如下图（与图1中输入序列顺序一致）

接下来，我以图1为例重点说明FFT运算的流程。首先规定，图1为8点FFT，共有L=3级蝶形运算，当前进行蝶形运算的级用m表示。

以第一级m=1为例，第一级需要 2 L − m = 2 3 − 1 = 4 2^{L-m}=2^{3-1}=4 2L−m=23−1=4个蝶形运算块，每个蝶形运算块中有 2 m − 1 = 2 1 − 1 = 1 2^{m-1}=2^{1-1}=1 2m−1=21−1=1次蝶形运算。从图1中也可以明显看出，第一级被分为4个蝶形运算块，每个蝶形运算块中有1次蝶形运算。

以第二级m=2为例，第二级需要 2 L − m = 2 3 − 2 = 2 2^{L-m}=2^{3-2}=2 2L−m=23−2=2个蝶形运算块，每个蝶形运算块中有 2 m − 1 = 2 2 − 1 = 2^{m-1}=2^{2-1}= 2m−1=22−1=次蝶形运算。从图1中也可以明显看出，第二级被分为2个蝶形运算块，每个蝶形运算块中有2次蝶形运算。

以此类推…

在进行第m级蝶形运算时，我们需要用到以下公式：
X m ( k ) = X m − 1 ( k ) + W N r X m − 1 ( k + 2 m − 1 ) X_m(k)=X_{m-1}(k)+W_N^rX_{m-1}(k+2^{m-1}) Xm(k)=Xm−1(k)+WNrXm−1(k+2m−1)

X m ( k + 2 m − 1 ) = X m − 1 ( k ) − W N r X m − 1 ( k + 2 m − 1 ) X_m(k+2^{m-1})=X_{m-1}(k)-W_N^rX_{m-1}(k+2^{m-1}) Xm(k+2m−1)=Xm−1(k)−WNrXm−1(k+2m−1)

从公式中我们可以看出，在进行蝶形运算时，旋转因子 W N r W_N^r WNr的计算是至关重要的，其中r的求法便是核心。具体步骤如下：

将地址k乘以 2 L − m 2^{L-m} 2L−m（即左移L-m位）：

将k变为L位二进制数
将此二进制数左移L-m位，右边空出的位置补零，结果即为r的值。

以第二级m=2为例，将k=0写成3位二进制数，然后左移（3-2）=1位，右边补零得到r=0。将k=1写成3位二进制数，然后左移（3-2）=1位，右边补零得到r=2。与图1中示意图结果相比，一致。

至此，我们便完成了FFT中最关键步骤的解释。在熟悉了FFT公式推导原理的基础上，再加上述我用8点FFT梳理的FFT算法关键点，你便可以写出FFT的代码了。

3、FFT核心代码

def fft_transform(Xn, N=None):if N == None:N = len(Xn)xx = Xn.copy()# 首先判断序列长度是否满足2的幂次，若不满足，补零if math.modf(math.log2(N))[0] != 0:N_ = int(math.pow(2,math.modf(math.log2(N))[1]+1))print(N_, N)xx = np.pad(xx,(0,N_-N),mode="constant",constant_values=0)N = N_print("xx:",len(xx))# 对输入序列进行倒序Xn_copy = reverse_order(xx, N)# 提前计算旋转因子，在进行蝶形运算时，根据r的值调用以计算的旋转因子即可rotation_factor = get_rotation_factor(N)sum = 0# 共有L级蝶形运算L = int(math.log2(N))for m in range(L):# 第m级蝶形运算m = m + 1# 每个蝶形运算块中，蝶形运算次数k_sum = int(math.pow(2,m-1))# 每一级蝶形运算块的总数总count来表示# 这里我并没有直接使用2^(L-m)来算出蝶形运算块的总数# 而是使用累计计数的方式，count表示当前蝶形运算时刻所在的是哪个蝶形运算块# 为什么要采用这种方式？因为我们每一级的蝶形运算结果都保存在一个数组内，# 为了方便数组位置寻址，所以采用这种方式，不同的蝶形运算块可以直接使用count计算得到# 该蝶形运算结果应该存放的位置# 如果对这个步骤依旧不是很理解，可以多花点时间仔细研究count = 0while count < N:for k in range(k_sum):# 蝶形运算公式后半部分tt = Xn_copy[count+k+k_sum]*rotation_factor[k<<(L-m)]t_1 = Xn_copy[count+k] + ttt_2 = Xn_copy[count+k] - ttXn_copy[count+k] = t_1Xn_copy[count+k+k_sum] = t_2count = count + int(math.pow(2,m))    # 在这里更新蝶形运算块#print(Xn_copy)return Xn_copy, running_time

4、实现效果

4.1 基2DIT-FFT与DFT效率对比

4.2 基2DIT-FFT与matlab fft运行结果对比

首先，我使用三个模拟频率分别为10HZ、100HZ、200HZ的正弦波序列作为输入序列，然后分别对这三个序列进行FFT变换，并将其时域曲线和频谱曲线可视化。如图7所示为我用Python实现的DIT-FFT算法变换的结果，图8为使用matlab提供的fft()函数对序列进行变换的结果。从最终结果来看，二者的效果相同。

之后，我将这三个频率不同，持续时间相同的正弦波信号混合得到x(n)混合序列，然后对其进行FFT变换。图9为使用DIT-FFT算法变换的结果，图10为使用matlab提供的fft()函数对其进行变换的结果。从最终呈现的效果来看，二者都可以很清晰地看到三个混合的信号在频谱中对应的谱线。

我将刚刚的混合信号x(n)时间长度设置为1。然后分别对其进行128点，256点，512点的FFT变换，并可视化其频率的幅度谱和相位谱。图11为使用DIT-FFT变换的结果，图12为使用matlab提供的fft()函数进行变换的结果。

从结果中可以看出，当使用128点FFT时，此时频率为100HZ和200HZ的频谱发生混叠。当使用256点FFT时，频率为200HZ的频谱发生混叠。当使用512点FFT时，频谱不混叠，可以在幅度谱中看到频率为50HZ、100HZ、200HZ的谱线。

当我使用对输入信号为600点的序列进行DIT-FFT变换时，此时不满足 N = 2 L N=2^L N=2L条件，需要先对输入序列进行补零，然后再进行变换。图13为使用我用Python实现的DIT-FFT算法对补零到到1024点的序列进行变换的结果。图14为将600点输入序列用matlab提供的fft()算法进行变换的结果。从结果中我们可以看到，我自己实现的DIT-FFT算法频谱有明显的拖尾现象，而matlab变换的结果频谱曲线却很干净。经过一番思考之后，我首先将输入matlab fft()算法的序列先进行1024点补零，然后输入fft()算法，最终得到了如图15所示的效果。由此我们可以发现，matlab提供的fft()算法是基于混合基的DIT-FFT。当输入序列满足条 N = 2 L N=2^L N=2L件时，进行基2DIT-FFT算法，如果不满足该条件则使用混合基DIT-FFT算法，而不是进行补零操作。

5、完整代码

如果需要可在码云获取，代码链接：FFT-Implement。（其中fft_implement.py为FFT实现，tools.py为计算旋转因子、输入序列倒序代码）