Теорема о свёртке утверждает, что преобразование Фурье от свёртки двух функций является произведением их Фурье образов:

Как следствие, свёртку двух функций можно вычислить путём обратного преобразования Фурье:

Когда свёртку не удаётся получить аналитически, использование теоремы позволяет вместо численного интегрирования воспользоваться алгоритмом быстрого преобразования Фурье и существенно ускорить процесс вычислений. В ходе решения возникшей у меня задачи мне не удалось найти в интернете достаточно подробного описания алгоритма действий, поэтому я решил изложить его здесь, во-первых, как памятку для себя и во-вторых, на случай если это окажется полезным кому-то ещё. Не претендуя на общность изложения, ниже я приведу примеры найденных решений для функций от одной и от двух переменных.

Функция от одной переменной

Допустим, нам нужно вычислить свёртку функций

причём имеется в виду, что нас интересует только действительная часть области значений функции, то есть область её определения лежит в интервале . Интеграл функции в этом интервале равен 1, а преобразование Фурье вычисляется аналитически:

где - функция Бесселя первого рода. Фурье образ функции Гаусса:

Дискретным преобразованием Фурье (ДПФ) называется метод вычислений, в котором и функция и ее преобразование заменяются дискретными аналогами. Быстрое преобразование Фурье (БПФ) относится к способу эффективного вычисления ДПФ с использованием симметрии в вычисляемых условиях. Симметрия наивысшая, когда число узлов, в которых вычисляется значение функции, и число гармоник,  являются степенью двойки. В библиотеке SciPy с открытым исходным кодом для работы с БПФ предназначен модуль scipy.fft. Для того, чтобы вычислить искомую свёртку методом БФП, нужно задать интересующую нас область определения результата и задать количество узлов разбиения, в которых мы хотим получить значение свёртки.

Выберем область определения , где и- положительные вещественные числа, которые мы будем считать размерными, например, пусть это будут сантиметры. Разобъём этот интервал на одинаковых отрезков, длина каждого из которых равна . Интервал называется шагом дискретизации и определяет расстояние между соседними узлами разбиения шкалы . Обратная величина называется частотой дискретизации и в нашем примере измеряется в обратных сантиметрах.

В шкале частот нашему диапазону шкалыбудет соответствовать набор из положительных и отрицательных гармоник разложения Фурье. Если - чётное, последовательность частот, которая генерируется функцией scipy.fft.fftfreq, будет такая:

Если функция, которую мы ищем, вещественна (не имеет мнимой части), положительные частоты в последовательности являются сопряженными величинами отрицательных частот, за счёт этого можно примерно в два раза сократить время вычислений. Для получения укороченной (без отрицательных частот) последовательности используется функция scipy.fft.rfftfreq. Для выполнения прямого и обратного преобразований Фурье для вещественного сигнала предназначены функции scipy.fft.rfft и scipy.fft.irfft соответственно. Функция scipy.fft.fftshift располагает компоненту с нулевой частотой в середине частотного спектра, а scipy.fft.ifftshift осуществляет обратное преобразование. Ниже приведён код для вычисления свёртки:

import numpy as np
from scipy.special import j0 as BesselJ0 
import scipy.fft as FFT
from matplotlib import pyplot as plt

π    = np.pi
σ    = 0.01
R    = 1.5    # интервал, в котором вычисляется свёртка
Nx   = 4096   # число разбиений по оси x
Dx   = 2*R/Nx # шаг дискретизации
x    = np.linspace(0.5*Dx-R, R-0.5*Dx, num = Nx)
u    = FFT.rfftfreq(Nx//16, d = 16*Dx) # цифровой ФНЧ
J0   = BesselJ0(2*π*u)*np.exp(-2*(π*σ*u)**2)
y    = FFT.ifftshift(FFT.irfft(J0, n=Nx, norm='forward'))

plt.subplot(211)
plt.plot(u, J0, 'r-', label='спектр частот')
plt.xlabel('u, [1/см]')
plt.grid(True); plt.legend()
plt.subplot(212)
plt.plot(x, y, 'r-', label='результат свёртки')
plt.xlabel('x, [см]')
plt.grid(True); plt.legend()
plt.show()

Результат представлен на картинке. Как видно из спектра частот, амплитуда гармоник экспоненциально падает с ростом частоты (из-за , см. выше) , поэтому число вычисляемых гармоник уменьшено в 16 раз с помощью "цифрового ФНЧ".

Функция от двух переменных

Перейдем к рассмотрению двумерного случая, основываясь на статье Natalie Baddour. Фурье преобразование функцииопределяется по аналогии с одномерным случаем:

где- круговые частоты в плоскости Фурье-образа функции. Обратное преобразование Фурье выглядит так:

где для краткости используются обозначения. Переход к полярным координатам осуществляется путём замены переменных в координатной плоскости и в частотной плоскости. В результате такой замены, преобразование Фурье в полярных координатах выглядит следующим образом:

Если функция не зависит от угла, т. е. , то её не нужно интегрировать по и предыдущее выражение упрощается:

Используя интегральное определение функции Бесселя нулевого порядка

наше Фурье преобразование можно записать так:

гденазывается преобразованием Ганкеля (Hankel transform) нулевого порядка.

Для примера вычислим свёртку двух функций:

одна из которых задана в полярных координатах, а другая - в декартовых. Фурье-образ функции в нашем случае легко вычисляется аналитически:

где . В иных случаях преобразование Ганкеля можно делать численно с помощью функции scipy.fft.fht. Пользуясь полученными выше формулами и рецептами, создадим код для вычисления свёртки:

import numpy as np
from scipy import fft as FFT
import matplotlib.pyplot as plt
π    = np.pi

σx   = 0.05 
σy   = 0.02
R    = 1.5
N    = 1024
D    = 2*R/N
s    = np.linspace(0.5*D-R, R-0.5*D, num = N)
f    = FFT.fftshift(FFT.fftfreq(N, d = D))
x,y  = np.meshgrid(s,s)
u,v  = np.meshgrid(f,f)
J    = 2*np.sinc(2*(u**2 + v**2)**0.5)*np.exp(-2*π*π*((σx*u)**2 + (σy*v)**2))
z    = np.abs(FFT.ifftshift(FFT.irfft2(J, s=[N,N], workers=-1, norm='forward')))

plt.figure(figsize=(16, 9))

levels = np.linspace(J.min(), J.max(), 16)
plt.subplot(121)
plt.contourf(u, v, np.real(J), levels)
plt.colorbar(); plt.title('Фурье образ')
plt.xlim(-5,5); plt.xlabel('u, [1/см]')
plt.ylim(-5,5); plt.ylabel('v, [1/см]')
plt.subplot(122)
levels = np.linspace(z.min(), z.max(), 16)
plt.contourf(x, y, z, levels)
plt.colorbar(); plt.title('Результат свёртки')
plt.xlabel('x, [см]')
plt.ylabel('y, [см]')
plt.show()

Опция workers=-1 в функции irfft2 позволяет параллельно использовать все имеющиеся в системе процессорные ядра. Результат вычислений представлен на картинке:

Заключение

Вот вроде бы и всё о чём я хотел написать. Буду рад критике, вопросам и предложениям. Всех с наступающим 2022!