***************************************************************
# 与scipy库里的fft作比较
***************************************************************
import numpy as np
import math
from scipy.fftpack import fft
def my_fft(x, N):
"""
实现FFT
x: 离散序列
N: x的长度
"""
# 反序排列
F = np.array(x)
LH = int(N / 2)
J = LH
# print(J)
M = N - 2
for I in range(1, M+1):
if I < J:
T = F[I]
F[I] = F[J]
F[J] = T
K = LH
while J >= K:
J = int(J - K)
K = int(K / 2)
J = J + K
# 蝶形运算
ls = [i for i in range(int(N/2))]
a = np.array(ls)
Wmn = np.exp(-1j * 2 * math.pi * a / N)
Wmn = list(Wmn)
F = list(F)
n1 = 1
n2 = 1
group = N
for k in range(1, K+1):
n2 = n1
n1 = n2 * 2
group = int(group / 2)
for j in range(1, group+1):
h = (j-1) * n1 + 1
for m in range(h, h+n2):
W1 = Wmn[int((m-h) * N / n1)]
W = W1 * F[m + n2 - 1]
F[m + n2 - 1] = F[m - 1] - W
F[m - 1] = complex(F[m - 1]) + W
return F
test_series = [3, 2, 6, 8, 10, 22, 1, 2]
N = len(test_series)
F1 = my_fft(test_series, N) # 我的fft
print("my_fft:",np.array(F1))
print("------------------------------------------------------------")
F2 = fft(test_series) # scipy的fft
print("scipy_fft:", F2)
"""
my_fft: [ 54. +0.j -25.38477631 +4.89949494j
6. -14.j 11.38477631+14.89949494j
-14. +0.j 11.38477631-14.89949494j
6. +14.j -25.38477631 -4.89949494j]
------------------------------------------------------------
scipy_fft: [ 54. -0.j -25.38477631 +4.89949494j
6. -14.j 11.38477631+14.89949494j
-14. -0.j 11.38477631-14.89949494j
6. +14.j -25.38477631 -4.89949494j]
"""
基2的FFT算法 | Python实现
最新推荐文章于 2022-10-28 09:55:25 发布