0%

数字信号处理之FFT终极理解

本节从离散傅里叶变换的深入理解说起,介绍了FFT的基本原理,其中包括基2、基4的按时间抽选FFT算法和按频率抽取FFT算法

从DFT说起

  • DFT的基本公式:

  • 看这个公式里面有个$\frac{2\pi}{N}$,有没有觉得很熟悉,哦没错,它其实就是数字角频率$w_0$(这个在我的其他笔记中有提到信号处理中的小知识点 | ssy的小天地 (ssy1938010014.github.io)

    • 那么该如何理解它呢?正如在之前笔记中提到的,它代表在复数平面内每个点转过的角度

      image-20231103170339511

    • 图中标出的$\theta$,其实就是数字角频率$w_0$,在对一个正弦波抽样后,得到的点,其实就没有时间的概念了,只是我们自己心里清楚,点与点之间的时间间隔是$T_s$(采样间隔)

    • 那么如何数字角频率和实际模拟信号的角频率是如何对应的呢?即在之前笔记中提到的:$\Omega_0=w_0f_s$,那对于频域来说,旋转因子中的角频率不能只有$w_0$,不然就只能测出$w_0$对应模拟信号的频率,类似DTFT的基频和倍频,这里我们也多取一些角频率$w_0,2w_0,3w_0,…(N-1)w_0$,也就是$kw_0,k=0,1,…,N-1$。理想情况下,如果在频域$kw_0$有个冲击,那么就说明对应时域信号中含有角频率为$kw_0f_s=k\frac{2\pi}{N}f_s=2\pi\cdot \frac{kf_s}{N}$的正弦信号,即频率为$k\frac{f_s}{N}$的正弦信号(其实,这也从公式本质的切入点上,解释了数字k域与实际时域频率的对应关系)

    • 上述一大段的解释可能也有误,但是我自己最极限的理解了,如果不理解也没关系,我觉得这里只需要记住:以角频率$w_0$每隔$T_s$时间在复数域内转动的N点(其在实轴上的投影就是$cos(\frac{2\pi}{N}n)$,其本身就是复数$e^{j\frac{2\pi}{N}n}$),在经过DFT之后,可以得到这个N点数字域对应模拟域的实际频率$f_0$,且$f_0=k\cdot \frac{f_s}{N}$

  • 下面重点解释为啥$x(n)$与$e^{-j\frac{2\pi}{N}nk}$这么个玩意相乘,就能挑选出对应频率(这些解释都偏个人的感性解释,也不一定完全正确)

    • 首先,假设$x(n)=cos(\frac{2\pi}{N}n)$,$x(k)$就等于$k$确定情况下,$n=0\sim N-1$个$x(n)$对应的点与旋转因子中$e^{-j\frac{2\pi}{N}n},n=0\sim N-1$对应元素相乘(例如$x(0)\times e^{-j\frac{2\pi}{N}\cdot 0\cdot k},x(1)\times e^{-j\frac{2\pi}{N}\cdot 1\cdot k},x(2)\times e^{-j\frac{2\pi}{N}\cdot 2\cdot k},…$)。简单来说,这一步骤就是:在$k$(也就是某一频率)下,求信号对这一频率幅频特性曲线的贡献(这么说可能不太规范,通俗来讲,就是如果信号中含有这个频率,那么它的幅频特性曲线中,这个频率对应的幅度就会格外高,也就是说这一步就是在求这个幅度到底有多高)

    • 此时,你会发现,有个离谱的事情发生了$x(n)$是个$cos$,$e^{j\frac{2\pi}{N}n}$是个复数,这咋乘,就现在,我发现了欧拉的伟大,应用欧拉公式,有:

    • 现在是不是cos就能和复数相乘了,为了便于解释,假设$N=8$,那么$e^{j\frac{2\pi}{N}n}$在复平面上的8个点(图1)为:

      Snipaste_2023-11-04_02-41-55

    • $e^{-j\frac{2\pi}{N}n}$在复平面上的8个点(图2)为:

      Snipaste_2023-11-04_02-43-08

    • 当$k=1$时,$W_N^{n}=e^{-j\frac{2\pi}{N}n}$,其在复平面上的8个点(图3)为:

      Snipaste_2023-11-04_02-43-08

      • 图1中8个复数乘以图3中对应的复数(例如图1中n=1的复数乘以图3中n=1中的复数,图1中n=2的复数乘以图3中n=2中的复数),显然,每个点对点的相乘,其均是模值相乘,由于角度相反,那么相乘后相位为0,最终这8个点对点相乘求和的结果就是8,可视化示意图如下:

        image-20231104030201498

      • 同样,图2中8个复数乘以图3中对应的复数,每个点对点结果是模值相乘,旋转角变为原来两倍,对结果求和后,最终为0

        image-20231104031359017

    • 其实,$W_N^{nk}=e^{-j\frac{2\pi}{N}nk}$就是将复平面均匀分成8份,由$nk$决定,具体对应是八份中的哪一份,当$k=2$时,$nk$为$0,2,4,6,8,10,12,14$,由于8个点就转了一圈,即周期时8,那么$nk$实际对应的是$0,2,4,6,0,2,4,6$,其可视化示意如下(图4):

      image-20231104153953253

      • 将图1与图4中n相同的点对应相乘,根据模值相乘,角度相加的原则,最终矢量相加的结果为0,图2乘以图4的结果也为0
    • 同理,当$k=3,4,5,6$时($k=7$的情况后面会解释,存在对称性),图1乘以图4,图2乘以图4的结果均为0(这里没有严格上的数学推导,但我觉得对于理解FFT,自己带数进去算一遍更能理解)

    • 总之,最后我们会发现:只有$e^{j\frac{2\pi}{N}n}$乘以$W_N^{nk}(k=1)$有值,其值为$e^{j\frac{2\pi}{N}n}$每个点的模值(这里都假设为1了)与$W_N^{nk}$对应点的模值(显然$W_N^{nk}$的模都是1)相乘后再求和(还要除以2,因为欧拉公式后$e^{j\frac{2\pi}{N}n}$前的系数为$\frac12$),其实,$e^{j\frac{2\pi}{N}n}$每个点的模值都一样,你想,$cos(\frac{2\pi}{N}n)$和$5cos(\frac{2\pi}{N}n)$用欧拉公式展开之后,不就是$e^{j\frac{2\pi}{N}n}$和$5e^{j\frac{2\pi}{N}n}$的区别么,那么实际上每个点的模值都是一样的,即$(1,1,1,1,1,1,1,1)$和$(5,5,5,5,5,5,5,5)$,那么实际上$cos(\frac{2\pi}{N}n)$乘以$W_N^{nk}$就是$\frac{A\times N}{2}$(其中$A$为正弦信号的幅度,$N$为采样的点数),这个峰值会出现在$k=1$的位置(这里的解释都先不讨论$k=7$的情况)

    • 哦对了,其实这里也很好的解释了采样点数$N$加倍,DFT结果幅值也会加倍

    • 有个细节是,当$k=0$时,很容易得到$cos(\frac{2\pi}{N}n)$乘以$W_N^{nk}$是16(因为不管是$e^{j\frac{2\pi}{N}n}$还是$e^{-j\frac{2\pi}{N}n}$都是一堆111去乘以$W_N^{nk}$的一堆1嘛),所以这也解释了零频对应的峰值是其他频率峰值的两倍

  • 下面解释DFT的共轭对称性,即$X(k)=X^*(N-k)$

    • 那么我们来讨论$k=7$的情况,$W_N^{7n}=e^{-j\frac{2\pi}{N}n\cdot7}$,其在复平面上的8个点(图5)为:

      image-20231104163326522

    • 不难发现,$k=7$时的情况刚好和$k=1$的情况是反过来的,此时是$e^{-j\frac{2\pi}{N}n}$乘以$W_N^{nk}$有值(且也为$\frac{A\times N}{2}$),$e^{j\frac{2\pi}{N}n}$乘以$W_N^{nk}$为0,所以频谱中还会在$k=7$时再出现一个峰值

    • 同时也解释了为什么复数信号的DFT是单边谱,因为复数信号只有$e^{j\frac{2\pi}{N}n}$部分,只存在$e^{j\frac{2\pi}{N}n}$乘以$W_N^{nk}$的结果,故复数信号只有单边谱(换句话来说,实数信号的双边谱多余的那边是$e^{-j\frac{2\pi}{N}n}$造成的)

    • 笑了,上面解释一大堆,其实就是用了旋转因子的对称性:$(W_N^{nk})^*=W_N^{-nk}=W_N^{(N-n)k}$


按时间抽选FFT算法

1.基2-FFT算法

  • 设序列$x(n)$点数为$N=2^L$,$L$为正整数

  • 将序列$x(n)$按$n$的奇偶分为以下两组:

  • 则其$N$点FFT为:

    • 其中有用到旋转因子的约分性$W_N^m=W_{\frac N2}^{\frac m2}$
    • 且$X_1(k)$和$X_2(k)$分别是$x_1(r)$和$x_2(r)$的$\frac N2$点DFT
  • 对于$X(k+\frac N2)$,有:

    • 我们知道,频谱也是有周期的,以N(参与傅里叶变换的点数)个点为周期,对于DFT只是取了最中心的一个周期,而$X_1(k)$和$X_2(k)$都是$\frac N2$点DFT,故:

    • 根据旋转因子的性质有:

      • 因为$W_N^{\frac N2}=-1$
    • 所以,最终$X(k+\frac N2)$的结果为:

  • 整理上述推导,可以得到:

    • 这个式子画成信号流图就是一个蝶形变换,每求一个点的频谱,都需要知道对应偶序列的一个点和奇序列的一个点
  • 根据上式可知,只要能求出$k=0,1,…,\frac N2-1$对应的频谱,就可以同时得到$k=\frac N2,\frac N2+1,\frac N2+2,…,N-1$的频谱

    • 那么,如何得到频谱中前一半的值呢?根据公式,只需要能求出$X_1(k)$和$X_2(k)$就搞定了,而$X_1(k)$和$X_2(k)$是偶序列和奇序列的傅里叶变换,貌似也求不出来,那么就继续分解$X_1(k)$和$X_2(k)$

    • 将$X_1(k)$对应的序列又分为偶序列和奇序列,$X_2(k)$对应的序列又分为偶序列和奇序列

    • 一直这么分解下去,直到时域上只有两个点,2点DFT(这也称为蝶形变换)我们还是很容易得到的,即:

      image-20231104184905802

    • 得到这两个点的DFT,我们就得到了4个点偶序列和奇序列的DFT,进而可以求出4个点的DFT,以此类推,从而可以求出更高点数的DFT

    • 也就是说,我们对时域序列进行偶序列和奇序列分类,再对这个偶序列进行偶序列和奇序列分类,对这个奇序列进行偶序列和奇序列分类,分到最后时域上的序列两两一组,然后依次往下求蝶形变换,就可以得到最终N个序列的DFT结果,举一个8点FFT的例子,先画一个时域上两两分组的示意图:

      image-20231104191548459

    • 一个标准的8点基2-FFT流图如下:

      image-20231104191810967

  • 所以简单概括来说,基2-FFT算法是,计算一个N点DFT时,先需要计算两个$\frac N2$点DFT,然后对每个$\frac N2$点DFT,继续往下划分,还需要计算2个$\frac N4$点DFT(总共这一层就有4个$\frac N4$点DFT结果),直至最后分解到$\frac N{2^k}$($k$表示第k次分解)=2时(容易通过这个等式得到$k=log_2^{\frac N2}$),我们就能迅速计算出第$k$次分解后$2^k$个2点DFT的结果(这里单纯用数学归纳法都能计算出第$k$次分解后需要计算的DFT个数为$2^k$),之后我们就可以反推回去,得到$N$点DFT的结果

  • N点基2-FFT的计算量
    • 蝶形级数:$M=log_2^N$
    • 每级的蝶形数量:$\frac N2$
    • 每个蝶形的计算量:
      • 加法2次
      • 乘法1次
    • 总共的计算量:
      • 加法:$Nlog_2^N$
      • 乘法:$\frac N2 log_2^N$

2.基4-FFT算法

  • 设序列$x(n)$点数为$N=4^L$,$L$为正整数

  • 将$N=4^L$分为以下4组:

  • 则其$N$点DFT为:

    • 其中$X_1(k),X_2(k),X_3(k),X_4(k)$分别是$x_1(r),x_2(r),x_3(r),x_4(r)$的$\frac N4$点DFT
  • 由于$W_N^k$有以下性质:

  • 故有:

  • 与之基2类似,计算一个N点DFT时,先需要计算4个$\frac N4$点DFT,然后对每个$\frac N4$点DFT,继续往下划分,还需要计算4个$\frac{N}{16}$点DFT(总共这一层就有16个$\frac N{16}$点DFT结果),直至最后分解到$\frac N{4^k}$=4时(容易通过这个等式得到$k=log_4^{\frac N4}$),我们就能比较容易计算出第$k$次分解后$4^k$个4点DFT的结果(这里单纯用数学归纳法都能计算出第$k$次分解后需要计算的DFT个数为$4^k$),之后我们就可以反推回去,得到$N$点DFT的结果(只是这里反推时乘以的旋转因子与基2时不太一样,这里的蝶形运算是以4个点为基本单位)。

    • 以16点FFT为例,其运算过程的可视化示意图如下:(如果是64点FFT,就把下面这个图看成一个整体,复制粘贴4份分析)

    image-20231104205651274

    • 图中第一层的4点DFT的算法和基2FF算法计算4点FFT一样

      image-20231104205958813

    • 将两幅图合起来理解就是:

      image-20231104210238454

    • 一个标准的16点基4-FFT流程图如下:

      image-20231104210134227

  • N点基4-FFT的计算量

    • 蝶形级数:$M=log_4^N$
    • 每级的蝶形数量:$\frac N4$
    • 每个蝶形运算有4点输入和4点输出,则:
      • 加法8次
      • 乘法3次
    • 总共的计算量:
      • 加法:$2Nlog_4^N$
      • 乘法:$\frac {3N}4 log_4^N$

按频率抽选FFT算法

1.基2-FFT算法

  • 设序列$x(n)$点数为$N=2^L$,$L$为正整数

  • 将序列按$n$的顺序分为前后两组,其$N$点DFT为:

  • 按$k$的奇偶可将$X(k)$分为两部分

  • 则有:

    • 令:
  • 则得到:

    • 仔细观察$\sum_{n=0}^{\frac N2-1}x_1(n)W_{\frac N2}^{nr}$这玩意不就是$x_1(n)$的$\frac N2$的DFT么
  • 简而言之,有:

  • 也就是$X(k)$的偶序列只跟$x_1(n)$(与时域序列对半分有关)的傅里叶变换有关,$X(k)$的奇序列只跟$x_2(n)$(与时域序列对半分有关)的傅里叶变换有关,那么问题的关键在于首先要知道$x_1(n)$和$x_2(n)$是啥,如果能将$x_1(n)$和$x_2(n)$继续往下分,使得最后只剩下两个点,那么做2点DFT就可以得到$X(k)$中2个点的DFT

    • 以8点FFT为例,整个运算过程的可视化示意图如下:

    image-20231105133359626

    • 下面是一个按频率抽选运算的完整的8点FFT流图:

      image-20231105133621770

2.基4-FFT算法

  • 设序列$x(n)$点数为$N=4^L$,$L$为正整数

  • 序列$x(n)$按$n$的顺序分为前后4组,其N点DFT为:

  • 按$k$将$X(k)$分为4部分

  • 则有:

    • 令:
  • 则得到:

  • 简而言之,有:

  • 这样我们就把一个N点DFT按$k$的抽取分解为了4个$\frac N4$点DFT

  • 下面是一个按频率抽选运算的16点FFT流图:

    image-20231105141302313


Reference

欢迎来到ssy的世界