numpy中的一维逆离散傅里叶变换是如何计算的？

Question

我有一个复数数组 (waves)，其排序方式与 fft 返回的方式相同。如果我在这个数组上调用 ifft，它 returns 是原始样本的近似值。

我想在 python 中自己实现 ifft。我找到了 the formula of IFFT. I implemented it, but it looks a bit different from the ifft result. I tried to fix it by looking at the ifft source，但这是一个高度优化的版本，我无法弄清楚它是如何工作的

这是我目前拥有的：

import numpy as np
import matplotlib.pyplot as plt

from scipy.fftpack import fft

# create samples
data = np.zeros(1024)
for k in range(0, 1024):
    data[k] = np.sin(k/20 * np.pi*2)

# plot original samples
#plt.plot(data)
#plt.show()

# apply FFT on samples
waves = fft(data)

# plot FFT result
#plt.plot(np.imag(waves),'r')
#plt.plot(np.real(waves),'b')
#plt.ylabel('value')
#plt.xlabel('period')
#plt.show()

# 
res = np.zeros(1024)
for k in range(0, 1024):
    val = 0.0
    for n in range(0, len(waves)-1):
        # https://dsp.stackexchange.com/a/510/25943
        val += waves[n]*np.exp(-1.j * 2*np.pi * n * k / len(waves)) / len(waves)
    res[k] = val.real


#np implementation
res2 = np.fft.ifft(waves)

plt.plot(data, 'b') # original
plt.plot(res,'g') # my implementation
plt.plot(res2,'y') # np implementation
plt.show()

也许必须以不同方式处理零频率项和负频率项。这个我不太确定，因为在任何傅里叶变换的描述中都没有提到

Answer 1

这里只有两个错误：

for n in range(0, len(waves)-1):

应该是

for n in range(0, len(waves)):

因为 range 不包括其上限（与基于 0 的索引一起，与 Matlab 相比，实现 FFT 类型的算法更容易一些）。

此外，

val += waves[n]*np.exp(-1.j * 2*np.pi * n * k / len(waves))

应该是

val += waves[n]*np.exp(1.j * 2*np.pi * n * k / len(waves))

签名约定不同；在 NumPy 中，直接变换有 -1j，逆变换有 1j。

当然，整个过程效率很低，但你可能想自己详细写出来。 NumPy 的矢量化操作将正常使用，以

开头

data = np.zeros(1024)
for k in range(0, 1024):
    data[k] = np.sin(k/20 * np.pi*2)

被

取代

data = np.sin(np.arange(1024)/20 * np.pi*2)

其他循环也是如此。

numpy中的一维逆离散傅里叶变换是如何计算的？

How is the one-dimensional inverse discrete Fourier Transform calculated in numpy?

python

numpy

fft

ifft