Loading... ## 1、FFT背景 快速傅里叶变换（FFT）是离散傅里叶变换（DFT）的快速算法，它是根据离散傅里叶的奇、偶、虚、实等特性，在DFT的基础上进行改进获得的。它对傅里叶变换的理论没有新的发现，但它的出现让离散傅里叶变换在计算机系统中得到了广泛的应用。 设x(n)为N项的复数序列，对其进行DFT，任一X(m)的计算都需要N次复数乘法和N-1次复数加法，因此要求出N项复序列的X(m)，大约需要$N^2$次运算。当N逐渐增大时，运算量将呈指数式增长，这在实际应用中是一场灾难。而在FFT中，利用$W_N$的周期性和对称性，把一个N项序列（设N=2k）分为两个N/2项的子序列，每个N/2点DFT变换需要$(\frac{N}{2})^2$次运算，再用N点运算把两个N/2点的DFT变换组合成一个N点的DFT变换。这样的变换总运算次数为$N+2(\frac{N}{2})^2=N+\frac{N^2}{2}$。如果将这种一分为二的思想一直继续下去，知道分成两两一组的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^{m-1}$次蝶形运算 首先，将序列输入FFT算法之前，首先要判断序列长度是否满足$N=2^L$，如果不满足情况，需要对输入序列进行补零操作，然后再送入算法。其次，要使最终输出FFT结果为正序排列，需要对输入的数据进行奇偶抽取（倒序）。倒序是指将指数n。倒序是指将数n的二进制码颠倒后的数$\hat{n}$。例如当N=8时，采用三位二进制码，$n=(110)=6$，则$\hat{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^{m-1}=2^{1-1}=1$次蝶形运算。从图1中也可以明显看出，第一级被分为4个蝶形运算块，每个蝶形运算块中有1次蝶形运算。 以第二级`m=2`为例，第二级需要$2^{L-m}=2^{3-2}=2$个蝶形运算块，每个蝶形运算块中有$2^{m-1}=2^{2-1}=$次蝶形运算。从图1中也可以明显看出，第二级被分为2个蝶形运算块，每个蝶形运算块中有2次蝶形运算。 以此类推... 在进行第m级蝶形运算时，我们需要用到以下公式： $$ X_m(k)=X_{m-1}(k)+W_N^rX_{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}) $$ 从公式中我们可以看出，在进行蝶形运算时，旋转因子$W_N^r$的计算是至关重要的，其中r的求法便是核心。具体步骤如下： 将地址k乘以$2^{L-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核心代码 ```python 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 = 0 while 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] + tt t_2 = Xn_copy[count+k] - tt Xn_copy[count+k] = t_1 Xn_copy[count+k+k_sum] = t_2 count = 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$条件，需要先对输入序列进行补零，然后再进行变换。图13为使用我用Python实现的DIT-FFT算法对补零到到1024点的序列进行变换的结果。图14为将600点输入序列用matlab提供的fft()算法进行变换的结果。从结果中我们可以看到，我自己实现的DIT-FFT算法频谱有明显的拖尾现象，而matlab变换的结果频谱曲线却很干净。经过一番思考之后，我首先将输入matlab fft()算法的序列先进行1024点补零，然后输入fft()算法，最终得到了如图15所示的效果。由此我们可以发现，matlab提供的fft()算法是基于混合基的DIT-FFT。当输入序列满足条$N=2^L$件时，进行基2DIT-FFT算法，如果不满足该条件则使用混合基DIT-FFT算法，而不是进行补零操作。    ## 5、完整代码 如果需要可在码云获取，代码链接：[FFT-Implement](https://gitee.com/lucasnan/python_advanced_learning.git)。（其中`fft_implement.py`为FFT实现，`tools.py`为计算旋转因子、输入序列倒序代码） Last modification：March 2nd, 2022 at 02:59 pm © 允许规范转载