我有一个复数数组(waves),其排列顺序与fft返回的顺序相同。如果在此数组上调用ifft,它将返回原始样本的近似值。

我想在python中自己实现ifft。我找到了the formula of IFFT。我实现了它,但是它与ifft结果看起来有些不同。我试图通过查看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()


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

最佳答案

这里只有两个错误:

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


应该

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


因为range不包括其上限(与基于0的索引一起,使得实现FFT类型算法比Matlab容易一些)。

也,

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)


和其他循环类似。

关于python - 如何以numpy计算一维逆离散傅立叶逆变换?,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/47001893/

10-13 07:28