作者热门文章
- c - 在位数组中找到第一个零
- linux - Unix 显示有关匹配两种模式之一的文件的信息
- 正则表达式替换多个文件
- linux - 隐藏来自 xtrace 的命令
我正在尝试使用 python+numpy 解决一个问题,其中我有一些 类型的函数我需要与另一个函数进行卷积 .为了优化代码,我对f和g进行了fft运算,将它们相乘,然后进行逆变换得到结果。
作为进一步的优化,我意识到,由于移位定理,我可以简单地计算一次 f(x,y,z) 的 fft,然后将其乘以一个取决于 的相位因子。获得的fft .特别是, ,其中 N 是 x 和 y 的长度。
我试图用 python+numpy 实现这个简单的公式,但由于某些目前我还不清楚的原因它失败了,所以我向 SO 社区寻求帮助,以找出我遗漏了什么.
我也提供了一个简单的例子。
In [1]: import numpy as np
In [2]: x = np.arange(-10, 11)
In [3]: base = np.fft.fft(np.cos(x))
In [4]: shifted = np.fft.fft(np.cos(x-1))
In [5]: w = np.fft.fftfreq(x.size)
In [6]: phase = np.exp(-2*np.pi*1.0j*w/x.size)
In [7]: test = phase * base
In [8]: (test == shifted).all()
Out[8]: False
In [9]: shifted/base
Out[9]:
array([ 0.54030231 -0.j , 0.54030231 -0.23216322j,
0.54030231 -0.47512034j, 0.54030231 -0.7417705j ,
0.54030231 -1.05016033j, 0.54030231 -1.42919168j,
0.54030231 -1.931478j , 0.54030231 -2.66788185j,
0.54030231 -3.92462627j, 0.54030231 -6.74850534j,
0.54030231-20.55390586j, 0.54030231+20.55390586j,
0.54030231 +6.74850534j, 0.54030231 +3.92462627j,
0.54030231 +2.66788185j, 0.54030231 +1.931478j ,
0.54030231 +1.42919168j, 0.54030231 +1.05016033j,
0.54030231 +0.7417705j , 0.54030231 +0.47512034j,
0.54030231 +0.23216322j])
In [10]: np.abs(shifted/base)
Out[10]:
array([ 0.54030231, 0.58807001, 0.71949004, 0.91768734,
1.18100097, 1.52791212, 2.00562555, 2.72204338,
3.96164334, 6.77009977, 20.56100612, 20.56100612,
6.77009977, 3.96164334, 2.72204338, 2.00562555,
1.52791212, 1.18100097, 0.91768734, 0.71949004, 0.58807001])
我希望通过 shifted/base
我可以获得相位因子的相应值,但可以看出,它不可能是相位因子,因为它的 np. abs
是 >= 1!
最佳答案
由于对 np.fft.fftfreq()
的返回值的不正确(我的错)解释,我的代码中的问题都在输入行 6 上,以及填充的必要性阵列以获得良好的结果。
以下代码效果很好,可以扩展到多维。
In [1]: import numpy as np
In [2]: shift = 1
In [3]: dx = 0.5
In [4]: pad = 20
In [5]: x = np.arange(-10, 11, dx)
In [6]: y = np.cos(x)
In [7]: y = np.pad(y, (0,pad), 'constant')
In [8]: y_shift = np.cos(x-shift)
In [9]: y_fft = np.fft.fft(y)
In [10]: w = np.fft.fftfreq(y.size, dx)
In [11]: phase = np.exp(-2.0*np.pi*1.0j*w*shift)
In [12]: test = phase * y_fft
In [13]: # we use np.real since the resulting inverse fft has small imaginary part values that are zero
In [14]: inv_test = np.real(np.fft.ifft(test))
In [15]: np.allclose(y[:-pad-2],inv_test[2:-pad])
Out[15]: True
关于python - 离散傅立叶变换中的移位定理,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/36885972/
我是一名优秀的程序员,十分优秀!