Radix-2 FFT 网上的资源很多,但是Radix-4 FFT的资源很少,我只找到一个C++版本的,而且网上几乎没有 IFFT 的代码。
我实现了一个python版本的Radix-4,包括正变换,和反变换,方便理解和学习。
其中第一个版本的复杂度更低,第二个版本更方便理解。
建议先从第二个版本看起。
最好是先理解Radix-2 FFT再学习Radix-4的版本,Radix-2 FFT是普通的碟式算法FFT,网上很多资料。
参考资料:
首先,定义位翻转函数:
def bit_reverse(X, Y, N):
j = 0
N2 = N//4
for i in range(N-1):
if(i<j):
X[i],X[j] = X[j],X[i]
Y[i],Y[j] = Y[j],Y[i]
N1 = N2
while j>= 3 * N1:
j = j - 3 * N1
N1 = N1 // 4
j = j + N1
return X, Y
第一个版本的 Radix-4 FFT
def FFT_Radix4(X, Y, N, M):
X_ = np.copy(X)
Y_ = np.copy(Y)
N2 = N
for k in range(M):
E = 2 * np.pi / N2
N1 = N2
N2 = N2//4
A = 0
for j in range(N2):
A1 = j * E
A2 = 2 * A1
A3 = 3 * A1
CO1, CO2, CO3 = np.cos(A1), np.cos(A2), np.cos(A3)
SI1, SI2, SI3 = np.sin(A1), np.sin(A2), np.sin(A3)
for i in range(j,N,N1):
i0 = i + 0 * N2
i1 = i + 1 * N2
i2 = i + 2 * N2
i3 = i + 3 * N2
R1 = X_[i0] + X_[i2]
R2 = X_[i1] + X_[i3]
R3 = X_[i0] - X_[i2]
R4 = X_[i1] - X_[i3]
S1 = Y_[i0] + Y_[i2]
S2 = Y_[i1] + Y_[i3]
S3 = Y_[i0] - Y_[i2]
S4 = Y_[i1] - Y_[i3]
X_[i0] = R1 + R2;
R2 = R1 - R2;
R1 = R3 - S4;
R3 = R3 + S4;
Y_[i0] = S1 + S2;
S2 = S1 - S2;
S1 = S3 + R4;
S3 = S3 - R4;
X_[i1] = CO1*R3 + SI1*S3<