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;
Y_[i1] = CO1*S3 - SI1*R3;
X_[i2] = CO2*R2 + SI2*S2;
Y_[i2] = CO2*S2 - SI2*R2;
X_[i3] = CO3*R1 + SI3*S1;
Y_[i3] = CO3*S1 - SI3*R1;
X_, Y_ = bit_reverse(X_, Y_, N)
return X_, Y_
写反变换函数,这里主要修改了两个地方,一个是 E = -2*pi /N2加了个负号,第二就是将R3,S3与R1,S1调换。
def IFFT_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*R1 + SI1*S1;
Y_[i1] = CO1*S1 - SI1*R1;
X_[i2] = CO2*R2 + SI2*S2;
Y_[i2] = CO2*S2 - SI2*R2;
X_[i3] = CO3*R3 + SI3*S3;
Y_[i3] = CO3*S3 - SI3*R3;
X_, Y_ = bit_reverse(X_, Y_, N)
return X_, Y_
另外,这里根据 反变换的共轭求法,给出一个共轭版本。
def IFFT_Radix4_conjugate(X, Y, N, M):
Y = - Y
res_X, res_Y = FFT_Radix4(X, Y, N, M)
return res_X, -res_Y
接下来是实验,可以看到,经过正反变换之后,得到的结果和 numpy自带函数一致。
import numpy as np
M = 2
N = 4 ** M
X = np.arange(N)
Y = np.arange(N, 2*N)
X = X.astype(np.float)
Y = Y.astype(np.float)
X_fft, Y_fft = FFT_Radix4(X, Y, N, M)
X_ifft, Y_ifft = IFFT_Radix4(X_fft, Y_fft, N, M)
tmp_fft = np.fft.fft(X+1j*Y)
tmp_ifft = np.fft.ifft(tmp_fft)
print('input:',X)
# print('fft2 :',np.real(tmp_fft))
# print('fft4 :',X_fft)
print('ifft2:',np.real(tmp_ifft))
print('ifft4:',X_ifft/16)
运行结果:
第二个版本更容易理解:
首先是正变换函数:
def FFT_Radix4_v2(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)
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
xa, xb, xc, xd = X_[i0], X_[i1], X_[i2], X_[i3]
ya, yb, yc, yd = Y_[i0], Y_[i1], Y_[i2], Y_[i3]
xa_ = xa + xb + xc + xd
ya_ = ya + yb + yc + yd
xb_ = (xa+yb-xc-yd)* co1 + (ya-xb-yc+xd)* (si1)
yb_ = (ya-xb-yc+xd)* co1 - (xa+yb-xc-yd)* (si1)
xc_ = (xa-xb+xc-xd)* co2 + (ya-yb+yc-yd)* (si2)
yc_ = (ya-yb+yc-yd)* co2 - (xa-xb+xc-xd)* (si2)
xd_ = (xa-yb-xc+yd)* co3 + (ya+xb-yc-xd)* (si3)
yd_ = (ya+xb-yc-xd)* co3 - (xa-yb-xc+yd)* (si3)
X_[i0], X_[i1], X_[i2], X_[i3] = xa_, xb_, xc_, xd_
Y_[i0], Y_[i1], Y_[i2], Y_[i3] = ya_, yb_, yc_, yd_
X_, Y_ = bit_reverse(X_, Y_, N)
return X_, Y_
然后是反变换函数,修改和上面的一样。
def IFFT_Radix4_v2(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)
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
xa, xb, xc, xd = X_[i0], X_[i1], X_[i2], X_[i3]
ya, yb, yc, yd = Y_[i0], Y_[i1], Y_[i2], Y_[i3]
xa_ = xa + xb + xc + xd
ya_ = ya + yb + yc + yd
xb_ = (xa-yb-xc+yd)* co1 - (ya+xb-yc-xd)* (si1)
yb_ = (ya+xb-yc-xd)* co1 + (xa-yb-xc+yd)* (si1)
xc_ = (xa-xb+xc-xd)* co2 - (ya-yb+yc-yd)* (si2)
yc_ = (ya-yb+yc-yd)* co2 + (xa-xb+xc-xd)* (si2)
xd_ = (xa+yb-xc-yd)* co3 - (ya-xb-yc+xd)* (si3)
yd_ = (ya-xb-yc+xd)* co3 + (xa+yb-xc-yd)* (si3)
X_[i0], X_[i1], X_[i2], X_[i3] = xa_, xb_, xc_, xd_
Y_[i0], Y_[i1], Y_[i2], Y_[i3] = ya_, yb_, yc_, yd_
X_, Y_ = bit_reverse(X_, Y_, N)
return X_, Y_
进行同样的实验:
import numpy as np
M = 2
N = 4 ** M
X = np.arange(N)
Y = np.arange(N, 2*N)
X = X.astype(np.float)
Y = Y.astype(np.float)
X_fft, Y_fft = FFT_Radix4_v2(X, Y, N, M)
X_ifft, Y_ifft = IFFT_Radix4_v2(X_fft, Y_fft, N, M)
tmp_fft = np.fft.fft(X+1j*Y)
tmp_ifft = np.fft.ifft(tmp_fft)
print('input :',X)
# print('fft2 :',np.real(tmp_fft))
# print('fft4_v2 :',X_fft)
print('ifft2 :',np.real(tmp_ifft))
print('ifft4_v2:',np.round(X_ifft/16))
下面是运行结果:
可以看到,两次结果均正确。
Good Job!