Как стать автором
Обновить

Опыт применения теоремы о свертке с использованием scipy.fft

Время на прочтение 4 мин
Количество просмотров 7K

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

\mathcal{F}\{f*g\} = \mathcal{F}\{f\}\cdot\mathcal{F}\{g\}.

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

(f*g)(x)\equiv\int\limits_{-\infty}^{\infty}f(t)g(x-t)dt=\mathcal{F}^{-1}\bigl\{\mathcal{F}\{f\}\cdot\mathcal{F}\{g\}\bigr\}.

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

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

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

f(x)=\frac{1}{\pi\sqrt{1-x^2}}\hspace{5mm}\text{и}\hspace{5mm}g(x)=\frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{x^2}{2\sigma^2}\right),

причём имеется в виду, что нас интересует только действительная часть области значений функцииf(x), то есть область её определения лежит в интервале x \in [-1,1]. Интеграл функции в этом интервале равен 1, а преобразование Фурье вычисляется аналитически:

\mathcal{F}(f(x))= \int\limits_{-1}^{+1}\frac{\exp(-2\pi i u x)}{\pi\sqrt{1-x^2}}dx = J_0 (2\pi u),

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

\mathcal{F}(g(x)) = \exp(-2\,\pi^2\sigma^2u^2).

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

Выберем область определения x\in[-R+D_x/2,R-D_x/2], где D_xиR- положительные вещественные числа, которые мы будем считать размерными, например, пусть это будут сантиметры. Разобъём этот интервал на N_xодинаковых отрезков, длина каждого из которых равна D_x = 2R/N_x. Интервал D_x\ll Rназывается шагом дискретизации и определяет расстояние между соседними узлами разбиения шкалы x. Обратная величина называется частотой дискретизации и в нашем примере измеряется в обратных сантиметрах.

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

u=\Bigl[0,\;\frac{1}{2R},\;\dots,\;\frac{(N_x-1)/2}{2R},\; -\frac{N_x/2}{2R},\;\dots,\; -\frac{1}{2R} \Bigr].

Если функция, которую мы ищем, вещественна (не имеет мнимой части), положительные частоты в последовательности являются сопряженными величинами отрицательных частот, за счёт этого можно примерно в два раза сократить время вычислений. Для получения укороченной (без отрицательных частот) последовательности используется функция 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()

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

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

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

\mathcal{F}(\vec{\omega}) = \mathcal{F}(\omega_x,\omega_y)=\int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty}f(x,y)\exp(-i(\omega_x x+\omega_y y))\,dx\,dy\,,

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

f(\vec{r}) = f(x,y) = \frac{1}{(2\pi)^2}\int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty}\mathcal{F}(\omega_x, \omega_y)\exp(i\,\vec{\omega}\cdot\vec{r})\,d\omega_x\,d\omega_y\,,

где для краткости используются обозначения\vec{\omega}=(\omega_x,\omega_y),\;\vec{r}=(x,y). Переход к полярным координатам осуществляется путём замены переменных x=r\cos\varphi,\,y=r\sin\varphiв координатной плоскости и  \omega_x=\rho\cos\psi,\,\omega_y=\rho\sin\psiв частотной плоскости. В результате такой замены, преобразование Фурье в полярных координатах выглядит следующим образом:

\mathcal{F}(\rho,\psi)=\int\limits_{0}^{\infty}\int\limits_{-\pi}^{\pi}f(r,\varphi)\exp(-i\,r\rho\cos(\psi-\varphi))\,r\,dr\,d\varphi.

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

\mathcal{F}(\rho,\psi)=\int\limits_{0}^{\infty}r\,f(r)\,dr \int\limits_{-\pi}^{\pi}\exp(-i\,r\rho\cos(\psi-\varphi))\,d\varphi.

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

J_0(x) = \frac{1}{2\pi}\int\limits_{-pi}^{\pi}\exp(-i\,x\,\cos\alpha)d\alpha,

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

\mathcal{F}(\rho) = 2\pi\int\limits_{0}^{\infty}f(r)J_0(\rho\,r)\,r\,dr = 2\pi\,\mathbb{H}_0\{f(r)\}\,,

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


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

f(r)=\frac{1}{\pi\sqrt{1-r^2}}\hspace{5mm}\text{и}\hspace{5mm}g(x,y)=\frac{1}{2\pi\sigma_x\sigma_y}\exp\left(-\frac{x^2}{2\sigma_x^2} - \frac{y^2}{2\sigma_y^2}\right),

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

\mathcal{F}(\rho) = 2\int\limits_{0}^{1}\frac{r\,J_0(r\rho)}{\sqrt{1-r^2}}dr = 2\frac{\sin(\rho)}{\rho},

где \rho = 2\pi\sqrt{\nu_x^2+\nu_y^2}. В иных случаях преобразование Ганкеля можно делать численно с помощью функции 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!

Теги:
Хабы:
+6
Комментарии 5
Комментарии Комментарии 5

Публикации

Истории

Работа

Python разработчик
136 вакансий
Data Scientist
66 вакансий

Ближайшие события

Московский туристический хакатон
Дата 23 марта – 7 апреля
Место
Москва Онлайн
Геймтон «DatsEdenSpace» от DatsTeam
Дата 5 – 6 апреля
Время 17:00 – 20:00
Место
Онлайн