Python实现DIT-FFT算法：递归与迭代两种方法对比
发布时间：2026-10-16 |
3人浏览过
本文详细介绍基于Python实现基2时域抽取快速傅里叶变换（DIT-FFT）的两种核心方法：递归实现与迭代流式计算。通过对比分析，指出递归法虽直观但内存开销大，而迭代法通过位反转和复用数组显著节省内存。文章包含完整代码示例及运行结果分析，帮助开发者根据实际需求选择最优实现方案。
手动编写函数实现FFT
采用递归策略进行计算
from math import log, ceil, cos, sin, pi
import matplotlib.pyplot as plt
import numpy as np
# 这两行代码解决 plt 中文显示的问题
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
def fft(x, N=None):
# DIT-FFT 函数说明
# x: 时域序列
# N: N点DFT, 理论上N=2**M
# 返回值为序列x的DFT
if N is None:
N = len(x)
elif N < len(x):
N = len(x)
if N == 2:
return [x[0]+x[1], x[0]-x[1]]
# 补0使得N=2**M
M = ceil(log(N, 2))
N = 2**M
x = x + [0] * (N-len(x))
# 递归地计算偶数项和奇数项的DFT
X1 = fft(x[0::2])
X2 = fft(x[1::2])
X = [0] * N
for i in range(N//2):
# 蝶形计算
tmp = (cos(2*pi/N*i)-1j*sin(2*pi/N*i))*X2[i]
X[i] = X1[i] + tmp
X[i+N//2] = X1[i] - tmp
return X
if __name__ == '__main__':
x = [1]*10
y = fft(x, 1024)
# print(y)
z = [abs(i) for i in y]
# print(z)
plt.plot(np.arange(len(z))*2/len(z), z, label='10点矩形窗函数的FFT')
plt.title("幅度谱")
plt.xlabel(r'单位：$\pi$')
plt.ylabel(r'$|H(j\omega)|$')
plt.grid(linestyle="-.")
plt.legend()
plt.show()
改用循环结构，通过流式计算大幅降低内存占用
from math import log, ceil, cos, sin, pi
import matplotlib.pyplot as plt
import numpy as np
# 这两行代码解决 plt 中文显示的问题
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
def fft(x, N=None):
# DIT-FFT 函数说明
# x: 时域序列
# N: N点DFT, 理论上N=2**M
# 返回值为序列x的DFT
"""
采用流式计算方法，只用了一个N(N=2**M)点的数组内存
"""
if N is None:
N = len(x)
elif N < len(x):
N = len(x)
# 补0使得：N=2**M
M = ceil(log(N, 2))
N = 2**M
x = x + [0] * (N-len(x))
fm = "{:0"+f"{M}"+"b}"
X = [0] * N
for i in range(N//2):
index1 = eval('0b'+fm.format(i*2)[::-1])
index2 = eval('0b'+fm.format(i*2+1)[::-1])
X[2*i] = x[index1] + x[index2]
X[2*i+1] = x[index1] - x[index2]
for i in range(1, M):
# 第i步表示将2**i点DFT合成2**(i+1)点的DFT
# 蝶形宽度width
width = 2**i
"""
将X(k)序列进行分组，每组2**(i+1)个点，
便于将每组中两组2**i点DFT合成一组2**(i+1)点的DFT
"""
# num=2*width=2**(i+1), 表示每组点数
num = 2*width
# 组数groups
groups = N//num
for j in range(groups):
# 对每组将2**i点DFT合成2**(i+1)=num点的DFT
for k in range(num//2):
# 旋转因子
W = cos(2*pi/num*k) - 1j * sin(2*pi/num*k)
# 第j组第k个
index = j*num + k
tmp = W * X[index+width] # 每个蝶形一次复数乘法
X[index], X[index+width] = X[index]+tmp, X[index]-tmp
return X
if __name__ == '__main__':
x = [1]*10
y = fft(x, 1024)
# print(y)
z = [abs(i) for i in y]
# print(z)
plt.plot(np.arange(len(z))*2/len(z), z, label='10点矩形窗函数的FFT')
plt.title("幅度谱")
plt.xlabel(r'单位：$\pi$')
plt.ylabel(r'$|H(j\omega)|$')
plt.grid(linestyle="-.")
plt.legend()
plt.show()程序执行输出如下：
# 建议优先采用第二种方法实现FFT。第一种递归方式在调用过程中会产生额外的栈成本，且内存消耗较大；相比之下，第二种方法仅依赖一个大小为N（其中N=2**M）的数组进行运算，实现了内存的有效复用。
利用Python第三方库执行FFT
import numpy as np
from numpy.fft import fft, ifft
# from scipy.fftpack import fft, ifft
import matplotlib.pyplot as plt
# 这两行代码解决 plt 中文显示的问题
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
if __name__ == '__main__':
x = 2*np.sin(np.pi/2*np.arange(100))+np.sin(np.pi/5*np.arange(100))
z = [abs(i) for i in fft(x, 2048)]
# print(z)
L = len(z)
plt.plot((np.arange(L)*2/L)[:L//2], z[:L//2], label='两个不同频率正弦信号相加的DFT')
plt.title("幅度谱")
plt.xlabel('$\pi$')
plt.ylabel('$|H(j\omega)|$')
plt.grid(linestyle="-.")
plt.legend()
plt.show()
print('max(abs(ifft(fft(x))-x)) = ', end='')
print(max(abs(ifft(fft(x))-x)))库函数调用结果如下：
max(abs(ifft(fft(x))-x)) = 9.01467522361575e-16以上即为基于Python实现DIT-FFT算法的完整解析，更多关于Python DIT-FFT算法的技术资料，敬请参阅本站其他相关文章。
本站声明
本文内容仅供参考，不代表本站观点。转载或引用请注明出处。如发现内容问题，请通过联系我们页面反馈。
