2016-04-27 17 views
4

私は別の関数$g(x,y,z,t)$と畳み込む必要があるタイプ$f(x-x_i,y-y_i,z)$のいくつかの関数を持つpython + numpyの問題を解決しようとしています。コードを最適化するために、fとfのfftを実行し、それらを乗算した後、逆変換を実行して結果を得ました。離散フーリエ変換のシフト定理

シフト定理のおかげで、fft(x、y、z)のfftを単純に計算してから、$(x_i, y_i)$に依存する位相係数を乗算してfft $f(x-x_i,y-y_i,z)$。特に、$\mathcal{F}(f(x-x_i,y-y_i,z)) = exp^{-2\pi j (x_i w_1 + y_i w_2)/N} \mathcal{F}(f(x,y,z))$(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]) 

Iは、shifted/baseによってIは位相係数の対応する値を得ることができることを期待するが、見ることができるように、そのnp.absは> = 1であるので、それは、位相因子であることはできません!

答えて

0

私のコードの問題は、返り値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 
関連する問題