А что, если функции БПФ в NumPy и SciPy на самом деле не вычисляют преобразование Фурье, как вы думаете?
Делиться
Чтобы извлечь максимальную пользу из этой статьи, вам необходимо иметь базовые знания об интегралах , преобразовании Фурье и его свойствах.
Если нет, мы рекомендуем прочитать Раздел 1 , где эти идеи рассматриваются.
Если вы уже с ними знакомы, вы можете перейти сразу к Разделу 2 .В разделе 2.1 мы реализуем преобразование Фурье с использованием метода численных квадратур .
В разделе 2.2 мы реализуем его с помощью алгоритма быстрого преобразования Фурье (БПФ) и интерполяционной формулы Шеннона .
Идея написать эту работу пришла нам в голову на последнем курсе обучения в инженерной школе.
В рамках курса по калибровке стохастических процессов наш преподаватель хотел проверить наше понимание материала. Для этого нам было предложено выбрать научную работу, подробно изучить её и воспроизвести полученные в ней результаты.
Мы выбрали работу Эль Колея (2013) , в которой предлагается параметрический метод оценки параметров стохастической модели волатильности , такой как модель GARCH. Подход основан на минимизации контраста и деконволюции .
Чтобы воспроизвести результаты, нам потребовалось реализовать метод оптимизации , включающий вычисление преобразования Фурье f̂ функции fθ:
Для вычисления преобразования Фурье функции fθ Эль Колей (2013) использовал метод левой суммы Римана (также известный как метод квадратуры прямоугольника), реализованный в MATLAB. В статье предполагалось, что использование алгоритма быстрого преобразования Фурье (БПФ) может ускорить вычисления.
Мы решили воспроизвести результаты с помощью Python и реализовали преобразование Фурье, чтобы проверить, действительно ли оно увеличивает скорость вычислений. Реализация метода квадратур левого прямоугольника была простой, и поначалу мы думали, что библиотеки scipy.fft и numpy.fft позволят нам вычислить преобразование Фурье fθ напрямую с помощью алгоритма БПФ.
Однако мы быстро обнаружили, что эти функции вычисляют не преобразование Фурье непрерывной функции, а дискретное преобразование Фурье (ДПФ) конечной последовательности.
Рисунок ниже показывает, что мы наблюдали.
import numpy as np import matplotlib.pyplot as plt from scipy.fft import fft, fftfreq, fftshift # Определим функцию f(t) = exp(-pi * t^2) def f(t): return np.exp(-np.pi * t**2) # Параметры N = 1024 T = 1.0 / 64 t = np.linspace(-N/2*T, N/2*T, N, endpoint=False) y = f(t) # БПФ с помощью scipy yf_scipy = fftshift(fft(y)) * T xf = fftshift(fftfreq(N, T)) FT_exact = np.exp(-np.pi * xf**2) # БПФ с помощью numpy yf_numpy = np.fft.fftshift(np.fft.fft(y)) * T xf_numpy = np.fft.fftshift(np.fft.fftfreq(N, T)) # График с подграфиком subplot_mosaic fig, axs = plt.subplot_mosaic([[«scipy», «numpy»]], figsize=(7, 5), layout=»constrained», sharey=True) # БПФ Scipy axs[«scipy»].plot(xf, FT_exact, 'k-', linewidth=1.5, label='Точное ПФ') axs[«scipy»].plot(xf, np.real(yf_scipy), 'r—', linewidth=1, label='БПФ (scipy)') axs[«scipy»].set_xlim(-6, 6) axs[«scipy»].set_ylim(-1, 1) axs[«scipy»].set_title(«БПФ Scipy») axs[«scipy»].set_xlabel(«Частота») axs[«scipy»].set_ylabel(«Амплитуда») axs[«scipy»].legend() axs[«scipy»].grid(False) # БПФ NumPy axs[«numpy»].plot(xf_numpy, ПФ_exact, 'k-', ширина линии=1.5, метка='Точное ПФ') axs[«numpy»].plot(xf_numpy, np.real(yf_numpy), 'b—', ширина линии=1, метка='БПФ (numpy)') axs[«numpy»].set_xlim(-6, 6) axs[«numpy»].set_title(«БПФ NumPy») axs[«numpy»].set_xlabel(«Частота») axs[«numpy»].legend() axs[«numpy»].grid(False) plt.suptitle(«Сравнение реализаций БПФ и точного преобразования Фурье», fontsize=14) plt.show()
Этот опыт вдохновил нас на написание данной статьи, в которой мы объясним, как вычислить преобразование Фурье функции в Python, используя два подхода: метод левой суммы Римана и алгоритм быстрого преобразования Фурье (БПФ).
В литературе имеется множество работ, посвященных аппроксимации преобразования Фурье и его реализации с использованием численных методов.
Однако нам не удалось найти столь же ясного и полного источника, как Балак (2011) . В этой работе представлен квадратурный подход к вычислению преобразования Фурье и объясняется, как использовать алгоритм быстрого преобразования Фурье (БПФ) для эффективного выполнения дискретного преобразования Фурье.
Алгоритм БПФ может быть применен как для вычисления преобразования Фурье интегрируемой функции , так и для вычисления коэффициентов Фурье периодической функции .
1. Определение и свойства преобразования Фурье
Мы принимаем структуру, предложенную Балаком (2011), для определения преобразования Фурье функции f и представления его свойств.
Мы рассматриваем функции, принадлежащие пространству интегрируемых функций, обозначаемому L¹(ℝ, 𝕂) . Это пространство включает все функции f: ℝ → 𝕂 , где 𝕂 представляет собой либо действительные числа (ℝ), либо комплексные числа (ℂ).
Эти функции интегрируемы в смысле Лебега, что означает, что интеграл от их абсолютной величины конечен.
Итак, чтобы функция f принадлежала L¹(ℝ, 𝕂) , произведение f(t) · e – 2iπνt также должно быть интегрируемым для всех ν ∈ ℝ .
В этом случае преобразование Фурье функции f , обозначаемое f̂ (или иногда 𝓕(f) ), определяется для всех ν ∈ ℝ следующим образом:
Отметим, что преобразование Фурье f̂ функции f является комплекснозначной линейной функцией, зависящей от частоты ν .
Если f ∈ L¹(ℝ) — вещественная и чётная, то f̂ также вещественная и чётная. И наоборот, если f — вещественная и нечётная, то f̂ также чисто мнимая и нечётная.
Для некоторых функций преобразование Фурье можно вычислить аналитически. Например, для функции
f : t ∈ ℝ ↦ 𝟙 [−a⁄2, a⁄2] (t) ,
Преобразование Фурье имеет вид:
f̂(ν) = a sinc(a π ν)
где sinc(t) = sin(t)/t или t ∈ ℝ * , а sinc(0) = 0 .
Однако для многих функций преобразование Фурье невозможно вычислить аналитически. В таких случаях мы используем численные методы для его аппроксимации. Мы рассмотрим эти численные подходы в следующих разделах статьи.
2. Как численно аппроксимировать преобразование Фурье?
В своей статье Балак (2011) показывает, что вычисление преобразования Фурье функции включает в себя ее аппроксимацию следующим интегралом по интервалу [−T⁄2, T⁄2] :
Где T — достаточно большое число, чтобы интеграл сходился. Приближенное значение интеграла S̃(ν) можно вычислить с помощью квадратурных методов .
В следующем разделе мы аппроксимируем этот интеграл, используя метод левых квадратур (также известный как левая сумма Римана ).
2.1 Метод квадратур левых прямоугольников
Чтобы вычислить интеграл S̃(ν) с помощью метода левых квадратур прямоугольника , действуем следующим образом:
- Дискретизация интервала
Разобьем интервал [−T⁄2, T⁄2] на N равномерных подынтервалов длиной ht = T⁄N .
Точки дискретизации, соответствующие левым концам прямоугольников, определяются как:
tk = −T⁄2 + ht , для k = 0, 1, …, N−1.
- Приближение интеграла
Используя соотношение Шаля , аппроксимируем интеграл S̃(ν) следующим образом:
Поскольку tₖ₊₁ − tₖ = hₜ и tₖ = −T⁄2 + k·hₜ = T(k⁄N − ½) , выражение принимает вид:
Мы называем этот подход методом квадратуры левого прямоугольника, поскольку он использует левую конечную точку tₖ каждого подинтервала для аппроксимации значения f(t) в пределах этого интервала.
- Окончательная формула
Окончательное выражение для аппроксимации преобразования Фурье имеет вид:
2.1.1 Реализация метода квадратур левого прямоугольника на Python.
Функция tfquad ниже реализует метод квадратуры левого прямоугольника для вычисления преобразования Фурье функции f на заданной частоте nu.
import numpy as np def tfquad(f, nu, n, T): «»» Вычисляет преобразование Фурье функции f на частоте nu, используя левую квадратурную сумму Римана на интервале [-T/2, T/2]. Параметры: ———- f : вызываемая Функция для преобразования. Должна принимать в качестве входных данных массив NumPy. nu : число с плавающей точкой Частота, на которой вычисляется преобразование Фурье. n : целое число Количество квадратурных точек. T : число с плавающей точкой Ширина временного окна [-T/2, T/2]. Возвращает: ——- tfnu : комплексное Приближенное значение преобразования Фурье на частоте 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) возвращает префактор * np.sum(f(t_k) * weights)
Мы также можем использовать функцию quad из SciPy для определения преобразования Фурье функции f на заданной частоте nu. Функция tf_integral, представленная ниже, реализует этот подход. Она использует численное интегрирование для вычисления преобразования Фурье функции f в интервале [-T/2, T/2].
из scipy.integrate import quad def tf_integral(f, nu, T): «»»Вычислить преобразование Фурье функции 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
Выведя формулу для метода квадратуры левого прямоугольника, мы теперь можем реализовать ее на Python, чтобы посмотреть, как она работает на практике.
Для этого рассмотрим простой пример, в котором функция fff определена как индикаторная функция на интервале [−1, 1] :
Для этой функции преобразование Фурье можно вычислить аналитически, что позволяет сравнить численное приближение с точным результатом .
Следующий скрипт Python реализует как аналитическое преобразование Фурье, так и его численные приближения с помощью:
- метод квадратуры левого прямоугольника и
- Интеграционная функция SciPy для справки.
import numpy as np import matplotlib.pyplot as plt # —— Определения функций —— def f(t): «»»Индикаторная функция на [-1, 1].»»» return np.where(np.abs(t) <= 1, 1.0, 0.0) def exact_fourier_transform(nu): """Аналитическое преобразование Фурье индикаторной функции на [-1, 1].""" # f̂(ν) = ∫_{-1}^{1} e^{-2πiνt} dt = 2 * sinc(2ν) return 2 * np.sinc(2 * nu) # ----- Вычисление ----- T = 2.0 n = 32 nu_vals = np.linspace(-6, 6, 500) exact_vals = exact_fourier_transform(nu_vals) tfquad_vals = np.array([tfquad(f, nu, n, T) for nu in nu_vals]) # Вычислить аппроксимацию с помощью интеграла scipy tf_integral_vals = np.array([tf_integral(f, nu, T) for nu in nu_vals]) # ----- Построение графика ----- fig, axs = plt.subplot_mosaic([["tfquad", "quad"]], figsize=(7.24, 4.07), dpi=100, layout="constrained") # Построить график с помощью реализации tfquad axs["tfquad"].plot(nu_vals, np.real(exact_vals), 'b-', linewidth=2, label=r'$hat{f}$ (exact)') axs["tfquad"].plot(nu_vals, np.real(tfquad_vals), 'r--', linewidth=1.5, label=r'approximation $hat{S}_n$') axs["tfquad"].set_title("TF с tfquad (прямоугольники)") axs["tfquad"].set_xlabel(r'$nu$') axs["tfquad"].grid(False) axs["tfquad"].set_ylim(-0.5, 2.1) # Построить график с помощью scipy.integrate.quad axs["quad"].plot(nu_vals, np.real(exact_vals), 'b-', linewidth=2, label=r'$hat{f}$ (quad)') axs["quad"].plot(nu_vals, np.real(tf_integral_vals), 'r--', linewidth=1.5, label=r'approximation $hat{S}_n$') axs["quad"].set_title("TF avec scipy.integrate.quad") axs["quad"].set_xlabel(r'$nu$') axs["quad"].set_ylabel('Amplitude') axs["quad"].grid(False) axs["quad"].set_ylim(-0.5, 2.1) # --- Глобальная легенда под графиками --- # Берем маркеры только из одного подграфика (предполагается, что метки согласованы) handles, labels = axs["quad"].get_legend_handles_labels() fig.legend(handles, labels, loc='lower center', bbox_to_anchor=(0.5, -0.05), ncol=3, frameon=False) plt.suptitle("Сравнение реализаций БПФ и точного преобразования Фурье", fontsize=14) plt.show()
2.1.2 Характеристика аппроксимации методом левых квадратур
- Аппроксимация преобразования Фурье с использованием метода квадратур левого прямоугольника по своей природе имеет колебательный характер .
Как отметил Балак (2011) , преобразование Фурье f̂ функции f по своей природе является осцилляционным. Это поведение обусловлено комплексным экспоненциальным членом e⁻²πiνt в интегральном определении преобразования.
Для иллюстрации этого поведения на рисунке ниже показана функция
ж : т € ℝ ↦ et^2 € ℝ
вместе с действительной и мнимой частями его преобразования Фурье
f̂ : ν ∈ ℝ ↦ f̂(ν) ∈ ℂ , оцененный как ν = 5⁄2 .
Хотя f(t) является гладкой, мы можем наблюдать сильные колебания f̂(ν) , которые выявляют влияние экспоненциального члена e-2 πiνt в интеграле Фурье.
import numpy as np import matplotlib.pyplot as plt nu = 5 / 2 t1 = np.linspace(-8, 8, 1000) t2 = np.linspace(-4, 4, 1000) f = лямбда t: np.exp(-t**2) phi = лямбда t: f(t) * np.exp(-2j * np.pi * nu * t) f_vals = f(t1) phi_vals = phi(t2) # Построить график fig, axs = plt.subplots(1, 2, figsize=(10, 4)) axs[0].plot(t1, f_vals, 'k', linewidth=2) axs[0].set_xlim(-8, 8) axs[0].set_ylim(0, 1) axs[0].set_title(r»$f(t) = e^{-t^2}$») axs[0].grid(True) axs[1].plot(t2, np.real(phi_vals), 'b', label=r»$Re(phi)$», linewidth=2) axs[1].plot(t2, np.imag(phi_vals), 'r', label=r»$Im(phi)$», linewidth=2) axs[1].set_xlim(-4, 4) axs[1].set_ylim(-1, 1) axs[1].set_title(r»$phi(t) = f(t)e^{-2ipinu t}$, $nu=5/2$») axs[1].legend() axs[1].grid(True) plt.tight_layout() plt.show()
- Приближение, полученное с помощью метода квадратур левого прямоугольника, имеет периодический характер.
Мы видим, что даже если функция f не является периодической , ее приближение с помощью преобразования Фурье f̂ оказывается периодическим .
На самом деле функция Ŝₙ , полученная методом квадратур, является периодической с периодом, определяемым по формуле:
Эта периодичность Ŝₙ означает, что невозможно вычислить преобразование Фурье для всех частот ν ∈ ℝ с помощью квадратурного метода, когда параметры T и n фиксированы.
На самом деле, становится невозможным точно вычислить f̂(ν) , когда
|ν| ≥ νмакс,
где v = n/T представляет собой максимальную частоту , которую можно разрешить из-за периодической природы Ŝₙ .
В результате на практике для вычисления преобразования Фурье для больших частот необходимо увеличить либо временное окно (T) , либо количество точек (n) .
Кроме того, оценивая погрешность аппроксимации f̂(ν) с использованием метода левых квадратур , можно показать, что аппроксимация надежна на частоте ν , когда выполняется следующее условие:
|ν| ≪ n/T или, что эквивалентно, (|ν|T)/n ≪ 1.
Согласно Эпштейну (2005) , при использовании алгоритма быстрого преобразования Фурье (БПФ) можно точно вычислить преобразование Фурье функции f для всех частот ν ∈ ℝ , даже когда (|ν|T)⁄n близко к 1, при условии, что f является кусочно-непрерывной и имеет компактный носитель .
2.2 Вычисление преобразования Фурье на частоте ν с использованием алгоритма БПФ
В этом разделе мы обозначаем через Ŝₙ(ν) приближение преобразования Фурье f̂(ν) функции f в точке ν ∈ [− νmax /2, νmax /2] , где νmax = n/T , то есть
Теперь я представлю алгоритм преобразования Фурье , используемый для аппроксимации f̂(ν) .
В этой статье я не буду вдаваться в подробности алгоритма быстрого преобразования Фурье (БПФ) .
Для упрощенного объяснения я отсылаю к работе Балака (2011) , а для более технического и всестороннего изучения — к оригинальной работе Кули и Тьюки (1965) .
Важно понимать, что использование алгоритма БПФ для аппроксимации преобразования Фурье функции f основано на результате, полученном Эпштейном (2005) . Этот результат утверждает, что при вычислении Ŝₙ(ν) на частотах vj = j/T , для j = 0, 1, …, n − 1 , он обеспечивает хорошее приближение непрерывного преобразования Фурье f̂(ν) .
Более того, известно, что Ŝₙ является периодическим . Эта периодичность придаёт индексам j ∈ {0, 1, …, n − 1} и k ∈ {−n/2, −n/2 + 1, …, −1} симметричную роль.
Фактически, значения преобразования Фурье функции f на интервале [− νmax /2, νmax /2] могут быть получены из значений Ŝₙ в точках νj = j/T для j = 0, 1, …, n − 1 следующим образом:
где мы использовали соотношение:
Это соотношение показывает, что преобразование Фурье Ŝₙ(j/T) можно вычислить для j = −n⁄2, …, n⁄2 − 1.
Более того, когда n — степень двойки, вычисления значительно ускоряются (см. Balac, 2011). Этот процесс известен как быстрое преобразование Фурье (БПФ) .
Подводя итог, я показал, что преобразование Фурье функции f можно аппроксимировать на интервале [−T⁄2, T⁄2] на частотах vj = j/T , для j = −n/2, …, n/2 − 1 , где n = 2ᵐ для некоторого целого числа m ≥ 0 , применяя алгоритм быстрого преобразования Фурье (БПФ) следующим образом:
- Шаг 1: Построить конечную последовательность F значений
f((2k − n)T/(2n)) , для k = 0, 1, …, n − 1 . - Шаг 2: Вычислить дискретное преобразование Фурье (ДПФ) последовательности F , используя алгоритм БПФ, заданный формулой:
- Шаг 3: Переиндексируйте и симметризуйте значения, чтобы охватить диапазон j = −n/2, …, −1 .
- Шаг 4: Умножьте каждое значение в массиве на (T/n)(−1) j-1, где j ∈ {1, …, n} .
Этот процесс дает массив, представляющий значения преобразования Фурье f̂(νⱼ) , где νj = j/T для j = −n⁄2, …, n⁄2 − 1 .
Функция Python tffft, представленная ниже, реализует эти шаги для вычисления преобразования Фурье заданной функции.
import numpy as np из scipy.fft import fft, fftshift def tffft(f, T, n): «»» Вычислить преобразование Фурье, подходящую функцию для поддержки в [-T/2, T/2], с использованием алгоритма FFT. Параметры ———- f : вызываемая функция преобразователя (doit être векторизуемый с числом). T : float Largeur de la fenêtre temporelle (intervalle [-T/2, T/2]). n : int Число точек дискретизации (doit être une puissance de 2 for FFT Retours ——- tf : np.ndarray Приблизительные значения). преобразование Фурье на дискретных частотах. freq_nu : np.ndarray Дискретные соответствующие частоты (de -n/(2T) à (n/2 — 1)/T). «»» h = T / nt = -0.5 * T + np.arange(n) * h # noeuds temporels F = f(t) # échantillonnage de f tf = h * (-1) ** np.arange(n) * fftshift(fft(F)) # TF Approimée freq_nu = -n / (2 * T) + np.arange(n) / T # частоты ν_j = j/T return tf, freq_nu, t
Следующая программа иллюстрирует, как вычислить преобразование Фурье гауссовой функции f(t) = e-10t^2 на интервале [-10, 10] с использованием алгоритма быстрого преобразования Фурье (БПФ).
# Параметры a = 10 f = лямбда t: np.exp(-a * t**2) T = 10 n = 2**8 # 256 # Вычислить преобразование Фурье с помощью БПФ tf, nu, t = tffft(f, T, n) # Построить график fig, axs = plt.subplots(1, 2, figsize=(7.24, 4.07), dpi=100) axs[0].plot(t, f(t), '-g', linewidth=3) axs[0].set_xlabel(«время») axs[0].set_title(«Рассматриваемая функция») axs[0].set_xlim(-6, 6) axs[0].set_ylim(-0.5, 1.1) axs[0].grid(True) axs[1].plot(nu, np.abs(tf), '-b', linewidth=3) axs[1].set_xlabel(«частота») axs[1].set_title(«Преобразование Фурье с использованием БПФ») axs[1].set_xlim(-15, 15) axs[1].set_ylim(-0.5, 1) axs[1].grid(True) plt.tight_layout() plt.show()
Метод, который я только что представил, позволяет вычислить и визуализировать преобразование Фурье функции f в дискретных точках vj = j/T , для j = −n/2, …, n/2 − 1 , где n — степень двойки.
Эти точки лежат в интервале [−n/(2T), n/(2T)] .
Однако этот метод не позволяет оценить преобразование Фурье на произвольных частотах ν, которые не имеют вида vj = j/T .
Для вычисления преобразования Фурье функции f на частоте ν , которая не совпадает ни с одной из выборочных частот νⱼ , можно использовать методы интерполяции, такие как линейная , полиномиальная или сплайн-интерполяция .
В этой статье мы используем подход, предложенный Балаком (2011) , который основан на интерполяционной теореме Шеннона для вычисления преобразования Фурье функции f на любой частоте ν .
Использование интерполяционной теоремы Шеннона для вычисления преобразования Фурье функции f в точке ν
О чем говорит нам теорема Шеннона?
В нем говорится, что для функции g с ограниченной полосой пропускания , то есть функции, преобразование Фурье которой ĝ равно нулю вне интервала [−B⁄2, B⁄2] , функция g может быть восстановлена по ее отсчетам gₖ = g(k⁄B) для k ∈ ℤ .
Если мы допустим, что ν c будет наименьшим положительным числом, таким, что ĝ(ν) равно нулю вне интервала [−2π ν c, 2π ν c] , то применима интерполяционная формула Шеннона .
Для любого t ∈ ℝ и любого положительного действительного числа α ≥ 1/(2 ν c) имеем:
где sinc(x) = sin(x)/x — функция sinc (кардинальный синус).
Балак (2011) показывает, что когда функция f имеет ограниченный носитель в интервале [−T⁄2, T⁄2] , интерполяционную теорему Шеннона можно использовать для вычисления ее преобразования Фурье f̂(ν) для любого ν ∈ ℝ .
Это делается с помощью значений дискретного преобразования Фурье Ŝₙ(νⱼ) для j = −n/2, …, (n/2 − 1 ).
Установив α = 1/T , Балак выводит следующую интерполяционную формулу Шеннона , верную для всех ν ∈ ℝ :
Следующая программа иллюстрирует, как использовать интерполяционную теорему Шеннона для вычисления преобразования Фурье функции f на заданной частоте ν .
Для вычисления преобразования Фурье на любой частоте ν можно воспользоваться интерполяционной теоремой Шеннона .
Идея состоит в том, чтобы оценить значение преобразования Фурье по его дискретным выборкам, полученным с помощью БПФ.
Функция ниже, shannon, реализует этот метод интерполяции в Python.
def shannon(tf, nu, T): «»» Аппроксимирует значение преобразования Фурье функции f на частоте 'nu', используя его дискретные значения, вычисленные из БПФ. Параметры: — tf : массив numpy, дискретные значения преобразования Фурье (центрированные с помощью fftshift) на частотах j/T для j = -n/2, …, n/2 — 1 — nu : число с плавающей точкой, частота, на которой аппроксимируется преобразование Фурье — T : число с плавающей точкой, ширина временного окна, используемого для БПФ Возвращает: — tfnu : аппроксимация преобразования Фурье на частоте 'nu' «»» n = len(tf) tfnu = 0.0 для j в диапазоне(n): k = j — n // 2 # соответствует индексу j в {-n/2, …, n/2 — 1} tfnu += tf[j] * np.sinc(T * nu — k) # np.sinc(x) = sin(pi x)/(pi x) en numpy return tfnu
Следующая функция, fourier_at_nu, объединяет два шага.
Сначала он вычисляет дискретное преобразование Фурье функции с использованием БПФ (tffft), затем применяет интерполяцию Шеннона для оценки преобразования Фурье на определенной частоте ν .
Это позволяет оценивать преобразование в любой произвольной точке ν , а не только на частотах дискретизации БПФ.
Теперь мы можем протестировать нашу реализацию на простом примере.
Здесь мы определяем функцию f(t) = ea|t| , преобразование Фурье которой известно аналитически.
Это позволяет нам сравнить точное преобразование Фурье с приближением , полученным с помощью нашего метода.
a = 0,5 f = лямбда t: np.exp(-a * np.abs(t)) # Функция преобразования fhat_exact = лямбда nu: (2 * a) / (a**2 + 4 * np.pi**2 * nu**2) # Точное преобразование Фурье T = 40 # Временное окно n = 2**10 # Количество точек дискретизации
Мы оцениваем преобразование Фурье на двух частотах: ν = 3/T и ν = π/T .
Для каждого случая мы выводим как точное значение , так и приближение, полученное с помощью нашей функции fourier_at_nu.
Это сравнение помогает проверить точность интерполяции Шеннона.
# Вычислить для nu = 3/T nu = 3 / T exact_value = fhat_exact(nu) approx_value = fourier_at_nu(f, T, n, nu) print(f»Точное значение при nu={nu}: {exact_value}») print(f»Приближение при nu={nu}: {np.real(approx_value)}») # Вычислить для nu = pi/T nu = np.pi / T exact_value = fhat_exact(nu) approx_value = fourier_at_nu(f, T, n, nu) print(f»Точное значение при nu={nu}: {exact_value}») print(f»Приближение при nu={nu}: {np.real(approx_value)}») Точное значение при ν = 3/T: 2.118347413776218 Аппроксимация при ν = 3/T: (2,1185707502943534) Точное значение при ν = π/T: 2,0262491352594427 Аппроксимация при ν = π/T: (2,0264201680784835)
Стоит отметить, что, хотя частота 3⁄T (значение преобразования Фурье которой появляется в выходных данных БПФ) близка к π⁄T , соответствующие значения преобразования Фурье f на этих двух частотах существенно различаются.
Это показывает, что интерполяционная формула Шеннона может быть очень полезна, когда желаемые значения преобразования Фурье не включены напрямую в результаты, полученные алгоритмом БПФ.
Заключение
В этой статье мы рассмотрели два метода аппроксимации преобразования Фурье функции.
Первым был метод численных квадратур (правило левого прямоугольника), а вторым — алгоритм быстрого преобразования Фурье (БПФ) .
Мы показали, что в отличие от встроенных функций БПФ в SciPy или NumPy, БПФ можно адаптировать для вычисления преобразования Фурье непрерывной функции с помощью правильной формулировки.
Наша реализация основана на работе Балака (2013) , который продемонстрировал, как воспроизвести вычисление БПФ на Python.
Мы также ввели интерполяционную теорему Шеннона , которая позволяет нам оценить преобразование Фурье любой функции на ℝ на произвольных частотах.
При необходимости эту интерполяцию можно заменить более традиционными подходами, такими как линейная или сплайновая интерполяция .
Наконец, важно отметить, что когда требуется преобразование Фурье на одной частоте , часто более эффективно вычислить его напрямую, используя метод численных квадратур .
Такие методы хорошо подходят для обработки колебательной природы интеграла Фурье и могут быть точнее, чем применение БПФ с последующей интерполяцией.
Ссылки
Балак, Стефан. 2011. «Преобразование Фурье под углом Du Calcul Numérique».
Кули, Джеймс В. и Джон В. Тьюки. 1965. «Алгоритм для машинного вычисления комплексных рядов Фурье». Математика вычислений 19 (90): 297–301.
Эль Колей, Салима. 2013. «Параметрическая оценка скрытой стохастической модели с помощью минимизации контраста и деконволюции: применение к модели стохастической волатильности». Метрика 76 (8): 1031–81.
Эпштейн, Чарльз Л. 2005. «Насколько хорошо конечное преобразование Фурье аппроксимирует преобразование Фурье?» Сообщения по чистой и прикладной математике: журнал, выпущенный Институтом математических наук Куранта 58 (10): 1421–35.
М. Крузе и А. Л. Миньо. Проанализируйте числовые значения различных уравнений. Коллекция математических аппликаций
pour la maîtrise. Массон, 1992.
А. М. Квартерони, Ж. Ф. Жербо, Р. Сакко и Ф. Салери. Числовые методы для научных расчетов: Программы
в МАТЛАБ. Коллекция ИРИС. Спрингер Париж, 2000.
А. Изерлес и С. Норсетт. О квадратурных методах для сильно осциллирующих интегралов и их реализации. BIT
Вычислительная математика, 44:755–772, 2004.
А. Изерлес и С. Норсетт. Эффективная квадратура сильно осциллирующих интегралов с использованием производных. Proc. R. Soc. A,
461 :1383–1399, 2005.
Кредиты изображений
Все изображения и визуализации в этой статье были созданы автором с использованием Python (pandas, matplotlib, seaborn и plotly) и Excel, если не указано иное.
Отказ от ответственности
Я пишу, чтобы учиться, поэтому ошибки — это нормально, хотя я и стараюсь изо всех сил. Пожалуйста, дайте мне знать, если заметите их. Я также открыт для любых предложений по новым темам!
Источник: towardsdatascience.com
