我有一个复数数组(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()
也许零频率项和负频率项的处理方法不同。我不确定这一点,因为在傅里叶变换的任何描述中都没有提到
这里只有两个错误:
应该是
^{pr2}$因为
range
不包括其上限(这与基于0的索引一起,使得实现FFT类型的算法比Matlab更容易一些)。在还有
应该是
符号约定不同;在NumPy中,直接变换有-1j,逆变换有1j
当然,整个过程效率很低,但你可能想自己详细地写出来。NumPy的向量化操作将正常使用,从
被取代的
其他回路也一样。在
相关问题 更多 >
编程相关推荐