本文参考视频见此处,此视频从更熟悉的问题入手介绍FFT思想,非常容易看的进去,但是由于该视频讲解节奏较快,中间有些步骤对我来说有些飘渺,所以借用python代码对所讲算法进行实现,并补充视频中未明确提到的边界条件和递推公式等等,根据代码结果观察效果。此博客记录了实现过程的代码,如有错误欢迎补充。
(补充:看到后面发现他有讲伪代码,实现逻辑确实差不多)
思路:将DFT问题看成一个多项式评估问题
一句话概括:从DFT的变换式来看,其实求
X
k
X_k
Xk就是求多项式函数
f
(
y
)
=
∑
n
=
0
N
−
1
x
[
n
]
⋅
y
n
f(y)=\sum_{n=0}^{N-1} x[n] \cdot y^n
f(y)=∑n=0N−1x[n]⋅yn 在
y
=
e
−
j
2
π
k
N
y=e^{\frac{-j 2\pi k} {N}}
y=eN−j2πk处的函数值。
详细分析:
- 多项式函数一般表达为 g ( x ) = a 0 ⋅ x 1 + a 1 ⋅ x 1 + . . . + a N − 1 ⋅ x N − 1 = ∑ n = 0 N − 1 a n ⋅ x n 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 g(x)=a0⋅x1+a1⋅x1+...+aN−1⋅xN−1=∑n=0N−1an⋅xn与上文中的f(y)形式完全一致。
- 在复平面上看, y = e − j 2 π k N , k = 0 , 1 , 2 , . . . , N − 1 y=e^{\frac{-j 2\pi k} {N}}, k=0,1,2,...,N-1 y=eN−j2πk,k=0,1,2,...,N−1这N个点,恰恰是平均分割单位圆的N个点。
- 所谓多项式评估式表达法,是对应于解析式表达法的。类似于直线的两点式,给定N个点,便可以确定一个N-1维的多项式。对于DFT来说,我们是已知解析式的,现在想要求评估式表达法。最朴素的算法就是一个一个顺序求和相加,复杂度为O(N)。但是由于这N个点已经指定了,并且非常的特殊(如以上2.所述),所以是存在优化空间的,经过优化的算法可以达到O(logN)的复杂度。(具体优化方法见原视频)
程序1:分治法解决多项式评估(仅一次分治)
给定一多项式P(x)的系数向量,求一个 x i , f ( x i ) {x_i,f(x_i)} xi,f(xi)对。
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 ) {x_i,f(x_i)} xi,f(xi)对
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]时域信号即可。