递归逆 FFT
Recursive Inverse FFT
我已经在递归模式下实现了两个函数 FFT 和 InverseFFT。
这些是函数:
def rfft(a):
n = a.size
if n == 1:
return a
i = 1j
w_n = e ** (-2 * i * pi / float(n))
w = 1
a_0 = np.zeros(int(math.ceil(n / 2.0)), dtype=np.complex_)
a_1 = np.zeros(n / 2, dtype=np.complex_)
for index in range(0, n):
if index % 2 == 0:
a_0[index / 2] = a[index]
else:
a_1[index / 2] = a[index]
y_0 = rfft(a_0)
y_1 = rfft(a_1)
y = np.zeros(n, dtype=np.complex_)
for k in range(0, n / 2):
y[k] = y_0[k] + w * y_1[k]
y[k + n / 2] = y_0[k] - w * y_1[k]
w = w * w_n
return y
def rifft(y):
n = y.size
if n == 1:
return y
i = 1j
w_n = e ** (2 * i * pi / float(n))
w = 1
y_0 = np.zeros(int(math.ceil(n / 2.0)), dtype=np.complex_)
y_1 = np.zeros(n / 2, dtype=np.complex_)
for index in range(0, n):
if index % 2 == 0:
y_0[index / 2] = y[index]
else:
y_1[index / 2] = y[index]
a_0 = rifft(y_0)
a_1 = rifft(y_1)
a = np.zeros(n, dtype=np.complex_)
for k in range(0, n / 2):
a[k] = (a_0[k] + w * a_1[k]) / n
a[k + n / 2] = (a_0[k] - w * a_1[k]) / n
w = w * w_n
return a
根据IFFT的定义,将FFT函数转换为IFFT函数,只需将2*i*pi
改为-2*i*pi
,再除以N
即可。 rfft()
函数工作正常,但 rifft()
函数在这些修改后不起作用。
我将我的函数的输出与 scipy.fftpack.fft
和 scipy.fftpack.ifft
函数进行比较。
我输入以下 NumPy 数组:
a = np.array([1, 0, -1, 3, 0, 0, 0, 0])
下框显示rfft()
函数和scipy.fftpack.fft
函数的结果。
//rfft(a)
[ 3.00000000+0.j -1.12132034-1.12132034j 2.00000000+3.j 3.12132034-3.12132034j -3.00000000+0.j 3.12132034+3.12132034j 2.00000000-3.j -1.12132034+1.12132034j]
//scipy.fftpack.fft(a)
[ 3.00000000+0.j -1.12132034-1.12132034j 2.00000000+3.j 3.12132034-3.12132034j -3.00000000+0.j 3.12132034+3.12132034j 2.00000000-3.j -1.12132034+1.12132034j]
此框显示 rifft()
函数和 scipy.fftpack.ifft
函数的结果。
//rifft(a)
[ 0.04687500+0.j -0.01752063+0.01752063j 0.03125000-0.046875j 0.04877063+0.04877063j -0.04687500+0.j 0.04877063-0.04877063j 0.03125000+0.046875j -0.01752063-0.01752063j]
//scipy.fftpack.ifft(a)
[ 0.37500000+0.j -0.14016504+0.14016504j 0.25000000-0.375j 0.39016504+0.39016504j -0.37500000+0.j 0.39016504-0.39016504j 0.25000000+0.375j -0.14016504-0.14016504j]
除以大小 N
是一个全局比例因子,应该对递归的结果执行,而不是像您所做的那样在递归的每个阶段进行除法(按递减因子在递归中更深入;总体上将结果缩小太多)。您可以通过删除原始实现的最终循环中的 / n
因子来解决此问题,该因子由执行缩放的另一个函数调用:
def unscaledrifft(y):
...
for k in range(0, n / 2):
a[k] = (a_0[k] + w * a_1[k])
a[k + n / 2] = (a_0[k] - w * a_1[k])
w = w * w_n
return a
def rifft(y):
return unscaledrifft(y)/y.size
或者,由于您正在执行基数 2 FFT,全局因子 N
将是 2 的幂,使得 N=2**n
,其中 n
是步数在递归中。因此,您可以在递归的每个阶段除以 2 以获得相同的结果:
def rifft(y):
...
for k in range(0, n / 2):
a[k] = (a_0[k] + w * a_1[k]) / 2
a[k + n / 2] = (a_0[k] - w * a_1[k]) / 2
w = w * w_n
return a
我已经在递归模式下实现了两个函数 FFT 和 InverseFFT。
这些是函数:
def rfft(a):
n = a.size
if n == 1:
return a
i = 1j
w_n = e ** (-2 * i * pi / float(n))
w = 1
a_0 = np.zeros(int(math.ceil(n / 2.0)), dtype=np.complex_)
a_1 = np.zeros(n / 2, dtype=np.complex_)
for index in range(0, n):
if index % 2 == 0:
a_0[index / 2] = a[index]
else:
a_1[index / 2] = a[index]
y_0 = rfft(a_0)
y_1 = rfft(a_1)
y = np.zeros(n, dtype=np.complex_)
for k in range(0, n / 2):
y[k] = y_0[k] + w * y_1[k]
y[k + n / 2] = y_0[k] - w * y_1[k]
w = w * w_n
return y
def rifft(y):
n = y.size
if n == 1:
return y
i = 1j
w_n = e ** (2 * i * pi / float(n))
w = 1
y_0 = np.zeros(int(math.ceil(n / 2.0)), dtype=np.complex_)
y_1 = np.zeros(n / 2, dtype=np.complex_)
for index in range(0, n):
if index % 2 == 0:
y_0[index / 2] = y[index]
else:
y_1[index / 2] = y[index]
a_0 = rifft(y_0)
a_1 = rifft(y_1)
a = np.zeros(n, dtype=np.complex_)
for k in range(0, n / 2):
a[k] = (a_0[k] + w * a_1[k]) / n
a[k + n / 2] = (a_0[k] - w * a_1[k]) / n
w = w * w_n
return a
根据IFFT的定义,将FFT函数转换为IFFT函数,只需将2*i*pi
改为-2*i*pi
,再除以N
即可。 rfft()
函数工作正常,但 rifft()
函数在这些修改后不起作用。
我将我的函数的输出与 scipy.fftpack.fft
和 scipy.fftpack.ifft
函数进行比较。
我输入以下 NumPy 数组:
a = np.array([1, 0, -1, 3, 0, 0, 0, 0])
下框显示rfft()
函数和scipy.fftpack.fft
函数的结果。
//rfft(a)
[ 3.00000000+0.j -1.12132034-1.12132034j 2.00000000+3.j 3.12132034-3.12132034j -3.00000000+0.j 3.12132034+3.12132034j 2.00000000-3.j -1.12132034+1.12132034j]
//scipy.fftpack.fft(a)
[ 3.00000000+0.j -1.12132034-1.12132034j 2.00000000+3.j 3.12132034-3.12132034j -3.00000000+0.j 3.12132034+3.12132034j 2.00000000-3.j -1.12132034+1.12132034j]
此框显示 rifft()
函数和 scipy.fftpack.ifft
函数的结果。
//rifft(a)
[ 0.04687500+0.j -0.01752063+0.01752063j 0.03125000-0.046875j 0.04877063+0.04877063j -0.04687500+0.j 0.04877063-0.04877063j 0.03125000+0.046875j -0.01752063-0.01752063j]
//scipy.fftpack.ifft(a)
[ 0.37500000+0.j -0.14016504+0.14016504j 0.25000000-0.375j 0.39016504+0.39016504j -0.37500000+0.j 0.39016504-0.39016504j 0.25000000+0.375j -0.14016504-0.14016504j]
除以大小 N
是一个全局比例因子,应该对递归的结果执行,而不是像您所做的那样在递归的每个阶段进行除法(按递减因子在递归中更深入;总体上将结果缩小太多)。您可以通过删除原始实现的最终循环中的 / n
因子来解决此问题,该因子由执行缩放的另一个函数调用:
def unscaledrifft(y):
...
for k in range(0, n / 2):
a[k] = (a_0[k] + w * a_1[k])
a[k + n / 2] = (a_0[k] - w * a_1[k])
w = w * w_n
return a
def rifft(y):
return unscaledrifft(y)/y.size
或者,由于您正在执行基数 2 FFT,全局因子 N
将是 2 的幂,使得 N=2**n
,其中 n
是步数在递归中。因此,您可以在递归的每个阶段除以 2 以获得相同的结果:
def rifft(y):
...
for k in range(0, n / 2):
a[k] = (a_0[k] + w * a_1[k]) / 2
a[k + n / 2] = (a_0[k] - w * a_1[k]) / 2
w = w * w_n
return a