快速傅里叶变换

1970年01月01日算法

快速傅里叶变换及逆变换

快速傅里叶变换是在离散傅里叶的基础上的优化型

本质是一样的

首先我们来看看离散傅里叶变换的公式：
$D(k) = \sum_{n=0}^{N-1}E(n)(\cos\frac{2\pi kn}{N} - j\sin\frac{2\pi kn}{N})$

方便期间我们设
$W^{n}_{N} = \cos\frac{2\pi n}{N} + j\sin\frac{2\pi n}{N}=e^{j\frac{2\pi n}{N}}$
即欧拉公式的变体

$W^{n}{N}$有如下特性
$W^{a}{N} * W^{b}{N} = W^{a+b}{N}$
$W^{a + \frac{N}{2}}{N} = \cos{(\frac{2\pi n}{N} + \pi)} + j\sin{(\frac{2\pi n}{N} + \pi)} = -W^{a}{N}$
$W^{a + N}{N} = \cos{(\frac{2\pi n}{N} + 2\pi)} + j\sin{(\frac{2\pi n}{N} + 2\pi)} = W^{a}{N}$
该特性也是快速傅里叶的关键

则
$D(k) = \sum_{n=0}^{N-1}E(n)W^{-kn}{N}$
$D(k) = E(0)W^{0}{N} + E(1)W^{-k}{N} + E(2)W^{-2k}{N} + E(3)W^{-3k}{N} + … + E(N-1)W^{-(N-1)k}{N}$

可以得到：(N必须为2的次方,则n最大为N-1)
$D(k) = (E(0) + E(2)W^{-2k}{N} + … + E(N-2)W^{-(N-2)k}{N}) + (E(1)W^{-1k}{N} + E(3)W^{-3k}{N} + … + E(N-1)W^{-(N-1)k}{N})$
$D(k) = (E(0) + E(2)W^{-2k}{N} + … + E(N-2)W^{-(N-2)k}{N}) + W^{-1k}{N}(E(1) + E(3)W^{-2k}{N} + … + E(N-1)W^{-(N-2)k}{N})$

设
$D1(k) = E(0) + E(2)W^{-2k}{N} + … + E(N-2)W^{-(N-2)k}{N}$
$D2(k) = E(1) + E(3)W^{-2k}{N} + … + E(N-1)W^{-(N-2)k}{N}$

则我们可以得到D(k)的递推式
$D(k) = D1(k) + W^{-1k}{N}D2(k)$
$D(k+\frac{N}{2}) = D1(k) - W^{-1k}{N}D2(k)$
$0<=k<=\frac{N}{2}-1$

当N=1时
$D(k) = E(0)$

由其递归式我们就可以写出快速傅里叶递归代码

def FFT(x):
    N = len(x)
    if N <= 1:
        return x
    even = FFT(x[0::2])#取出偶数项
    odd = FFT(x[1::2])#取出奇数项
    #计算W*D2
    T = [np.exp(-2j * math.pi * k / N) * odd[k] for k in range(N//2)]
    #合并前面一半和后面一半
    return [even[k] + T[k] for k in range(N//2)] + [even[k] - T[k] for k in range(N//2)]

由离散傅里叶逆变换公式
$ID(k) = \frac{1}{N}\sum_{n=0}^{N-1}D(n)(cos\frac{2\pi kn}{N} + jsin\frac{2\pi kn}{N})$

同理可得快速傅里叶的逆变换
推导过程与上面相似
与正变换的唯一区别在于
$W^{-nk}_{N}$中负号去掉了


def IFFT(x, n=0):
    N = len(x)
    if N <= 1:
        return x
    even = IFFT(x[0::2], n+1)
    odd = IFFT(x[1::2], n+1)
    #这里与傅里叶变换的唯一区别就在于符号变了
    T = [np.exp(2j * math.pi * k / N) * odd[k] for k in range(N//2)]
    if n == 0:
        return [(even[k] + T[k])/N for k in range(N//2)] + [(even[k] - T[k])/N for k in range(N//2)]
    return [(even[k] + T[k]) for k in range(N//2)] + [even[k] - T[k] for k in range(N//2)]