Оценка свертки двух непрерывных функций с использованием fftconvolve

Я пытаюсь оценить свертку двух непрерывных функций, используя scipy.signal.fftconvolve. Сценарий кода следующий:

Я пытаюсь аппроксимировать следующий двойной интеграл:

введите здесь описание изображения

, то есть в области C_1(x',y'), представляющей круг радиусом 1 с центром в (x', y'). Это можно аппроксимировать следующим интегралом:

введите здесь описание изображения

где функция K выбрана как непрерывная интегрируемая функция, скажем, exp(-x^2-y^2), форма которой приблизительно соответствует форме окружности радиуса 1. Если я возьму функцию K'(x,y)=K(-x,-y ), то интеграл есть в точности свертка двух функций:

введите здесь описание изображения

Поэтому я пытаюсь разбить эти две функции на массивы, а затем выполнить свертку.

Следующий код будет написан на Julia, а функция fftconvolve будет импортирована с использованием PyCall.jl.

using PyCall
using Interpolations

r = 1
xc = -10:0.05:10
yc = -10:0.05:10

K(x, y) = exp(-(x^2+y^2)/r^2)
rho(x, y) = x^2+y^3    # Try some arbitrary function

ss = pyimport("scipy.signal")  # Import scipy.signal module from Python
a = [rho(x,y) for x in xc, y in yc]
b = [K(-x,-y) for x in xc, y in yc]
c = ss.fftconvolve(a,b,mode="same")  # zero-paddings beyond boundary, unimportant since rho is near zero beyond the boundary anyway
c_unscaled = interpolate(c', BSpline(Cubic(Line(OnCell())))) 
# Adjoint because the array comprehension switched x and y, then interpolate the array
c_scaled = Interpolations.scale(c_unscaled, xc, yc) # Scale the interpolated function w.r.t. xc, yc
print(c_scaled(0.0,0.0)) # The result of the integral for (x', y') = (0, 0)

Результат равен 628.3185307178969, а результат численного интегрирования равен 0.785398. В чем проблема?


python scipy julia convolution fft
person Sato    schedule 20.10.2019    source источник
comment
Возможно, вы захотите нормализовать K, чтобы его интеграл был равен 1. Затем вам нужно разделить на dx * dy, которое в вашем коде не равно 1.   -  person Cris Luengo    schedule 20.10.2019
comment
Там, где я сказал разделить, я имел в виду умножить, конечно. :)   -  person Cris Luengo    schedule 21.10.2019

Ответы (1)


arrow_upward
0
arrow_downward

Вероятно, вы могли бы попытаться использовать scipy.signal.convolve, который будет свертываться из двух N-мерных массивов, но не с использованием быстрого преобразования Фурье. Он использует прямой метод для вычисления свертки. Здесь я имею в виду, что свертка определяется непосредственно из сумм.

Таким образом, вы могли бы попробовать заменить строку, в которой вы вычисляете c, на эту:

c = ss.convolve(a,b,mode="same", method='direct')
person Nikola Zubic    schedule 20.10.2019
comment
Но я не понимаю, почему использование БПФ должно отличаться для получения точного результата? Насколько я слышал, БПФ должен иметь лучшую производительность для больших массивов? Когда я запускаю ss.convolve, кажется, что он работает очень медленно, что, возможно, даже использование пакета численного интегрирования может быть даже лучшей альтернативой. Я думаю, не заключается ли проблема в моем выборе xc и yc. Чем меньше интервал, тем больше результат. Я думаю, следует ли мне, скажем, умножить результат на 0,02, если я выберу xc и yc будет -10:0.05:10. - person Sato; 20.10.2019
comment
Мне кажется, что полученный результат совпадает с результатом численного интеграла, если я умножу свертку на 0.02 ** 2 и выберу более узкую область (rho(x, y) = x^2+y^3, которую я выбрал, расходится, когда x и y становятся больше, поэтому rho(x,y) вносит вклад в свертку, даже когда x и y далеко от центра) - person Sato; 20.10.2019
comment
Однако, когда я попробовал некоторые другие K и rho, скажем, g(x, y) = sin(x)*cos(2*y) и K(x,y) = x*y*exp(-x^2-y^2)*(1-exp(-10*y))^(-1) (логистическая сигмовидная функция предназначена для эмуляции ступенчатой ​​функции), она больше не работает. Поэтому я хотел бы знать, является ли мое использование fftconvolve концептуально неправильным. - person Sato; 20.10.2019