scipy.signal.fftconvolve가 필요한 결과를 제공하지 않습니다.

Question

python의 fftconvolve과 관련하여 질문이 있습니다. 나의 현재 연구에서는 두 함수 사이의 회선을 계산해야했습니다. 이렇게하려면 푸리에 변환을 사용하여 계산합니다 (이는 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) 

당신은 또한 w 축에 변화를 볼 수 있습니다 JH의 실제/허수 부분을 플롯하면 어떻게 든 얻기 위해 J의 결과를 곱했다 : 여기 내 코드입니다 정확한 결과를 닫을 수는 있습니다.

제안 사항?

감사합니다.

Answer 1

컨디션을 계산할 때 경계 조건이 중요합니다.

두 개의 신호를 컨볼 루션하면 결과의 가장자리는 외부의 어떤 값으로 입력 가장자리에 종속됩니까? fftconvolve은 0으로 채워진 경계를 가정하고 컨볼 루션을 계산합니다.

source code of 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 제로 패드 주기적 경계 조건을 방지하기 위해 각각의 신호 : 순방향의 FFT의 곱의 역 FFT가 패딩 된 결과를 제공

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

:

을

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' 옵션을 사용합니다).