快速傅里叶变换(FFT)

前言

傅里叶级数(FS)

傅里叶变换(FT)

离散时间傅里叶级数(DFS)

离散时间傅里叶变换(DTFT)

离散傅里叶变换(DFT)

建议先看以上文章

FFT是DFT的一种快速算法而不是一种新的变换,它可以在数量级的意义上提高运算速度。

直接计算DFT的问题

DFT的运算量

设有限长序列x(n),非零值长度为N,若对x(n)进行次DFT运算,共需多大的运算工作量?

X(k)=DFT[x(n)]=\sum_{n=0}^{N-1}x(n)W_{N}^{nk},k=0,1,\cdots N-1 
x(n)=IDFT[x(n)]=\frac{1}{N}\sum_{n=0}^{N-1}x(n)W_{N}^{-nk},n=0,1,\cdots N-1
  • x(n)为复数,W_{N}^{nk}=e^{-jk\frac{2\pi}{N} n}也为复数。
  • X(k)和x(n)的计算量基本相同,只差一个因子\frac{1}{N}

复数乘法次数N^{2},复数加法次数N(N-1),若N远远大于1,则这两者都近似为N^{2},随N增大而急速增大。

改善DFT运算效率的基本途径

1、利用W_{N}^{nk}=e^{-jk\frac{2\pi}{N} n}的性质

W_{N}^{nk}(单位周期复指数序列、旋转因子)具有以下性质

  • 周期性 W_{N}^{nk}=W_{N}^{(n+rN)k}=W_{N}^{(k+rN)n},r为整数

        W_{N}^{(n+rN)k}= W_{N}^{nk}\cdot W_{N}^{rNk}=W_{N}^{nk}\cdot 1=W_{N}^{nk}

        同理可得W_{N}^{(k+rN)n}=W_{N}^{nk}

  • 共轭对称性 \left ( W_{N}^{nk} \right ) ^{*}=W_{N}^{-nk}=W_{N}^{(rN-n)k}=W_{N}^{(rN-k)n},r为整数
  • 可约性 W_{N}^{nk}=W_{mN}^{mnk},W_{N}^{nk}=W_{\frac{N}{m}}^{\frac{nk}{m}},m为整数,N/m为整数

        W_{mN}^{mnk}=e^{-jk\frac{2\pi}{Nm} mn}=e^{-jk\frac{2\pi}{N} n}=W_{N}^{nk}

        同理可得W_{\frac{N}{m}}^{\frac{nk}{m}}=W_{N}^{nk}

  • 由以上性质,得出一些特殊点 

        W_{N}^{0}=e^{0}=1

        W_{N}^{N}=e^{-2\pi j}=1

        W_{N}^{\frac{N}{2}}=e^{-j\pi}=-1

        W_{N}^{k+\frac{N}{2}}= W_{N}^{\frac{N}{2}}\cdot W_{N}^{k}=-W_{N}^{k}

FFT算法的基本思想与分类

基本思想

利用DFT系数的特性,合并DFT运算中的某些项;

把长序列DFT分解为短序列DFT,从而减少运算量。

FFT算法分类

  • 时间抽选法(DIT):Decimation-In-Time
  • 频率抽选法(DIF):Decimation-In-Frequency

按时间抽选(DIT)的基-2 FFT算法 

算法原理

设输入序列长度为N=2^{M}(M为正整数),将该序列按时间顺序的奇偶分解为越来越短的子序列,称为基2按时间抽取的FFT算法。也称为Coolkey-Tukey(库利-图基)算法。

其中基2表示:N=2^{M},M为整数,若不满足这个条件,可以人为地加上若干零值(加零补长)使其达到N=2^{M}

N=2^{M}的序列x(n),按奇偶分成两组:

\left.\begin{matrix} x(2r) = x_{1}(r) \\ x(2r+1) = x_{2}(r) \end{matrix}\right\},r=0,1,\cdots ,\frac{N}{2}-1                                                                          (1)

则其N点DFT为

X(k)=DFT[x(n)]=\sum_{n=0}^{N-1}x(n)W_{N}^{nk}

           =\sum_{r=0}^{\frac{N}{2}-1} x(2r)W_{N}^{2rk}+\sum_{r=0}^{\frac{N}{2}-1} x(2r+1)W_{N}^{(2r+1)k}

           =\sum_{r=0}^{\frac{N}{2}-1} x_{1}(r)W_{N}^{2rk}+W_{N}^{k}\sum_{r=0}^{\frac{N}{2}-1} x_{2}(r)W_{N}^{2rk}

利用W_{N}^{nk}的可约性

X(k)=\sum_{r=0}^{\frac{N}{2}-1} x_{1}(r)W_{\frac{N}{2}}^{rk}+W_{N}^{k} \sum_{r=0}^{\frac{N}{2}-1} x_{2}(r)W_{\frac{N}{2}}^{rk}=X_{1}(k) + W_{N}^{k} X_{2}(k)                      (2)

式中X_{1}(k)X_{2}(k)分别是x_{1}(r)x_{2}(r)\frac{N}{2}点DFT

X_{1}(k)=\sum_{r=0}^{\frac{N}{2}-1} x_{1}(r)W_{\frac{N}{2}}^{rk}=\sum_{r=0}^{\frac{N}{2}-1} x(2r)W_{\frac{N}{2}}^{rk}                                                                 (3)

X_{2}(k)=\sum_{r=0}^{\frac{N}{2}-1} x_{2}(r)W_{\frac{N}{2}}^{rk}=\sum_{r=0}^{\frac{N}{2}-1} x(2r+1)W_{\frac{N}{2}}^{rk}                                                         (4)

由(2)式可看出,一个N点DFT已分解成两个\frac{N}{2}点DFT,但是x_{1}(r)x_{2}(r)以及X_{1}(k)X_{2}(k)都是 \frac{N}{2}点序列,即r,k 满足r,k=0,1,\cdots ,\frac{N}{2}-1。而X(k)却有 N点,用(2)式计算得到的只是 X(k)的前一半项数的结果,要用X_{1}(k)X_{2}(k)表达全部的X(k)值,还必须应用W_{N}^{nk}的周期性,即

W_{\frac{N}{2}}^{rk}= W_{\frac{N}{2}}^{r(k+\frac{N}{2})}

这样可得到

X_{1}(k+\frac{N}{2})=\sum_{r=0}^{\frac{N}{2}-1} x_{1}(r)W_{\frac{N}{2}}^{r(k+\frac{N}{2})}=\sum_{r=0}^{\frac{N}{2}-1} x_{1}(r)W_{\frac{N}{2}}^{rk}=X_{1}(k)                               (5)

同理可得

X_{2}(k+\frac{N}{2})=X_{2}(k)                                                                                                                 (6)

(5)式、(6)式说明了后半部分k值 \left ( \frac{N}{2}\leqslant k\leq N-1 \right ) 所对应的 X_{1}(k)X_{2}(k)等于前半部分k值\left ( 0\leqslant k\leq \frac{N}{2}-1 \right )所对应的 X_{1}(k)X_{2}(k)

再利用特殊点

W_{N}^{k+\frac{N}{2}}= W_{N}^{\frac{N}{2}}\cdot W_{N}^{k}=-W_{N}^{k}                                                                                                 (7)

整理得

前半部分 X(k) \left ( k=0,1\cdots ,\frac{N}{2}-1 \right )可表示为

X(k)=X_{1}(k) + W_{N}^{k} X_{2}(k),    k=0,1,\cdots ,\frac{N}{2}-1                                                          (8)

后半部分 X(k) \left ( k=\frac{N}{2},\cdots ,N-1 \right )可表示为

X(k+\frac{N}{2})=X_{1}(k) - W_{N}^{k} X_{2}(k),    k=0,1,\cdots ,\frac{N}{2}-1                                                 (9)

这样,只要求出0\sim \left ( \frac{N}{2} -1\right )区间的所有X_{1}(k)X_{2}(k)值,即可求出0~(N-1)区间内的所有 X(k)值,大大节省了运算。

(8)式和(9)式的运算可以用下图所示的蝶形运算流图符号表示。

N=2^{3}=8为例子

继续分解,把\frac{N}{2}点DFT,分解成\frac{N}{4}点DFT

先将x_{1}(r)进行分解:

\left.\begin{matrix} x_{1}(2l) = x_{3}(l) \\ x_{1}(2l+1) = x_{4}() \end{matrix}\right\},l=0,1,\cdots ,\frac{N}{4}-1                                                                        (10)

同样可得

X_{1}(k)=\sum_{l=0}^{\frac{N}{4}-1} x_{1}(2l)W_{\frac{N}{2}}^{2lk}+\sum_{l=0}^{\frac{N}{4}-1} x_{1}(2l+1)W_{\frac{N}{2}}^{(2l+1)k}

            =\sum_{l=0}^{\frac{N}{4}-1} x_{3}(l)W_{\frac{N}{2}}^{2lk}+W_{\frac{N}{2}}^{k}\sum_{l=0}^{\frac{N}{4}-1} x_{4}(l)W_{\frac{N}{2}}^{2lk}

            =\sum_{l=0}^{\frac{N}{4}-1} x_{3}(l)W_{\frac{N}{4}}^{lk}+W_{\frac{N}{2}}^{k}\sum_{l=0}^{\frac{N}{4}-1} x_{4}(l)W_{\frac{N}{4}}^{lk}

           =X_{3}(k) + W_{\frac{N}{2}}^{k} X_{4}(k),    k=0,1,\cdots ,\frac{N}{4}-1                                                       

X_{1}(k+\frac{N}{4})=X_{3}(k) - W_{\frac{N}{2}}^{k} X_{4}(k),    k=0,1,\cdots ,\frac{N}{4}-1                                             

其中

X_{3}(k)=\sum_{l=0}^{\frac{N}{4}-1} x_{1}(2l)W_{\frac{N}{4}}^{lk}=\sum_{l=0}^{\frac{N}{4}-1} x_{3}(l)W_{\frac{N}{4}}^{lk}                                                                (11)

X_{4}(k)=\sum_{l=0}^{\frac{N}{4}-1} x_{1}(2l+1)W_{\frac{N}{4}}^{lk}=\sum_{l=0}^{\frac{N}{4}-1} x_{4}(l)W_{\frac{N}{4}}^{lk}                                                        (12)

同理x_{2}(r)也可以分解得到

\left.\begin{matrix} X_{2}(k+)=X_{5}(k) + W_{\frac{N}{2}}^{k} X_{6}(k) \\ X_{2}(k+\frac{N}{4})=X_{5}(k) - W_{\frac{N}{2}}^{k} X_{6}(k)\end{matrix}\right\},k=0,1,\cdots ,\frac{N}{4}-1

其中

X_{5}(k)=\sum_{l=0}^{\frac{N}{4}-1} x_{2}(2l)W_{\frac{N}{4}}^{lk}=\sum_{l=0}^{\frac{N}{4}-1} x_{5}(l)W_{\frac{N}{4}}^{lk}                                                                (13)

X_{6}(k)=\sum_{l=0}^{\frac{N}{4}-1} x_{2}(2l+1)W_{\frac{N}{4}}^{lk}=\sum_{l=0}^{\frac{N}{4}-1} x_{6}(l)W_{\frac{N}{4}}^{lk}                                                        (14)

通常会把W_{\frac{N}{2}}^{k} 写成 W_{N}^{2k},这样,一个N点DFT就可以分解成4个\frac{N}{4}点DFT

N=2^{3}=8为例子,进行两次时域抽取分解

如果上图的2点DFT也能化成蝶形运算,就很完美了

X_{3}(k)=\sum_{l=0}^{\frac{N}{4}-1} x_{1}(2l)W_{\frac{N}{4}}^{lk}=\sum_{l=0}^{\frac{N}{4}-1} x_{3}(l)W_{\frac{N}{4}}^{lk},    k=0,1

展开X_{3}(k)有:

X_{3}(0)=x_{3}(0)+W_{2}^{0}x_{3}(1)=x(0)+W_{2}^{0}x(4)=x(0)+W_{8}^{0}x(4)

X_{3}(1)=x_{3}(0)+W_{2}^{1}x_{3}(1)=x(0)+W_{2}^{1}x(4)=x(0)-x(4)

因为W_{8}^{0}=1,所以X_{3}(1)=x(0)-W_{8}^{0}x(4)

由此可以看出2点DFT原本就可以看成一个蝶形运算

对比直接DFT和FFT的运算耗时

编程方法

未完待续。。。

编程实例

未完待续。。。

参考资料

《数字信号处理 华东理工大学 万永菁》课件、视频

《数字信号处理教程》(第5版)程佩青

  • 16
    点赞
  • 15
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值