Чтобы извлечь максимум пользы из этой статьи, вам следует иметь базовое представление об интегралах и преобразовании Фурье, а также о его свойствах.
Если у вас нет таких знаний, мы рекомендуем ознакомиться с разделом 1, где эти понятия рассматриваются. Если вы уже знакомы с ними, вы можете перейти сразу к разделу 2.
В разделе 2.1 мы реализуем преобразование Фурье, используя метод численной квадратуры. В разделе 2.2 мы реализуем его с помощью алгоритма быстрого преобразования Фурье (БПФ) и формулы интерполяции Шеннона.
Эта статья попала к нам на последнем курсе инженерного факультета.
В рамках курса по калибровке стохастических процессов наш профессор хотел проверить понимание материала. Для этого нам было предложено выбрать научную статью, подробно изучить её и воспроизвести результаты.
Мы выбрали статью Эль Колеи (2013), в которой предлагается параметрический метод оценки параметров модели стохастической волатильности, такой как модель GARCH. Подход основан на минимизации контраста и деконволюции.
Чтобы воспроизвести результаты, нам нужно было реализовать метод оптимизации, который включал вычисление преобразования Фурье функции $f_\theta$:
Для вычисления преобразования Фурье функции $f_\theta$ Эль Колеи (2013) использовал метод левой римановой суммы (также известный как прямоугольная квадратура), реализованный в MATLAB. В статье предполагалось, что использование алгоритма быстрого преобразования Фурье (БПФ) может ускорить вычисления.
Мы решили воспроизвести результаты, используя Python, и реализовали метод левой прямоугольной квадратуры, чтобы проверить, действительно ли это ускоряет вычисления.
1. Определение и свойства преобразования Фурье
Мы принимаем рамки, представленные Балаком (2011), чтобы определить преобразование Фурье функции $f$ и ввести его свойства.
Мы рассматриваем функции, которые принадлежат пространству интегрируемых функций, обозначаемому $L^1(\mathbb{R}, \mathbb{K})$. Это пространство включает все функции $f: \mathbb{R} \to \mathbb{K}$, где $\mathbb{K}$ представляет либо действительные числа ($\mathbb{R}$), либо комплексные числа ($\mathbb{C}$).
Эти функции интегрируемы в смысле Лебега, что означает, что интеграл от их абсолютного значения конечен.
2. Как численно аппроксимировать преобразование Фурье?
В своей статье Балак (2011) показывает, что вычисление преобразования Фурье функции включает в себя аппроксимацию его следующим интегралом по интервалу $[-T/2, T/2]$:
Где $T$ — достаточно большое число, чтобы интеграл сходился. Приближённое значение интеграла $S(\nu)$ можно вычислить с помощью методов квадратуры.
В следующем разделе мы аппроксимируем этот интеграл с помощью метода левой прямоугольной квадратуры (также известного как метод левой римановой суммы).
2.1 Метод квадратуры левых прямоугольников
Чтобы вычислить интеграл $S(\nu)$ с помощью метода левой прямоугольной квадратуры, мы поступаем следующим образом:
Дискретизация интервала. Мы делим интервал $[-T/2, T/2]$ на $N$ равномерных подотрезков длины $h_t = T/N$.
Точки дискретизации, соответствующие левым концам подотрезков, определяются как:
$t_k = -T/2 + h_t, \text{ для } k = 0, 1, \ldots, N-1.$
Аппроксимация интеграла. Используя отношение Шале, мы аппроксимируем интеграл $S(\nu)$ следующим образом:
$S(\nu) \approx \sum_{k=0}^{N-1} f(t_k) h_t e^{-2i\pi \nu t_k}.$
Поскольку $t_{k+1} - t_k = h_t$ и $t_k = -T/2 + k \cdot h_t = T(k/N - 1/2)$, выражение принимает вид:
$S(\nu) \approx T \sum_{k=0}^{N-1} f(t_k) e^{-2i\pi \nu T k/N}.$
Мы называем этот подход методом левой прямоугольной квадратуры, потому что он использует левый конец $t_k$ каждого подотрезка для аппроксимации значения $f(t)$ в этом интервале.
Окончательная формула. Окончательное выражение для аппроксимации преобразования Фурье имеет вид:
$\hat{f}(\nu) \approx T \sum_{k=0}^{N-1} f(t_k) e^{-2i\pi \nu T k/N}.$
2.1.1 Реализация метода левой прямоугольной квадратуры в Python
Функция tfquad ниже реализует метод левой прямоугольной квадратуры для вычисления преобразования Фурье функции $f$ при заданной частоте $\nu$.
import numpy as np
def tfquad(f, nu, n, T):
"""
Вычисляет преобразование Фурье функции f при частоте nu
с использованием левой римановой суммы квадратуры
над интервалом [-T/2, T/2].
Параметры:
----------
f : callable
Функция для преобразования. Должна принимать массив NumPy в качестве входных данных.
nu : float
Частота, при которой вычисляется преобразование Фурье.
n : int
Количество точек квадратуры.
T : float
Ширина временного окна [-T/2, T/2].
Возвращает:
-------
tfnu : complex
Приближённое значение преобразования Фурье при частоте nu.
"""
k = np.arange(n)
t_k = (k / n - 0.5) * T
weights = np.exp(-2j * np.pi * nu * T * k / n)
prefactor = (T / n) * np.exp(1j * np.pi * nu * T)
return prefactor * np.sum(f(t_k) * weights)
Мы также можем использовать функцию quad из SciPy для определения преобразования Фурье функции $f$ при заданной частоте $\nu$. Функция tf_integral ниже реализует этот подход. Она использует численное интегрирование для вычисления преобразования Фурье $f$ по интервалу $[-T/2, T/2]$.
from scipy.integrate import quad
def tf_integral(f, nu, T):
"""Вычисляет FT функции f при частоте nu над [-T/2, T/2] с помощью scipy quad."""
real_part = quad(lambda t: np.real(f(t) * np.exp(-2j * np.pi * nu * t)), -T/2, T/2)[0]
imag_part = quad(lambda t: np.imag(f(t) * np.exp(-2j * np.pi * nu * t)), -T/2, T/2)[0]
return real_part + 1j * imag_part
2.2 Вычисление преобразования Фурье при частоте $\nu$ с помощью алгоритма БПФ
В этом разделе мы обозначаем приближение преобразования Фурье $f(\nu)$ функции $f$ при точке $\nu \in [-\nu_{max}/2, \nu_{max}/2]$ как $\hat{S}n(\nu)$, где $\nu{max} = n/T$.
Теперь я представляю алгоритм преобразования Фурье, используемый для аппроксимации $f(\nu)$. Я не буду вдаваться в подробности алгоритма быстрого преобразования Фурье (БПФ) в этой статье. Для упрощённого объяснения я ссылаюсь на Балака (2011), а для более технического и всестороннего рассмотрения — на оригинальную работу Кули и Туки (1965).
Важно понимать, что использование алгоритма БПФ для аппроксимации преобразования Фурье функции $f$ основано на результате, установленном Эпштейном (2005). Этот результат гласит, что когда $\hat{S}_n(\nu)$ оценивается на частотах $v_j = j/T$, для $j = 0, 1, \ldots, n - 1$, он обеспечивает хорошую аппроксимацию непрерывного преобразования Фурье $f(\nu)$.
Более того, $\hat{S}_n$ известно как периодическое. Эта периодичность присваивает симметричную роль индексам $j \in {0, 1, \ldots, n - 1}$ и $k \in {-n/2, -n/2 + 1, \ldots, -1}$.
В следующем разделе мы рассмотрим, как использовать интерполяцию Шеннона для вычисления преобразования Фурье функции $f$ при любой частоте $\nu$.
Заключение
В этой статье мы рассмотрели два метода аппроксимации преобразования Фурье функции. Первый — это метод численной квадратуры (метод левой прямоугольной квадратуры), а второй — алгоритм быстрого преобразования Фурье (БПФ).
Мы показали, что, в отличие от встроенных функций fft в SciPy или NumPy, БПФ можно адаптировать для вычисления преобразования Фурье непрерывной функции путём правильной формулировки.
Наша реализация основана на работе Балака (2013), который продемонстрировал, как воспроизвести вычисление БПФ в Python. Мы также представили теорему интерполяции Шеннона, которая позволяет нам оценить преобразование Фурье любой функции на $\mathbb{R}$ при произвольных частотах.
Эта интерполяция может быть заменена более традиционными подходами, такими как линейная или сплайн-интерполяция, когда это уместно.
Наконец, важно отметить, что когда требуется преобразование Фурье на одной частоте, часто более эффективно вычислять его напрямую с помощью метода численной квадратуры. Такие методы хорошо подходят для обработки колебательной природы интеграла Фурье и могут быть более точными, чем применение БПФ, а затем интерполяция.
Список литературы
Balac, Stéphane. 2011. “La Transformée de Fourier Vue Sous l’angle Du Calcul Numérique.”
Cooley, James W, and John W Tukey. 1965. “An Algorithm for the Machine Calculation of Complex Fourier Series.” Mathematics of Computation 19 (90): 297–301.
El Kolei, Salima. 2013. “Parametric Estimation of Hidden Stochastic Model by Contrast Minimization and Deconvolution: Application to the Stochastic Volatility Model.” Metrika 76 (8): 1031–81.
Epstein, Charles L. 2005. “How Well Does the Finite Fourier Transform Approximate the Fourier Transform?” Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences 58 (10): 1421–35.
M. Crouzeix and A.L. Mignot. Analyse numérique des équations différentielles. Collection mathématiques appliquées pour la maîtrise. Masson, 1992.
A.M. Quarteroni, J.F. Gerbeau, R. Sacco, and F. Saleri. Méthodes numériques pour le calcul scientifique : Programmes en MATLAB. Collection IRIS. Springer Paris, 2000.
A. Iserles and S. Norsett. On quadrature methods for highly oscillatory integrals and their implementation. BIT Numerical Mathematics, 44 :755–772, 2004.
A. Iserles and S. Norsett. Efficient quadrature of highly oscillatory integrals using derivatives. Proc. R. Soc. A, 461 :1383–1399, 2005.