3

パイソンについて質問fftconvolveです。私の現在の研究では、2 つの関数間の畳み込みを計算する必要がありました。そうするために、フーリエ変換を使用して計算しています（これを使用numpy.fftして正規化しました）。問題は、パッケージを使用して比較したい場合fftconvolve、正しい結果が得られないことです。これが私のコードです：

#!/usr/bin/python
import numpy as np
from scipy.signal import fftconvolve , convolve 

def FFT(array , sign):
  if sign==1:
    return np.fft.fftshift(np.fft.fft(np.fft.fftshift(array))) * dw / (2.0 * np.pi)
  elif sign==-1:
    return np.fft.fftshift(np.fft.ifft(np.fft.fftshift(array))) * dt * len(array)


def convolve_arrays(array1,array2,sign):
  sign = int(sign)
  temp1 = FFT(array1 , sign,)
  temp2 = FFT(array2 , sign,)
  temp3 = np.multiply(temp1 , temp2)
  return  FFT(temp3 , -1 * sign) / (2. * np.pi) 

""" EXAMPLE """ 

dt    = .1
N     = 2**17
t_max = N * dt / 2
time  = dt * np.arange(-N / 2 , N / 2 , 1)

dw    = 2. * np.pi / (N * dt)
w_max = N * dw / 2.
w     = dw * np.arange(-N / 2 , N / 2 , 1)

eta_fourier = 1e-10




Gamma   = 1.
epsilon = .5
omega   = .5


G    = zeros(N , complex)
G[:] = 1. / (w[:] - epsilon + 1j * eta_fourier)

D    = zeros(N , complex)
D[:] = 1. / (w[:] - omega + 1j * eta_fourier) - 1. / (w[:] + omega + 1j * eta_fourier)

H    = convolve_arrays(D , G , 1)     
J    = fftconvolve(D , G , mode = 'same') * np.pi  / (2. * N)

の実数部/虚数部をプロットするとH、軸にJシフトが見られます。また、の結果を乗算して、何らかの方法で正しい結果に近づける必要がありました (ただし、まだそうではありません)。wJ

助言がありますか？

ありがとう！

4

1 に答える 1

6

畳み込みを計算する場合、境界条件は重要です。

2 つの信号をたたみ込む場合、結果のエッジは、入力のエッジの外側で仮定する値によって異なります。fftconvolveゼロパディングされた境界を想定して畳み込みを計算します。

fftconvolveのソース コードを見てください。ゼロでパディングされた境界条件を達成するために彼らが経験する悪ふざけに注意してください。特に、次の行です。

size = s1 + s2 - 1

...

fsize = 2 ** np.ceil(np.log2(size)).astype(int) #For speed; often suboptimal!
fslice = tuple([slice(0, int(sz)) for sz in size])

...

ret = ifftn(fftn(in1, fsize) * fftn(in2, fsize))[fslice].copy()

...

return _centered(ret, s1) #strips off padding

これは良いものです！fftconvolveフーリエに基づく畳み込みを理解したい場合は、 のコードを注意深く読む価値があります。

簡単なスケッチ

順方向 FFT は、周期的な境界条件を防ぐために各信号をゼロ パディングします。

a = np.array([3, 4, 5])
b = np.fft.ifftn(np.fft.fftn(a, (5,))).real
print b #[ 3.  4.  5.  0.  0.]

フォワード FFT の積の逆 FFT は、パディングされた結果を与えます。

a = np.array([3, 4, 5])
b = np.array([0., 0.9, 0.1])
b = np.fft.ifftn(np.fft.fftn(a, (5,))*
                 np.fft.fftn(b, (5,))
                 ).real
print b #[ 0.   2.7  3.9  4.9  0.5]

関数は_centered最後に余分なパディング ピクセルを取り除きます (mode='same'オプションを使用すると仮定します)。

于 2014-01-04T23:33:20.670 に答える