用代码尝试入门FFT

本文参考视频见此处，此视频从更熟悉的问题入手介绍FFT思想，非常容易看的进去，但是由于该视频讲解节奏较快，中间有些步骤对我来说有些飘渺，所以借用python代码对所讲算法进行实现，并补充视频中未明确提到的边界条件和递推公式等等，根据代码结果观察效果。此博客记录了实现过程的代码，如有错误欢迎补充。
(补充：看到后面发现他有讲伪代码，实现逻辑确实差不多）

思路：将DFT问题看成一个多项式评估问题

一句话概括：从DFT的变换式来看，其实求$X_k$就是求多项式函数$f(y)={\sum }_{n=0}^{N-1}x[n]\cdot y^n$ 在$y=e^{\frac{-j2\pi k}{N}}$处的函数值。
详细分析：

1. 多项式函数一般表达为$g(x)=a_0\cdot x^1+a_1\cdot x^1+...+a_{N-1}\cdot x^{N-1}={\sum }_{n=0}^{N-1}{a}_n\cdot x^n$与上文中的f(y)形式完全一致。
2. 在复平面上看，$y=e^{\frac{-j2\pi k}{N}},k=0,1,2,...,N-1$这N个点，恰恰是平均分割单位圆的N个点。
3. 所谓多项式评估式表达法，是对应于解析式表达法的。类似于直线的两点式，给定N个点，便可以确定一个N-1维的多项式。对于DFT来说，我们是已知解析式的，现在想要求评估式表达法。最朴素的算法就是一个一个顺序求和相加，复杂度为O(N)。但是由于这N个点已经指定了，并且非常的特殊（如以上2.所述），所以是存在优化空间的，经过优化的算法可以达到O(logN)的复杂度。（具体优化方法见原视频）

程序1：分治法解决多项式评估（仅一次分治）

给定一多项式P(x)的系数向量，求一个$x_i,f(x_i)$对。

import numpy as np
P_in = np.random.random(2**5).astype(np.complex128) # 分别为0, 1, ..., n-1, n项系数
cnt = 0 # 程序时间复杂度计数

# 硬算，输入一个x应该返回一个值P(x)
def calc_hard(P, xi):
    # 程序时间复杂度计数
    global cnt
    cnt += len(P) # 加权求和运算复杂度为O(len(P))
    # 返回硬计算的结果
    return np.sum([P[i]*(xi**i) for i in range(len(P))])

# 对优化后的分治算法来说，输入一个x应该返回两个值，一个是P(x)，另一个是P(-x)
def calc_part(P, x):
    global cnt
    # P(x) = Po(x^2)+x*Pe(x^2)
    Po = P[1::2] # 取奇数项系数
    Pe = P[::2] # 取偶数项系数
    # 分别硬算
    vo = calc_hard(Po, x**2)
    ve = calc_hard(Pe, x**2)
    # 返回P(x)和P(-x)
    cnt += 2 # 仅增加常量级的两次计算，也可以忽略不计
    return x*vo+ve, ve-x*vo
    
i = 1
cnt = 0
print('hard calculation: value =', (calc_hard(P_in, i), calc_hard(P_in, -i)), ', cnt =', cnt)

cnt = 0
print('parted calculatio: value =', (calc_part(P_in, i)), ', cnt =', cnt)

print('可以看到结算出来的P(x)和P(-x)值是一样的，但是分治法仅用了一半的计算量')

程序2：分治法解决多项式评估（递归）

给定一n-1次多项式f(x)的系数向量，求n个$x_i,f(x_i)$对

import numpy as np
import matplotlib.pyplot as plt
from timeit import timeit


########################### 函数功能定义 ##########################

# 对分治算法来说，应该返回n个值，R_i代表P(x_i), x_i=np.exp((1j*2*np.pi/n)*i)
def calc_part(P):
    if len(P) == 1: return np.array([P[0]])
    # P(x) = Po(x^2)+x*Pe(x^2)
    Po = P[1::2] # 取奇数项系数
    Pe = P[::2] # 取偶数项系数
    # 分治
    list_o = calc_part(Po)
    list_e = calc_part(Pe)

    n = len(P)
    # 先计算单位圆上半部分的点
    l_1 = [ve+vo*np.exp((1j*2*np.pi/n)*i) for i, (ve, vo) in enumerate(zip(list_e, list_o))]
    # 再计算单位圆下半部分的点
    l_2 = [ve-vo*np.exp((1j*2*np.pi/n)*i) for i, (ve, vo) in enumerate(zip(list_e, list_o))]
    # 合并起来返回
    return np.concatenate([l_1, l_2], 0)

# 硬算，输入一个x应该返回一个值P(x)
def calc_hard(P, xi):
    # 返回硬计算的结果
    return np.sum([P[i]*(xi**i) for i in range(len(P))])

# 输入输出格式同calc_part
def calc_hard_n_complex(P):
    return np.array([calc_hard(P_in, np.exp((1j*2*np.pi/len(P_in))*i)) for i in range(len(P_in))])


########################### 结果正确性验证 ##########################
P_in = np.random.random(2**5).astype(np.complex128) # 分别为0, 1, ..., n-1, n项系数

res = calc_part(P_in)
plt.scatter(np.real(res), np.imag(res), c='b')
plt.pause(0.5)

res = calc_hard_n_complex(P_in)
plt.scatter(np.real(res), np.imag(res), c='r')
plt.pause(0.5)


########################## 程序耗时测试 ##########################
P_in = np.random.random(2**11).astype(np.complex128)
hard_time = timeit('calc_hard_n_complex(P_in)', 'from __main__ import calc_hard_n_complex, P_in', number=1)
print('hard_time =', hard_time, 's')
part_time = timeit('calc_part(P_in)', 'from __main__ import calc_part, P_in', number=1)
print('part_time =', part_time, 's')

print('可以看到结算出来的P(x)和P(-x)值是一样的，但是分治法仅用不到十分之一的运行时间，若进一步加大n，则差距会更大')

如何将多项式评估程序修改为FFT程序？

将系数替换为x[n]时域信号即可。