快速傅里叶变换和快速卷积算法

1. 快速傅里叶变换

有限长序列的傅里叶变换公式为
$$X(k)=DFT[x(n)]=\sum _{n=0}^{N-1}x(n)W_N^{kn},W_N=e^{-j\frac{2\pi }{N}}\dots \dots （1）$$
将序列x(n)按奇偶分为两个系列
$$\{\begin{array}{l}{x}_1(r)=x(2r)\\ x_2(r)=x(2r+1)\end{array},r=0,1,...,\frac{N}{2}-1$$
于是
$$X(k)=\sum _{r=0}^{N/2-1}x(2r)W_N^{k2r}+\sum _{r=0}^{N/2-1}x(2r+1)W_N^{k(2r+1)}$$$$=\sum _{r=0}^{N/2-1}{x}_1(r)W_{N/2}^{kr}+W_N^k\sum _{r=0}^{N/2-1}{x}_2(r)W_{N/2}^{kr}k=0,1,...,N-1\dots \dots （2）$$
设x₁(n)和x₂(n)的N/2点DFT分别为X₁(k)和X₂(k)，带入（2）式中有
$$X(k)=X_1(k)+W_N^k{X}_2(k),k=0,1,...,\frac{N}{2}-1\dots \dots （3）$$
为得到另一半k值的表达式，将（2）式中的k用k+N/2代替，即
$$X(k+\frac{N}{2})=\sum _{r=0}^{N/2-1}{x}_1(r)W_{N/2}^{(k+N/2)r}+W_N^{k+N/2}\sum _{r=0}^{N/2-1}{x}_2(r)W_{N/2}^{(k+N/2)r}$$$$=\sum _{r=0}^{N/2-1}{x}_1(r)W_{N/2}^{kr}-W_N^k\sum _{r=0}^{N/2-1}{x}_2(r)W_{N/2}^{kr}$$
即
$$X(k)=X_1(k)-W_N^k{X}_2(k),k=\frac{N}{2},\frac{N}{2}+1,...,N-1\dots \dots （4）$$
由（3）（4）两式可知，一个N点DFT可以分解为两个N/2点DFT计算，这两个N/2点DFT计算通过（3）（4）两式的运算组合成N点DFT。这种运算的信号流图如下，因为形似蝴蝶，故这种运算被称为蝶形运算。

若$N=2^L$，则这个分解可以一只持续下去，直到分解得到的子序列式2点序列为止，易得，两点DFT表示为如下形式
$$\{\begin{array}{l}{X}_b(0)=x_b(0)+W_2^0{x}_b(1)\\ X_b(1)=x_b(0)-W_2^0{x}_b(1)\end{array}$$
与（3）（4）两式一致，故2点DFT也可化为蝶形运算，于是整个N点DFT可以分解为L级蝶形运算，这样就大大降低了DFT的运算量。

import numpy as np

def dft1d(x):
    N=len(x)
    x=np.mat(x)
    W=np.mat(range(N)).transpose()*np.mat(range(N))
    W=np.exp(-2j*np.pi/N*W)
    return x*W
    
def binaryInverse(n,L):
    b=bin(n)[2:].zfill(L)
    b=b[-1::-1]
    return int(b,2)
fftResort=np.frompyfunc(binaryInverse,2,1)

def fft1d(x,oversample=1):
    N0=len(x)
    N=2**(len(bin(N0-1)[2:]))*2**(oversample-1)
    L=len(bin(N0-1)[2:])+oversample-1
    x+=[0]*(N-N0)
    x=np.array(x,dtype=np.complex)
    x=x[np.int16(fftResort(range(N),L))]
    for n in range(L):
        for k in range(int(2**n)):
            W=np.exp(-2j*np.pi/2**(n+1)*k)
            x[k::int(2**(n+1))],x[k+int(2**n)::int(2**(n+1))]=\
            x[k::int(2**(n+1))]+W*x[k+int(2**n)::int(2**(n+1))],\
            x[k::int(2**(n+1))]-W*x[k+int(2**n)::int(2**(n+1))]
    return x 

2. 快速卷积算法

2.1 基本步骤

快速卷积运算步骤如下：

2.2 运算量

若直接按照$y=\sum _{n=0}^{M-1}h(m)x(n-m)$计算卷积，则总的乘法次数为
$$m_d=ML\dots \dots （5）$$
快速卷积运算需要做3次FFT和N次复数乘法，计算量为
$$m_f=N(\frac{3}{2}lo{g}_{2}N+1)\dots \dots （6）$$
比值为
$$K_m=\frac{m_d}{m_f}=\frac{ML}{N(\frac{3}{2}lo{g}_{2}N+1)}=\frac{ML}{(M+L-1)(\frac{3}{2}lo{g}_2(M+L-1)+1)}$$
当$L\approx M$时，
$$K_m\approx \frac{L}{3lo{g}_{2}L+5}$$
比值随L增加而增大，说明序列越长，快速卷积方法效率越高。
但是当$L>>M$时
$$K_m\approx \frac{M}{\frac{3}{2}lo{g}_{2}L+1}$$
当L足够大时，快速卷积法反而比直接计算卷积效率低。

2.3 分段卷积

实际信号处理中，序列x(n)的长度大大超过h(n)是很普遍的，为了提高计算效率，应先将序列x(n)分段，使每段长度都与h(n)长度相当，再进行卷积，分段卷积后，还需将卷积结果拼合为长序列的卷积输出，根据分段方式的不同，合成算法也不一样，常用的有重叠相加法和重叠保留法。

2.3.1 重叠相加法

将x(n)分为若干长度为L的子序列，其中每个子序列的表达式为
$$x_i(n)=\{\begin{array}{ll}x(n+iL)& n=0,1,...,L-1\\ 0& 其它\end{array}\dots \dots （7）$$
$$y(n)=x(n)\ast h(n)=[\sum _{i=0}^{\infty }{x}_i(n-iL)]\ast h(n)$$$$=\sum _{i=0}^{\infty }{x}_i(n-iL)\ast h(n)$$
由于LSI系统的时不变特性，有
$$y_i(n-iL)=x_i(n-iL)\ast h(n)$$
故
$$y(n)=\sum _{i=0}^{\infty }{y}_i(n-iL)$$
在按照（7）式将卷积结果叠加过程中，y_i(n)的长度是L+M-1，而每一段的位移是iL，前一段y_i-1(n)值的后M-1个样会与后一段y_i(n)的前M-1个样值重叠，故这种算法称为重叠相加法。

2.3.2 重叠保留法

将y(n)均匀分为长度为L没有重叠的子序列，即令
$$y_i(n)=\{\begin{array}{ll}y(n+iL)& n=0,1,...L-1\\ 0& 其它\end{array}$$
由图可见，y_i(n-iL)不仅与第i段输入x_i(m-iL)有关，还与前一段的后M-1个值有关，故对x(n)分段时，每一段x_i(n)要包含前一段的后M-1个值，分为长度为N=L+M-1相互之间有重叠的子序列x_i(n)，即
$$x_i(n)=\{\begin{array}{ll}x(n+iL-M+1)& n=0,1,...,N-1\\ 0& 其它\end{array}\dots \dots （8）$$
设x_i(n)与h(n)的卷积结果为y_i^’(n)，有由于x_i(n)前面增加了M-1个样值，y_i^’(n)的起始样值与输出分段y_i(n)的起始样值不一样，y_i(n)应是y_i^’(n)的第M-1号样值，故
$$y_i(n)=y_i^{{}^{\prime }}(n+M-1),n=0,1,...,L-1$$