Дискретное преобразование Фурье

Напоминание про ряд Фурье (непрерывный случай)

Рассмотрим разложение произвольной действительной функции в следующем виде

\[ x(t)=\sum_{i=1}^n A_i \cos \left(2 \pi \nu_i t+\phi_i\right) \]

Здесь мы определяем набор частот \(\left\{\nu_i\right\}\) некоторым способом, который мы обсудим позже. Значения амплитуд и сдвиг фаз оцениваются для каждого данного \(\nu_i\). Ряд Фурье описывает сигнал \(x(t)\) набором пар \(\left\{A_i^2, \nu_i\right\}_{i=1}^n\), где \(A_i^2\) представляет собой интенсивность \(i\)-й косинусоиды. Этот набор пар называется спектром.

Как вычислить \(A_i\)? Раскроем к общепринятому виду.

\[\begin{split} \begin{aligned} A_i \cos \left(2 \pi \nu_i t+\phi_i\right) & =A_i \cos \left(\phi_i\right) \cos \left(2 \pi \nu_i t\right)-A_i \sin \left(\phi_i\right) \sin \left(2 \pi \nu_i t\right)= \\ & =a_i \cos \left(2 \pi \nu_i t\right)+b_i \sin \left(2 \pi \nu_i t\right) \end{aligned} \end{split}\]

А тогда разложение,

\[ x(t)=\sum_{i=1}^n\left(a_i \cos \left(2 \pi \nu_i t\right)+b_i \sin \left(2 \pi \nu_i t\right)\right) = \sum_{i=1}^n\left(a_i - i \cdot b_i \right) \cdot e^{2 \pi \nu_i t} \]

где все коэффициенты считаются как

\[\begin{split} \begin{aligned} a_i & =\frac{2}{P} \int_P x(t) \cos \left(2 \pi \nu_i t\right) d t \\ b_i & =\frac{2}{P} \int_P x(t) \sin \left(2 \pi \nu_i t\right) d t \end{aligned} \end{split}\]

и

\[ A_i^2=a_i^2+b_i^2=\left(a_i - i \cdot b_i \right)^2 \]

Дискретный Фурье (DFT)

На компьютере нельзя задать континуальный набор точек - любой сигнал это дискретизованный конечный набор точек \(\left\{x_n\right\}_{n=0}^{N-1}\). Его вывод возьмём с википедии

Рассмотрим некоторый периодический сигнал \(x(t)\) с периодом, равным Т. Разложим его в ряд Фурье:

\[ x(t)=\sum_{k=-\infty}^{+\infty} c_k e^{i \omega_k t}, \quad \omega_k=\frac{2 \pi k}{T} . \]

Проведем дискретизацию сигнала так, чтобы на периоде было \(\mathrm{N}\) отсчетов. Дискретный сигнал представим в виде отсчетов: \(x_n=x\left(t_n\right)\), где \(t_n=\frac{n}{N} T\), тогда эти отсчеты через ряд Фурье запишутся следующим образом:

\[ x_n=\sum_{k=-\infty}^{+\infty} c_k e^{i \omega_k t_n}=\sum_{k=-\infty}^{+\infty} c_k e^{\frac{2 \pi i}{N} k n} \]

Используя соотношение \(e^{\frac{2 \pi i}{N}(k+m N) n}=e^{\frac{2 \pi i}{N} k n}\), получаем:

\[ x_n=\sum_{k=0}^{N-1} X_k e^{\frac{2 \pi i}{N} k n}, \quad \text { где } \quad X_k=\sum_{l=-\infty}^{+\infty} c_{k+l N} . \]

Таким образом мы получили обратное дискретное преобразование Фурье. Умножим теперь скалярно выражение для \(x_n\) на \(e^{-\frac{2 \pi i}{N} m n}\) и получим:

\[ \sum_{n=0}^{N-1} x_n e^{-\frac{2 \pi i}{N} m n}=\sum_{k=0}^{N-1} \sum_{n=0}^{N-1} X_k e^{\frac{2 \pi i}{N}(k-m) n}=\sum_{k=0}^{N-1} X_k \frac{1-e^{2 \pi i(k-m)}}{1-e^{\frac{2 \pi i(k-m)}{N}}}=\sum_{k=0}^{N-1} X_k N \delta_{k m} . \]

Здесь использованы: а) выражение для суммы конечного числа членов (экспонент) геометрической прогрессии, и б) выражение символа Кронекера как предела отношения функций Эйлера для комплексных чисел. Отсюда следует, что:

\[ X_k=\frac{1}{N} \sum_{n=0}^{N-1} x_n e^{-\frac{2 \pi i}{N} k n} . \]

Эта формула описывает прямое дискретное преобразование Фурье. В литературе принято писать множитель \(\frac{1}{N}\) в обратном преобразовании, и поэтому обычно пишут формулы преобразования в следующем виде:

\[ X_n=\sum_{m=0}^{N-1} x_m e^{-\frac{2 \pi i}{N} n \cdot m}, \quad x_m=\frac{1}{N} \sum_{n=0}^{N-1} X_n e^{\frac{2 \pi i}{N} n \cdot m} \]

Не столь важно, где этот коэффициент находится. Наиболее ествественным, вообще говоря, является \(\frac{1}{\sqrt {N}}\) в обоих выражениях.

Дискретный Фурье как оператор

Заметим, что дискретный Фурье можно задать в операторной форме:

\[ \vec{X}=M \cdot \vec{x} \]

с матрицей \(M\) заданной

\[ M_{k n}=e^{-i 2 \pi \cdot k \cdot n / N} \]

Итого,

\[ X_k = M_{kn} \cdot x_{n} \]

Для обратного преобразования Фурье можно применять алгоритм прямого преобразования Фурье — нужно лишь использовать \(\varepsilon^{-1}\) вместо \(\varepsilon\) (или применить операцию комплексного сопряжения вначале к входным данным, а затем к результату, полученному после прямого преобразования Фурье), и окончательный результат поделить на \(N\).

Аналогично можно определить \(m\)-мерный Фурье. Пример для Фурье преобразования матрицы:

\[ X_{kl} = M_{kn} \cdot M_{lm} \cdot x_{nm} \]

Связь оператора дискретного Фурье с быстрым умножением многочленов*

Дискретное преобразование Фурье является линейным преобразованием, которое переводит вектор временных отсчётов \(\vec x\) в вектор спектральных отсчётов той же длины. Таким образом преобразование может быть реализовано как умножение симметричной квадратной матрицы на вектор:

\[\vec X = \mathcal F \vec x,\]

матрица \(\mathcal F\) имеет вид:

\[\begin{split} \mathcal F = \frac{1}{\sqrt n}\begin{pmatrix} 1 & 1 & 1 & 1 & \ldots & 1 \\ 1 &\omega_n & \omega_n^2 & \omega_n^3 & \ldots & \omega_n^{n-1} \\ 1 &\omega_n^2 & \omega_n^4 & \omega_n^6 & \ldots & \omega_n^{2(n-1)} \\ 1 &\omega_n^3 & \omega_n^6 & \omega_n^9 & \ldots & \omega_n^{3(n-1)} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \omega_n^{3(n-1)} & \ldots &\omega_n^{(n-1)^2} \end{pmatrix} \quad \text{или} \quad \mathcal F_{jk} = \omega_n^{(j-1)(k-1)} \end{split}\]

где \(\omega_n = e^{-\frac{2\pi i}{n}}\)

Тут мы видим так называемую матрицу Вандермонда. Она тесно связана с полиномиальной интерполяцией.

А тогда дискретное преобразование Фурье вектора \((a_0,\,a_1,\,\dots,\,a_{n-1})\) может быть интерпретировано как вычисление значений многочлена

\[p(x)=a_0+a_1 x + \dots + a_{n-1}x^{n-1}\]

в следующих точках:

\[ x_0=1, \quad x_1=\omega_n, \quad x_2 = \omega_n^2, \quad …, \quad x_{n-1}=\omega_n^{n-1}. \]

Иными словами, Фурье-образ \((A_0,\,A_1,\,\dots,\,A_{n-1})\) вектора \((a_0,\,a_1,\,\dots,\,a_{n-1})\) можно страктовать как

\[ A_k = p(x_k) \]

Значения многочлена \(n\)-й степени в \(n+1\) точках однозначно определяют сам многочлен. Таким образом, если мы уже знаем преобразование Фурье \(\vec{A}\) (вычислили по формулам выше или через fft (формулы ниже)), мы имеем всю информацию для задания многочлена p(x).

Теорема. Пусть есть два многочлена каких-то своих степеней

\[p(x)=a_0+a_1 x + \dots + a_{n-1}x^{n-1}\]

\[q(x)=b_0+b_1 x + \dots + b_{m-1}x^{m-1}\]

Тогда для быстрого нахождения коэффициентов многочлена \((q\cdot p) \,(x)=c_0+c_1 x + \dots + c_{n+m-2}x^{n+m-2}\) достаточно применить обратное дискретное преобразование Фурье для вектора

\[ C_k = q(x_k) \cdot p(x_k), \quad x_k \stackrel{\operatorname{def}}{=} e^{-\frac{2\pi i \cdot k}{n+m-1}}, \quad k \in \overline{0, \ldots, n+m-2} \]

Примечание. Алгоритм выше эффективнее простого перемножения коэффициентов, если использовать fft (о нём ниже, работает за \(O(n log n)\)) для вычисления дискретного преобразования Фурье.

Следствие. Так как любое натуральное число представимо в виде многочлена от основания системы счисления:

\[N=a_{n-1} \dots a_1 a_0 \stackrel{\operatorname{def}}{=} a_0 + a_1 \cdot 10 + \dots + a_{n-1} \cdot 10^{n-1}\]

\[N(x)\stackrel{\operatorname{def}}{=} a_0 + a_1 \cdot x + \dots + a_{n-1} \cdot x^{n-1}\]

умножение двух чисел может быть в свою очередь сведено к умножению двух многочленов и нормализации результата. Работать такое умножение будет за \(O(N log N)\), где \(N\) - суммарное количество цифр в десятичной записи обоих чисел.