УказательРазделыОбозначенияАвторО проекте


§

Вспомогательная страница к разделу ИНТЕРПОЛЯЦИЯ. Cущественно используются материалы разделов КОМПЛЕКСНЫЕ ЧИСЛА и ПОЛИНОМ ОДНОЙ ПЕРЕМЕННОЙ.


Тригонометрическая интерполяция

Этот раздел посвящен периодическим функциям. Они рассматриваются для случая периода равного 2\pi_{}, к которому можно свести общий случай соответствующей заменой переменной («растяжением» шкалы).

Тригонометрический полином

Тригонометрическим полиномом от переменной x_{} называется выражение вида

\begin{matrix} g_n(x)=a_0 & + & a_1 \cos x + a_2 \cos 2x+\dots + a_n \cos nx + \\ \ &+&b_1 \sin x + b_2 \sin 2x+\dots + b_n \sin nx = \end{matrix}
=a_0 + \sum_{k=1}^n (a_k \cos\, kx+ b_k \sin \, kx) \ .

при фиксированных постоянных a_0,a_1,a_2,\dots,a_n,b_1,\dots,b_n — вещественных или комплексных. Если a_{n} или b_n отличны от нуля, то говорят, что полином g_n(x) имеет порядок n_{}. В этом случае числа a_0,\{a_j,b_j\}_{j=1}^n называются коэффициентами тригонометрического полинома, а числа a_{n} и b_nстаршими коэффициентами. Если все числа \{b_j\}_{j=1}^n равны нулю, то о полиноме g_n(x) говорят как о полиноме по косинусам, если все числа \{a_j\}_{j=0}^n равны нулю — то полином называют полиномом по синусам.

§

Чтобы отличать тригонометрический полином от «обычного» a_0+a_1x+\dots + a_nx^n мы будем называть последний степенным или алгебраическим.

П

Пример. Функции \cos^n x и \sin^n x являются тригонометрическими полиномами n_{}-го порядка; первый из них является полиномом по косинусам:

\cos^3 x= \frac{3}{4}\cos x + \frac{1}{4}\cos 3\,x \ ;

второй же — в зависимости от четности n_{} — либо полиномом по косинусам, либо — по синусам:

\sin^4 x=\frac{3}{8}-\frac{1}{2}\cos 2\,x+\frac{1}{8}\cos 4\, x \ , \ \sin^5 x=\frac{5}{8}\sin x-\frac{5}{16} \sin 3\,x+\frac{1}{16}\sin 5\,x \ .

Для общего случая доказано ☞ ЗДЕСЬ.

Т

Теорема. Произвольный тригонометрический полином n_{}-го порядка может быть представлен в виде

g_n(x) \equiv e^{-\mathbf i n x} G_{2n} (e^{\mathbf i x}) \ ,

где

G_{2n}(z)=\sum_{k=0}^{2n} u_kz^k

— некоторый алгебраический полином степени 2n_{}.

Доказательство. Перейдем к комплексной переменной: в предположении, что x_{} переменная вещественная , обозначим

z= \cos x + \mathbf i \sin x \ ,

Тогда, на основании свойств операций над комплексными числами (в том числе формулы Муавра ), имеем

z^{-1} = \cos x - \mathbf i \sin x \ , z^j= \cos j\,x + \mathbf i \sin j\,x \ , z^{-j}= \cos j\,x - \mathbf i \sin j\,x \quad npu \ \in\{1,\dots,n\} \ ,

и, следовательно,

\cos j\,x = \frac{1}{2}\left(z^j+z^{-j}\right),\quad \sin j\, x=\frac{1}{2\mathbf i} \left(z^j-z^{-j}\right) \ .

Подставим эти выражения в g_n(x) и домножим получившееся на z^n:

G_{2n}(z)=\underbrace{\frac{a_n+\mathbf i b_n}{2}}_{u_0}+\underbrace{\frac{a_{n-1}+\mathbf i b_{n-1}}{2}}_{u_1}z+\dots+\underbrace{\frac{a_1+\mathbf i b_1}{2}}_{u_{n-1}}z^{n-1}+ \underbrace{a_0}_{u_n}z^n+\underbrace{\frac{a_1-\mathbf i b_1}{2}}_{u_{n+1}}z^{n+1}+\dots+\underbrace{\frac{a_{n-1}-\mathbf i b_{n-1}}{2}}_{u_{2n-1}}z^{2n-1} + \underbrace{\frac{a_n-\mathbf i b_n}{2}}_{u_{2n}}z^{2n} \ .

Если переменная x предполагается комплексной, то придется использовать экспоненту от комплексного числа: z=e^{\mathbf i x} и представления через нее \sin x и \cos x (см. ☞ЗДЕСЬ ).

=>

Если все коэффициенты a_0,\{a_j,b_j\}_{j=1}^n тригонометрического полинома g_n(x) вещественны, то полином G_{2n}(z) из теоремы удовлетворяет тождеству:

G_{2n}(z) \equiv z^{2n} \overline{G_{2n}}(z^{-1}) \ .

Иными словами, для коэффициентов такого полинома

G_{2n}(z)=u_0+u_1z+u_2z^2+\dots+u_{2n-1}z^{2n-1} + u_{2n}z^{2n}

справедлива «симметрия с сопряжением» набора его коэффициентов (u_0,u_1,\dots,u_{2n-1},u_{2n}) относительно середины этого набора:

u_0=\overline{u_{2n}},u_1=\overline{u_{2n-1}},\dots, u_j=\overline{u_{2n-j}},\dots

Корнем или нулем тригонометрического полинома g_n(x) называется комплексное число \lambda такое, что g_n(\lambda)=0. Договариваются не различать между собой корни различающиеся на целое кратное 2\pi_{} (т. е. сравнимые по модулю 2\pi_{}). Определение кратного корня тригонометрического полинома вводится аналогично случаю полинома алгебраического.

Следующая теорема является аналогом основной теоремы высшей алгебры для степенных полиномов.

Т

Теорема. Тригонометрический полином n_{}-го порядка, у которого старшие коэффициенты a_{n} и b_n удовлетворяют условиям a_n+ \mathbf i b_n \ne 0, a_n- \mathbf i b_n \ne 0, имеет в точности 2n комплексных корней с учетом их кратностей.

ДоказательствоЗДЕСЬ.

=>

Если числа \{\lambda_1,\dots, \lambda_{2n} \} при \{ \mathfrak{Re} (\lambda_k)\}_{k=1}^{2n} \subset [0,2\pi[ — корни тригонометрического полинома g_n(x), у которого a_n+ \mathbf i b_n \ne 0, a_n- \mathbf i b_n \ne 0, то этот полином можно представить в виде

g_n(x) \equiv A \prod_{k=1}^{2n} \sin \frac{x-\lambda_k}{2} \ npu \ A=\pm 2^{2n-1}\sqrt{a_n^2+b_n^2} \ .

=>

Имеет место тождество:

\sin \, nx \equiv 2^{n-1} \prod_{j=0}^{n-1} \sin\left(x-\frac{\pi j}{n} \right) \equiv 2^{n-1} \sin \, x \sin \left(x-\frac{\pi}{n}\right) \sin \left(x-\frac{2\pi}{n} \right)\times \dots \times \sin \left(x-\frac{(n-1)\pi}{n} \right) \ .

Доказательство обоих следствий ☞ ЗДЕСЬ.

Тригонометрическая интерполяция

Задача. Построить тригонометрический полином порядка n_{}, принимающий значения по следующей таблице

\begin{array}{c|ccccc} x & x_1 & x_2 & \dots & x_{2\,n+1} \\ \hline y & y_1 & y_2 &\dots & y_{2\,n+1} \end{array}

т.е. g_n(x_j) = y_j при j\in \{1,\dots,2\,n+1\}. Здесь узлы интерполяции \{x_j\}_{j=1}^{2\,n+1} предполагаются вещественными, различными и расположенными в интервале [0, 2\pi [.

Т

Теорема. Функция

\begin{matrix} g_n(x)&=&y_1\frac{\displaystyle \sin\frac{x-x_2}{2}\times\dots \times\sin \frac{x-x_{2n+1}}{2}}{\displaystyle \sin \frac{x_1-x_2}{2}\times \dots \times \sin \frac{x_1-x_{2n+1}}{2}} + \\ &+&y_2\frac{\displaystyle \sin\frac{x-x_1}{2}\sin\frac{x-x_3}{2}\times \dots \times \sin\frac{x-x_{2n+1}}{2}}{\displaystyle \sin\frac{x_2-x_1}{2}\sin\frac{x_2-x_3}{2}\times \dots \times \sin\frac{x_2-x_{2n+1}}{2}} + \\ & + & \dots + \\ &+& y_{2n+1}\frac{\displaystyle \sin \frac{x-x_1}{2} \times \dots \times \sin \frac{x-x_{2n}}{2}}{ \displaystyle \sin \frac{x_{2n+1}-x_1}{2} \times \dots \times \sin \frac{x_{2n+1}-x_{2n}}{2}} \end{matrix}

удовлетворяет условиям g_n(x_j)=y_j при j\in \{1,\dots,2n+1\}. Эта функция является тригонометрическим полиномом порядка не выше n_{}. При таком ограничении на порядок, полином, удовлетворяющий заданным интерполяционным условиям, определяется единственным образом.

Доказательство факта, что функция g_{n}(x) удовлетворяет интерполяционной таблице очевидно: имеем полный аналог интерполяционного степенного полинома в форме Лагранжа. Образно говоря, каждое слагаемое y_jW_j(x) в этой сумме «отвечает» исключительно только за свой узел интерполяции (т.е. обеспечивает в нем нужное значение) — и при этом «не портит» остальные узлы (обращается в них в нуль).

Теперь нужно показать, что функция g_{n}(x) является именно тригонометрическим полиномом, т.е. допускает представление вида \displaystyle a_0 + \sum_{k=1}^n (a_k \cos \, kx + b_k \sin \, kx). Это достигается путем преобразования произведений синусов, стоящих в числителях слагаемых. В каждом таком произведении содержится четное число сомножителей, объединяем их попарно и используем формулы приведения:

\sin \frac{x-x_1}{2} \sin \frac{x-x_2}{2} \equiv \frac{1}{2} \left( \cos \frac{x_2-x_1}{2} - \cos(x-\frac{x_1+x_2}{2} \right)\equiv \frac{1}{2} \cos \frac{x_2-x_1}{2} - \frac{1}{2} \cos \frac{x_1+x_2}{2} \cos\, x - \frac{1}{2} \sin \frac{x_1+x_2}{2} \sin \,x \ .

В результате, мы избавляемся от половинного аргумента под синусами. Дальнейшие преобразования аналогичны. Окончательная формула ☞ ЗДЕСЬ.

Единственность этого полинома при условии ограничения на порядок следует из теоремы предыдущего пункта: если h_n(x) — другой полином, удовлетворяющий условию теоремы, то полином g_n(x)-h_n(x) порядка \le 2\,n имеет \ge (2n+1)-го корней на интервале [0,2\, \pi[, что невозможно.

Итак, теоретическая база решения задачи тригонометрической интерполяции построена. Осталась практическая часть: интересует представление полинома из теоремы в каноническом виде

g_n(x)\equiv a_0 + \sum_{k=1}^n (a_k \cos\, kx+ b_k \sin \, kx) \ ,

т.е. его коэффициенты. Стандартный способ их поиска — метод неопределенных коэффициентов — приводит к системе из (2n+1)-го линейного уравнения

\left(\begin{array}{llllllll} 1 & \cos x_1 & \sin x_1 & \cos \, 2\, x_1 & \sin \, 2\, x_1 & \dots & \cos \, n\, x_1 & \sin \, n\, x_1 \\ 1 & \cos x_2 & \sin x_2 & \cos \, 2\, x_2 & \sin \, 2\, x_2 & \dots & \cos \, n\, x_2 & \sin \, n\, x_2 \\ \dots & & & & & & & \dots \\ 1 & \cos x_{2n+1} & \sin x_{2n+1} & \cos \, 2\, x_{2n+1} & \sin \, 2\, x_{2n+1} & \dots & \cos \, n\, x_{2n+1} & \sin \, n\, x_{2n+1} \end{array} \right) \left(\begin{array}{l} a_0 \\ a_1 \\ \vdots \\ b_n \end{array} \right)= \left(\begin{array}{l} y_1 \\ y_2 \\ \vdots \\ y_{2n+1} \end{array} \right) \ .

относительно (2n+1)-го коэффициента a_0,a_1,b_1,a_2,b_2,\dots,a_n,b_n. Согласно теореме, эта система должна иметь решение относительно столбца коэффициентов, причем решение единственное. Из этого факта, в частности, следует, что определитель матрицы этой системы должнен быть отличен от нуля (см. ☞ формулы Крамера ). Так оно и есть.

Т

Теорема. Определитель матрицы системы равен

2^{2n^2} \prod_{0\le k < j \le 2n} \sin \frac{x_k-x_j}{2} \ .

§

Очень похоже на выражение для определителя Вандермонда (к которому, собственно, и сводится вычисление).

Теперь обсудим возможные упрощения представления тригонометрического полинома из теоремы, которые получатся за счет удачного выбора узлов интерполяции.

Т

Теорема. Если в качестве узлов интерполяции взять равноотстоящие

x_j=\frac{2(j-1)\pi}{2n+1} \quad npu \quad j\in \{1,2,\dots,2n+1\} \ ,

то коэффициенты интерполяционного тригонометрического полинома находятся по формулам

a_0=\frac{1}{2n+1}\sum_{j=1}^{2n+1}y_j,
a_k= \frac{2}{2n+1} \sum_{j=1}^{2n+1} y_j \cos k\,x_j ,\quad b_k= \frac{2}{2n+1} \sum_{j=1}^{2n+1} y_j \sin k\,x_j \quad npu \ k\in \{1,\dots,n \} .

ДоказательствоЗДЕСЬ.

Формулы последней теоремы остаются справедливыми и в случае «сдвига» системы узлов интерполяции на величину, кратную \pi_{}, например, при выборе системы узлов

x_j=\frac{2\pi j}{2n+1} \quad npu \quad j\in \{-n,-n+1,\dots,-1,0,1,2,\dots,n\}

на интервале [-\pi,\pi[:

a_0=\frac{1}{2n+1}\sum_{j=-n}^{n}y_j,\quad a_k= \frac{2}{2n+1} \sum_{j=-n}^n y_j \cos k\,x_j ,\quad b_k= \frac{2}{2n+1} \sum_{j=-n}^n y_j \sin k\,x_j \quad npu \ k\in \{1,\dots,n \} .

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

Задачу тригонометрической интерполяции можно интерпретировать как построение некоторого отображения. «На вход» этого отображения подается набор значений \{y_1,\dots,y_{2n+1}\} интерполируемой функции, «на выходе» должны получить коэффициенты \{a_0,a_1,b_1,\dots,a_n,b_n\} интерполирующего тригонометрического полинома. Обсудим свойства этого отображения. Первым делом условимся, что узлы интерполяции \{x_1,\dots,x_{2n+1} \} фиксированы, а значения в них интерполируемых функций могут меняться от эксперимента к эксперименту. Запишем каждый набор значений в виде вектора, и в том же стиле запишем коэффициенты тригонометрического полинома. Для определенности будем рассматривать эти векторы как столбцы в пространстве \mathbb C^n (т.е. и значения функции и коэффициенты тригонометрического полинома могут быть комплексными), получаем некоторое соответствие

\left( \begin{array}{l} y_1\\ y_2 \\ \vdots \\ y_{2n+1} \end{array} \right) \mapsto \left( \begin{array}{l} a_0\\ a_1 \\ b_1 \\ \vdots \\ a_{n} \\ b_n \end{array} \right) \ .

Итак, вектор пространства \mathbb C^n отображается в вектор того же пространства. Это соответствие является именно функцией, причем функцией взаимно-однозначной: поскольку, по доказанному выше, каждому набору значений соответствует единственный набор коэффициентов соответствующего тригонометрического полинома. Таким образом, мы имеем дело с преобразованием пространства \mathbb C^n.

Т

Теорема. Это преобразование линейно.

Доказательство этого факта разбивается на проверку двух свойств линейности преобразования, которые мы переведем «на язык интерполяции». Умножение всех значений y_j на одну и ту же константу \alpha_{} приводит к умножению на ту же константу интерполяционного полинома

\begin{array}{c} \begin{array}{c|ccccc} x & x_1 & x_2 & \dots & x_{2\,n+1} \\ \hline y & y_1 & y_2 &\dots & y_{2\,n+1} \end{array} \\ \downarrow \\ g_n(x) \end{array} \qquad \qquad \begin{array}{c} \begin{array}{c|ccccc} x & x_1 & x_2 & \dots & x_{2\,n+1} \\ \hline y & \alpha y_1 & \alpha y_2 &\dots & \alpha y_{2\,n+1} \end{array} \\ \downarrow \\ \alpha f_n(x) \end{array}

а сложение двух интерполяционных таблиц приводит к сложению соответствующих интерполяционных полиномов:

\begin{array}{c} \begin{array}{c|ccccc} x & x_1 & x_2 & \dots & x_{2\,n+1} \\ \hline y & y_1 & y_2 &\dots & y_{2\,n+1} \end{array} \\ \downarrow \\ g_n(x) \\ \searrow \end{array} \qquad \qquad \begin{array}{c} \begin{array}{c|ccccc} x & x_1 & x_2 & \dots & x_{2\,n+1} \\ \hline y & \tilde y_1 & \tilde y_2 &\dots & \tilde y_{2\,n+1} \end{array} \\ \downarrow \\ \tilde g_n(x) \\ \swarrow \end{array}
\begin{array}{c} \begin{array}{c|ccccc} x & x_1 & x_2 & \dots & x_{2\,n+1} \\ \hline y & y_1+ \tilde y_1 & y_2 + \tilde y_2 &\dots & y_{2\,n+1} + \tilde y_{2\,n+1} \end{array} \\ \downarrow \\ g_n(x)+\tilde g_n(x) \end{array}

Оба утверждения следуют непосредственно из доказательства теоремы о существовании тригонометрического полинома:

g_n(x) = \sum_{j=1}^{2n+1} y_j W_j(x)

где тригонометрические полиномы \{ W_j(x) \}_{j=1}^{2n+1} не зависят от \{y_1,y_2,\dots,y_{2n+1}\}.

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

x_j=\frac{2(j-1)\pi}{2n+1} \quad npu \quad j\in \{1,2,\dots,2n+1\} \ ,

для которого мы получили компактные формулы в предыдущем пункте:

a_0=\frac{1}{2n+1}\sum_{j=1}^{2n+1}y_j,
a_k= \frac{2}{2n+1} \sum_{j=1}^{2n+1} y_j \cos k\,x_j ,\quad b_k= \frac{2}{2n+1} \sum_{j=1}^{2n+1} y_j \sin k\,x_j \ npu \ k\in \{1,\dots,n \} .

Обратим внимание на симметрию вхождения косинусов и синусов в две последние формулы. Вспоминая тригонометрическую форму комплексного числа, хочется объединить эти формулы в одну, с использованием мнимой единицы \mathbf i:

a_k+ \mathbf i b_k = \frac{2}{2n+1} \sum_{j=1}^{2n+1} y_j \left( \cos k\,x_j + \mathbf i \sin k\,x_j \right) =

Теперь еще одно приятное воспоминание о замечательной формуле Муавра позволит нам привести выражение к виду

=\frac{2}{2n+1} \sum_{j=1}^{2n+1} y_j \left( \cos \, x_j + \mathbf i \sin \, x_j \right)^k \ .

Но это ещё не всё, мы можем продолжить и дальше, если обратим внимание на величины x_{j}:

\cos \, x_j + \mathbf i \sin \, x_j = \cos \frac{2(j-1)\pi}{2n+1} + \mathbf i \sin \frac{2(j-1)\pi}{2n+1}

и в правой части формулы стоит объект, известный как корень из единицы степени 2n+1. Если обозначить его традиционным для теории комплексных чисел обозначением

\varepsilon_{j-1} = \cos \frac{2\pi (j-1)}{2n+1} + \mathbf i \sin \frac{2\pi (j-1)}{2n+1} \quad npu \quad j\in \{1,2,\dots,2n+1 \}

то мы получаем окончательную формулу

a_k+ \mathbf i b_k = \frac{2}{2n+1} \sum_{j=1}^{2n+1} y_j \varepsilon_{j-1}^k \quad npu \quad k\in \{1,2,\dots,n \} \ .

Теперь объединим формулы для a_{k} и b_{k} по-другому и преобразуем выражение a_k- \mathbf i b_k:

a_k- \mathbf i b_k = \frac{2}{2n+1} \sum_{j=1}^{2n+1} y_j \left( \cos k\,x_j - \mathbf i \sin k\,x_j \right) = \frac{2}{2n+1} \sum_{j=1}^{2n+1} y_j \left( \cos (-k\,x_j) + \mathbf i \sin (-k\,x_j) \right)=
= \frac{2}{2n+1} \sum_{j=1}^{2n+1} y_j \left( \cos (-\,x_j) + \mathbf i \sin (-\,x_j) \right)^k=\frac{2}{2n+1} \sum_{j=1}^{2n+1} y_j \varepsilon_{-j+1}^k \quad npu \quad k\in \{1,2,\dots,n \} \ .

Здесь комплексное число

\varepsilon_{-j+1} =\cos \frac{2(j-1)\pi}{2n+1} - \mathbf i \sin \frac{2(j-1)\pi}{2n+1}

также является корнем степени 2n+1 из единицы.

Числа (a_k+ \mathbf i b_k)/2 и (a_k- \mathbf i b_k)/2 уже встречались нам ВЫШЕ и являлись коэффициентами алгебраического полинома

G_{2n}(z)=\underbrace{\frac{a_n+\mathbf i b_n}{2}}_{u_0}+\underbrace{\frac{a_{n-1}+\mathbf i b_{n-1}}{2}}_{u_1}z+\dots+\underbrace{\frac{a_1+\mathbf i b_1}{2}}_{u_{n-1}}z^{n-1}+ \underbrace{a_0}_{u_n}z^n+\underbrace{\frac{a_1-\mathbf i b_1}{2}}_{u_{n+1}}z^{n+1}+\dots+\underbrace{\frac{a_{n-1}-\mathbf i b_{n-1}}{2}}_{u_{2n-1}}z^{2n-1} + \underbrace{\frac{a_n-\mathbf i b_n}{2}}_{u_{2n}}z^{2n} \ ,

связанного с тригонометрическим полиномом g_n(x) соотношением

g_n(x) \equiv e^{-\mathbf i n x} G_{2n} (e^{\mathbf i x}) \ .

Оказывается, что наша задача тригонометрической интерполяции удобнее всего записывается именно в терминах коэффициентов полинома G_{2n}(z). Собственно, мы почти дошли до цели: мы собираемся рассматривать отображение вектора значений интерполяционной таблицы в вектор коэффициентов полинома G_{2n}(z). Осталось сделать последнее упрощение: с не очень понятными пока последствиями мы переставим местами эти коэффициенты:

\left( \begin{array}{l} y_1\\ y_2 \\ \vdots \\ y_{2n+1} \end{array} \right) \mapsto \left( \begin{array}{l} u_n\\ u_{n+1} \\ \vdots \\ u_{2n} \\ u_0 \\ u_1 \\ \vdots \\ u_{n-1} \end{array} \right) \ .

Теперь выпишем формулы, связывающие эти столбцы. Первые n_{} компонент второго столбца

u_n=a_0=(y_1+\dots+y_{2n+1})/(2n+1)
2 u_{n+1}=a_1-\mathbf i b_1=(y_1+y_2\varepsilon_{-1}+y_3\varepsilon_{-2}+\dots+y_{2n+1}\varepsilon_{-2n})/(2n+1)
2 u_{n+2}=a_2-\mathbf i b_2=(y_1+y_2\varepsilon_{-1}^2+y_3\varepsilon_{-2}^2+\dots+y_{2n+1}\varepsilon_{-2n}^2)/(2n+1)
\dots
2 u_{2n}=a_n-\mathbf i b_n=(y_1+y_2\varepsilon_{-1}^n+y_3\varepsilon_{-2}^n+\dots+y_{2n+1}\varepsilon_{-2n}^n)/(2n+1)

Теперь рассмотрим следующую компоненту столбца коэффициентов

2u_0=a_n+\mathbf i b_n=(y_1+y_2\varepsilon_{1}^n+y_3\varepsilon_{2}^n+\dots+y_{2n+1}\varepsilon_{2n}^n)/(2n+1) \ .

Заметим, что корни степени 2n+1 из единицы удовлетворяют легко проверяемым соотношениям

\varepsilon_{1}^n=\varepsilon_{-1}^{n+1}, \varepsilon_{2}^n=\varepsilon_{-2}^{n+1},\dots \varepsilon_{2n}^n=\varepsilon_{-2n}^{n+1} \ ,

и, таким образом, последняя формула эквивалентна

2u_0=a_n+\mathbf i b_n=(y_1+y_2\varepsilon_{-1}^{n+1}+y_3\varepsilon_{-2}^{n+1}+\dots+y_{2n+1}\varepsilon_{-2n}^{n+1})/(2n+1) \ .

Пойдем дальше аналогичными преобразованиями

2u_1=a_{n+1}+\mathbf i b_{n+1}=(y_1+y_2\varepsilon_{-1}^{n+2}+y_3\varepsilon_{-2}^{n+2}+\dots+y_{2n+1}\varepsilon_{-2n}^{n+2})/(2n+1) \ ,
\dots
2u_{n-1}=a_{2n}+\mathbf i b_{2n}=(y_1+y_2\varepsilon_{-1}^{2n}+y_3\varepsilon_{-2}^{2n}+\dots+y_{2n+1}\varepsilon_{-2n}^{2n})/(2n+1) \ .

Теперь переписываем формулы в матричном виде

(2n+1)\left( \begin{array}{rl} & u_n \\ 2 & u_{n+1} \\ 2 & u_{n+2} \\ & \vdots \\ 2 & u_{2n} \\ 2 & u_0 \\ 2 & u_1 \\ & \vdots \\ 2 & u_{n-1} \end{array} \right) = \left( \begin{array}{lllll} 1 & 1 & 1 & \dots & 1 \\ 1 & \varepsilon_{-1} & \varepsilon_{-2} & \dots & \varepsilon_{-2n} \\ 1 & \varepsilon_{-1}^2 & \varepsilon_{-2}^2 & \dots & \varepsilon_{-2n}^2 \\ & \dots & & & \dots \\ 1 & \varepsilon_{-1}^n & \varepsilon_{-2}^n & \dots & \varepsilon_{-2n}^n \\ 1 & \varepsilon_{-1}^{n+1} & \varepsilon_{-2}^{n+1} & \dots & \varepsilon_{-2n}^{n+1} \\ & \dots & & & \dots \\ 1 & \varepsilon_{-1}^{2n} & \varepsilon_{-2}^{2n} & \dots & \varepsilon_{-2n}^{2n} \end{array} \right) \left( \begin{array}{l} y_1 \\ y_2 \\ \vdots \\ y_{2n+1} \end{array} \right)

и видим, что в правой части возникла матрица Вандермонда так характерная для задачи алгебраической интерполяции.

Пришло время давать самое главное определение.


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

Пусть N_{} — произвольное натуральное число. Говорят, что cистема (последовательность) чисел \{ \mathfrak X_0,\dots, \mathfrak X_{N-1} \} получается из системы (последовательности) чисел \{ \mathfrak x_0,\dots, \mathfrak x_{N-1} \} в результате дискретного преобразования Фурье (ДПФ)1), если

\left( \begin{array}{l} \mathfrak X_0 \\ \mathfrak X_1 \\ \vdots \\ \mathfrak X_{N-1} \end{array} \right) = \left( \begin{array}{lllll} 1 & 1 & 1 & \dots & 1 \\ 1 & \varepsilon_{-1} & \varepsilon_{-2} & \dots & \varepsilon_{-N+1} \\ 1 & \varepsilon_{-1}^2 & \varepsilon_{-2}^2 & \dots & \varepsilon_{-N+1}^2 \\ & \dots & & & \dots \\ 1 & \varepsilon_{-1}^{N-1} & \varepsilon_{-2}^{N-1} & \dots & \varepsilon_{-N+1}^{N-1} \\ \end{array} \right)_{N\times N} \left( \begin{array}{l} \mathfrak x_0 \\ \mathfrak x_1 \\ \vdots \\ \mathfrak x_{N-1} \end{array} \right) \quad npu \quad \varepsilon_{-j}= \cos \frac{2\pi j}{N} - \mathbf i \sin \frac{2\pi j}{N} \ .

или, в эквивалентном виде,

\mathfrak X_k=\sum_{j=0}^{N-1} \mathfrak x_j \varepsilon_{-j}^k \quad npu \quad k\in\{0,\dots, N-1\} \ .

Будем обозначать ДПФ буквой \mathcal F:

\mathcal F(\mathfrak x_0,\dots, \mathfrak x_{N-1}) = (\mathfrak X_0,\dots, \mathfrak X_{N-1}) \ ;

оно, очевидно, является линейным преобразованием пространства \mathbb C^N. Естественно считать наборы чисел (\mathfrak x_0,\dots, \mathfrak x_{N-1}) и (\mathfrak X_0,\dots, \mathfrak X_{N-1}) векторами этого пространства; но в ТЕОРИИ СИГНАЛОВ принято говорить о последовательности \{ \mathfrak x_j \}_{j=0}^{N-1} как о входной последовательности; а о последовательности \{ \mathfrak X_k \}_{k=0}^{N-1} — как о выходной последовательности. Еще раз подчеркнем смысл чисел \{\varepsilon_{-j}\}_{j=0}^{N-1}: это — система всех различных корней N_{}-й степени из 1_{} ; только в отличие от стандартного определения, они пронумерованы в обратном порядке (на единичной окружности комплексной плоскости движемся по часовой стрелке). Поскольку имеет место равенство

\varepsilon_{-j}=\varepsilon_{-1}^j = \varepsilon_{1}^{-j} \ ,

то дискретное преобразование Фурье можно переписать еще в двух эквивалентных видах:

\mathfrak X_k=\sum_{j=0}^{N-1} \mathfrak x_j \varepsilon_{-1}^{jk}= \sum_{j=0}^{N-1} \mathfrak x_j \varepsilon_{1}^{-jk} \quad npu \quad k\in\{0,\dots, N-1\} \ .

Наконец, использование экспоненциального представления корня из единицы \varepsilon_{-j} =e^{-\mathbf i (2\pi j/N)} дает еще одну часто используемую формулу:

\mathfrak X_k=\sum_{j=0}^{N-1} \mathfrak x_j e^{-\mathbf i(2\pi jk/N)} \quad npu \quad k\in\{0,\dots, N-1\} \ .
П

Найти ДПФ набора (-1, 2,1,-3).

Решение. Здесь N_{}=4, а корни четвертой степени из единицы хорошо известны: \{\varepsilon_{-j}\}_{j=0}^3=\{1,-\mathbf i,-1, \mathbf i\}.

\left( \begin{array}{rrrr} 1 & 1 & 1 & 1 \\ 1 & -\mathbf i & -1 & \mathbf i \\ 1 & -1 & 1 & -1 \\ 1 & \mathbf i & -1 & -\mathbf i \end{array} \right) \left( \begin{array}{r} -1 \\ 2 \\ 1 \\ -3 \end{array} \right) = \left( \begin{array}{c} -1 \\ -2-5\, \mathbf i \\ 1 \\ -2+5\, \mathbf i \end{array} \right) .

Ответ. ( -1,\ -2-5\, \mathbf i, \ 1, \ -2+5\, \mathbf i ).


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

\mathcal F(y_1,\dots, y_{2n+1})=(A_0,\dots, A_{2n}) \ ,

то интерполяционный тригонометрический полином строится по формуле

g_n(x) \equiv \frac{1}{2n+1} \left[ A_0 +A_1 z+A_2z^2+\dots+A_nz^n+A_{n+1}z^{-n}+A_{n+2}z^{-n+1}+\dots+A_{2n}z^{-1}\right] \quad npu \quad z=\cos \, x + \mathbf i \sin \, x \ .

Такая нумерация коэффициентов полинома позволяет развить альтернативный подход к изложению теории ДПФ. Вернемся к исходной точке: будем искать коэффициенты полинома из системы линейных уравнений. Если x_j=2\pi j /(2n+1), то соответствующее значение z_{j} будет снова корнем из 1_{} степени 2n+1, как и любая положительная степень z_{j}:

z_j= \cos \frac{2\pi j}{2n+1} +\mathbf i \sin \frac{2\pi j}{2n+1}=\varepsilon_j, z_j^n=\varepsilon_j^n \ ,

а отрицательные степени z_{j} переписываются через положительные степени \varepsilon_j:

z_j^{-n}=\varepsilon_j^{-n}=\varepsilon_j^{n+1}, z_j^{-n+1}= \varepsilon_j^{n+2},\dots, z_j^{-1}= \varepsilon_j^{2n} \ .

Следовательно, система линейных уравнений имеет вид:

\frac{1}{2n+1} \left( \begin{array}{lllll} 1 & 1 & 1 & \dots & 1 \\ 1 & \varepsilon_{1} & \varepsilon_{1}^2 & \dots & \varepsilon_{1}^{2n} \\ 1 & \varepsilon_{2} & \varepsilon_{2}^2 & \dots & \varepsilon_{2}^{2n} \\ \vdots & & & & \vdots \\ 1 & \varepsilon_{2n} & \varepsilon_{2n}^2 & \dots & \varepsilon_{2n}^{2n} \\ \end{array} \right) \left( \begin{array}{l} A_0 \\ A_1 \\ \vdots \\ A_{2n} \end{array} \right)= \left( \begin{array}{l} y_1 \\ y_2 \\ \vdots \\ y_{2n+1} \end{array} \right) \ .

При таком подходе матрица Вандермонда возникает в еще более естественной манере: фактически мы сводим задачу тригонометрической интерполяции к задаче интерполяции алгебраической. Теперь надо решить получившуюся систему. Однако ответ нам уже известен! И этот ответ — то есть ДПФ — записан снова с помощью матрицы Вандермонда, которая отличается от только что полученной только сменой знаков у индексов. Если обозначить эту последнюю через V(1,\varepsilon_{1},\varepsilon_{2},\dots, \varepsilon_{2n}), а первую — через V(1,\varepsilon_{-1},\varepsilon_{-2},\dots, \varepsilon_{-2n}) то получаем два соотношения

\frac{1}{2n+1} V(1,\varepsilon_{1},\varepsilon_{2},\dots, \varepsilon_{2n}) \left( \begin{array}{l} A_0 \\ A_1 \\ \vdots \\ A_{2n} \end{array} \right) =\left( \begin{array}{l} y_1 \\ y_2 \\ \vdots \\ y_{2n+1} \end{array} \right) \quad u \quad V(1,\varepsilon_{-1},\varepsilon_{-2},\dots, \varepsilon_{-2n}) \left( \begin{array}{l} y_1 \\ y_2 \\ \vdots \\ y_{2n+1} \end{array} \right) = \left( \begin{array}{l} A_0 \\ A_1 \\ \vdots \\ A_{2n} \end{array} \right) \ ,

откуда выводим изящную формулу:

V(1,\varepsilon_{1},\varepsilon_{2},\dots, \varepsilon_{2n})^{-1} = \frac{1}{2n+1} V(1,\varepsilon_{-1},\varepsilon_{-2},\dots, \varepsilon_{-2n}) \ .

Обращение матрицы интерполяционной системы оказывается крайне простым!


Опять-таки, отвлекаясь на время от задачи интерполяции, введем для такой удивительной матрицы специальное название. Для произвольного натурального N_{} назовем матрицу

F=\left[ \varepsilon_j^{k} \right]_{j,k=0}^{N-1}= \left( \begin{array}{lllll} 1 & 1 & 1 & \dots & 1 \\ 1 & \varepsilon_1 & \varepsilon_1^2 & \dots & \varepsilon_1^{N-1} \\ 1 & \varepsilon_2 & \varepsilon_2^2 & \dots & \varepsilon_2^{N-1} \\ 1 & \varepsilon_3 & \varepsilon_3^2 & \dots & \varepsilon_3^{N-1} \\ \vdots & & & & \vdots \\ 1 & \varepsilon_{N-1} & \varepsilon_{N-1}^{2} & \dots & \varepsilon_{N-1}^{N-1} \end{array} \right)_{N\times N} \quad npu \quad \varepsilon_j = \cos \frac{2 \pi j}{N} + {\mathbf i} \, \sin \frac{2 \pi j}{N}

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

F^{-1}=\frac{1}{N}\left[ \varepsilon_{-j}^{k} \right]_{j,k=0}^{N-1} \ .
§

Непосредственное доказательство этого свойства, а также вывод других свойств матрицы ДПФЗДЕСЬ.


§

Итак, хороший выбор узлов интерполяции и удачное переобозначение коэффициентов существенно упростили жизнь аналитическое решение задачи. В общем случае произвольных узлов для нахождения коэффициентов полинома пришлось бы решать систему линейных уравнений с матрицей, состоящей из синусов и косинусов от значений узлов (см. ☞ ЗДЕСЬ ).

Быстрое преобразование Фурье

Дискретное преобразование Фурье, рассмотренное в предыдущем пункте, представляет собой только лишь альтернативную форму записи решения задачи тригонометрической интерполяции с равноотстоящими узлами, приведенного в конце ☞ ПУНКТА — его введение не дает ни какого вычислительного преимущества.

Подсчитаем количество элементарных операций, необходимых для вычисления ДПФ в общем случае:

\mathfrak X_k=\sum_{j=0}^{N-1} \mathfrak x_j \varepsilon_{-j}^k \quad npu \quad k\in\{0,\dots, N-1\} ;

числа \{\mathfrak x_j\}_{j=0}^{N-1} предполагаются комплексными. Вычисление \mathfrak X_k «в лоб» требует N_{} умножений и N-1 сложений. Следовательно, вычисление всех чисел \{\mathfrak X_k\}_{k=0}^{N-1} требует N^2 умножений и N(N-1) сложений.

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

Прореживание по времени

П

Пример. Для случая N_{}=6 сначала перегруппируем суммы, составляющие ДПФ, отделив слагаемые, содержащие \mathfrak x_j с четными индексами от слагаемых, содержащих \mathfrak x_j с нечетными индексами; потом воспользуемся свойствами корней 6-й степени из единицы:

\varepsilon_{-3}=-1, \varepsilon_{-4}=-\varepsilon_{-1}, \varepsilon_{-5}=-\varepsilon_{-2} \ .
\begin{array}{lccccc} \mathfrak X_0=& (\mathfrak x_0+\mathfrak x_2+\mathfrak x_4) & +(\mathfrak x_1+\mathfrak x_3+\mathfrak x_5) & & \\ \mathfrak X_1=& (\mathfrak x_0+\varepsilon_{-1}^2 \mathfrak x_2+\varepsilon_{-1}^4\mathfrak x_4) & +\varepsilon_{-1}(\mathfrak x_1+\varepsilon_{-1}^2\mathfrak x_3+\varepsilon_{-1}^4\mathfrak x_5) &= (\mathfrak x_0+\varepsilon_{-2} \mathfrak x_2+\varepsilon_{-2}^2\mathfrak x_4) & + \varepsilon_{-1}(\mathfrak x_1+\varepsilon_{-2} \mathfrak x_3+\varepsilon_{-2}^2\mathfrak x_5) \\ \mathfrak X_2=& (\mathfrak x_0+\varepsilon_{-2}^2 \mathfrak x_2+\varepsilon_{-2}^4\mathfrak x_4) & +\varepsilon_{-2}(\mathfrak x_1+\varepsilon_{-2}^2\mathfrak x_3+\varepsilon_{-2}^4\mathfrak x_5) &= (\mathfrak x_0+\varepsilon_{-4} \mathfrak x_2+\varepsilon_{-4}^2\mathfrak x_4) &+ \varepsilon_{-2}(\mathfrak x_1+\varepsilon_{-4} \mathfrak x_3+\varepsilon_{-4}^2\mathfrak x_5) \\ \mathfrak X_3=& (\mathfrak x_0+\mathfrak x_2+\mathfrak x_4) & - (\mathfrak x_1+\mathfrak x_3+\mathfrak x_5) & & \\ \mathfrak X_4=& (\mathfrak x_0+\varepsilon_{-4}^2 \mathfrak x_2+\varepsilon_{-4}^4\mathfrak x_4) & +\varepsilon_{-4}(\mathfrak x_1+\varepsilon_{-4}^2\mathfrak x_3+\varepsilon_{-4}^4\mathfrak x_5) &= (\mathfrak x_0+\varepsilon_{-2} \mathfrak x_2+\varepsilon_{-2}^2\mathfrak x_4) &- \varepsilon_{-1}(\mathfrak x_1+\varepsilon_{-2} \mathfrak x_3+\varepsilon_{-2}^2\mathfrak x_5) \\ \mathfrak X_5=& (\mathfrak x_0+\varepsilon_{-5}^2 \mathfrak x_2+\varepsilon_{-5}^4\mathfrak x_4) &+ \varepsilon_{-5}(\mathfrak x_1+\varepsilon_{-5}^2\mathfrak x_3+\varepsilon_{-5}^4\mathfrak x_5) &= (\mathfrak x_0+\varepsilon_{-4} \mathfrak x_2+\varepsilon_{-4}^2\mathfrak x_4) &- \varepsilon_{-2}(\mathfrak x_1+\varepsilon_{-4} \mathfrak x_3+\varepsilon_{-4}^2\mathfrak x_5) \end{array}

В получившихся формулах наблюдается следующая закономерность: вычисления \mathfrak X_0 и \mathfrak X_3 фактически дублируют друг друга (с различием лишь в знаке одного слагаемого); то же самое можно утверждать относительно пар \mathfrak X_1 и \mathfrak X_4, а также \mathfrak X_2 и \mathfrak X_5. Иными словами, грамотная организация процесса вычисления \mathfrak X_0, \mathfrak X_1, \mathfrak X_2 позволит фактически «даром» получить значения \mathfrak X_3, \mathfrak X_4, \mathfrak X_5.

Рассмотрим теперь повнимательней каждую из скобок ( \ ). Все скобки, содержащие \mathfrak x_0,\mathfrak x_2,\mathfrak x_4, имеют одинаковую структуру: они представляют суммы этих величин, домноженных на степени

\varepsilon_{-2} = \cos \frac{4\pi}{6} - \mathbf i \, \sin \frac{4\pi}{6} = \cos \frac{2\pi}{3} - \mathbf i \, \sin \frac{2\pi}{3} \ .

Последняя величина, с одной стороны является корнем 6_{}-й степени из единицы, а, с другой стороны, является корнем 3_{}-й степени из единицы. Если ввести обозначение этих корней:

\tilde \varepsilon_{0}=1,\ \tilde \varepsilon_{-1}=\cos \frac{2\pi}{3} - \mathbf i \, \sin \frac{2\pi}{3},\ \tilde \varepsilon_{-2}=\cos \frac{4\pi}{3} - \mathbf i \, \sin \frac{4\pi}{3} \ ,

то скобки,содержащие \mathfrak x_0,\mathfrak x_2,\mathfrak x_4, переписываются в виде

\begin{array}{rl} \mathfrak Y_0 = & \mathfrak x_0 +\mathfrak x_2\, \tilde \varepsilon_{0}+\mathfrak x_4\, \tilde \varepsilon_{0}^2, \\ \mathfrak Y_1 = & \mathfrak x_0 +\mathfrak x_2\, \tilde \varepsilon_{-1}+\mathfrak x_4\, \tilde \varepsilon_{-1}^2, \\ \mathfrak Y_2 = & \mathfrak x_0 +\mathfrak x_2\, \tilde \varepsilon_{-2}+\mathfrak x_4\, \tilde \varepsilon_{-2}^2, \end{array}

т.е. они представляют ДПФ набора \mathfrak x_0,\mathfrak x_2,\mathfrak x_4. Отличие от исходной задачи заключается в уменьшении размерности: исходная размерность была шестой, получившаяся — третьей. Совершенно аналогично расправляемся со скобками, содержащими \mathfrak x_1,\mathfrak x_3,\mathfrak x_5:

\begin{array}{rl} \mathfrak Z_0 = & \mathfrak x_1 +\mathfrak x_3\, \tilde \varepsilon_{0}+\mathfrak x_5\, \tilde \varepsilon_{0}^2, \\ \mathfrak Z_1 = & \mathfrak x_1 +\mathfrak x_3\, \tilde \varepsilon_{-1}+\mathfrak x_5\, \tilde \varepsilon_{-1}^2, \\ \mathfrak Z_2 = & \mathfrak x_1 +\mathfrak x_3\, \tilde \varepsilon_{-2}+\mathfrak x_5\, \tilde \varepsilon_{-2}^2; \end{array}

они также представляют ДПФ набора \mathfrak x_1,\mathfrak x_3,\mathfrak x_5. Окончательно:

\begin{array}{ccc} \mathfrak X_0=\mathfrak Y_0+\mathfrak Z_0, & \mathfrak X_1= \mathfrak Y_1+\varepsilon_{-1}\mathfrak Z_1, & \mathfrak X_2= \mathfrak Y_2+\varepsilon_{-2}\mathfrak Z_2, \\ \mathfrak X_3=\mathfrak Y_0-\mathfrak Z_0, & \mathfrak X_4= \mathfrak Y_1-\varepsilon_{-1}\mathfrak Z_1, & \mathfrak X_5= \mathfrak Y_2-\varepsilon_{-2}\mathfrak Z_2. \end{array}

Полученные формулы сводят вычисление ДПФ от шестого порядка к третьему.

Идея очевидно распространяется на случай произвольного четного числа N_{}:

\mathfrak X_k = \left\{ \begin{array}{llll} \mathfrak Y_k & +\varepsilon_{-k} & \times \mathfrak Z_k & npu \ k\in \{0,1,\dots,N/2-1\}, \\ \mathfrak Y_{k-N/2}& -\varepsilon_{-(k-N/2)}& \times \mathfrak Z_{k-N/2} & npu \ k\in \{N/2,N/2+1,\dots,N-1\}; \end{array} \right.

при

\varepsilon_{-k}=\cos \frac{2\pi k}{N}- \mathbf i \sin \frac{2\pi k}{N}=e^{-\mathbf i (2\pi k/N)} \ .

Здесь набор \mathfrak Y_0,\mathfrak Y_1, \dots , \mathfrak Y_{N/2-1} представляет собой ДПФ от набора \mathfrak x_0,\mathfrak x_2,\dots, \mathfrak x_{N-2}, а набор \mathfrak Z_0,\mathfrak Z_1, \dots , \mathfrak Z_{N/2-1}ДПФ от набора \mathfrak x_1,\mathfrak x_3,\dots, \mathfrak x_{N-1}; число \varepsilon_{-k} называется дополнительным множителем2). Если число N/2 само является четным, то процедура повторяется. Если число N_{} представляет собой степень двойки: N=2^{\eta}, то рекурсивный алгоритм сведет вычисление ДПФ порядка N_{} к вычислению ДПФ порядка 2_{}. Этот алгоритм вычисления ДПФ называется бинарным прореживанием по времени3).

П

Пример. Для случая N_{}=8 схема вычислений будет следующей:

Подсчитаем теперь стоимость бинарного прореживания в элементарных операциях. Если считать, что ДПФ порядка N/2 уже вычислены, то для вычисления ДПФ порядка N_{} требуется 2\, N операций ( N_{} умножений величин \mathfrak Z_{\ell} на дополнительные сомножители и N_{} сложений). Аналогично, для вычисления каждого ДПФ порядка N/2 потребуется 2(N/2) операций; и нам требуется вычисление двух таких ДПФ (см. предыдущий пример). Таким образом, в случае N=2^{\eta} общее количество операций будет равно

(2\,N)+2\left(2\frac{N}{2}\right)+4\left(2\frac{N}{4}\right)+\dots \ ,

и сумма завершается через \eta шагов когда мы добираемся до ДПФ вторых порядков, каждое из которых стоит две операции4). В результате общая стоимость вычислений 2\,N \, \eta = 2\, N \log_2 N, что существенно меньше приведенной в начале пункта стоимости вычислений «в лоб». Так, при N=2^{10} имеем

N^2=1\, 048\, 576, \quad N \log_2 N = 10\, 240 \ .

Прореживание по частоте

Другой подход к уменьшению стоимости вычислений ДПФ проиллюстрируем снова на примере.

П

Пример. Для случая N_{}=6 сначала сгруппируем отдельно входящие в состав ДПФ равенства для \mathfrak X_k с четными индексами и равенства для \mathfrak X_k с нечетными индексами; потом объединим скобками слагаемые, содержащие пары \mathfrak x_0 и \mathfrak x_3, \mathfrak x_1 и \mathfrak x_4, \mathfrak x_2 и \mathfrak x_5:

\left\{\begin{array}{lrrrr} \mathfrak X_0& &=(\mathfrak x_0+\mathfrak x_3) &+(\mathfrak x_1+\mathfrak x_4) &+(\mathfrak x_2 +\mathfrak x_5) \\ \mathfrak X_2&=\mathfrak x_0+\varepsilon_{-2} \mathfrak x_1+\varepsilon_{-2}^2 \mathfrak x_2+\varepsilon_{-2}^3 \mathfrak x_3+\varepsilon_{-2}^4 \mathfrak x_4+\varepsilon_{-2}^5 \mathfrak x_5 &= (\mathfrak x_0+\mathfrak x_3)&+\varepsilon_{-2}(\mathfrak x_1+\mathfrak x_4)&+ \varepsilon_{-2}^2(\mathfrak x_2+\mathfrak x_5)\\ \mathfrak X_4&=\mathfrak x_0+\varepsilon_{-4} \mathfrak x_1+\varepsilon_{-4}^2 \mathfrak x_2+\varepsilon_{-4}^3 \mathfrak x_3+\varepsilon_{-4}^4 \mathfrak x_4+\varepsilon_{-4}^5 \mathfrak x_5 & = (\mathfrak x_0+\mathfrak x_3)&+\varepsilon_{-4}(\mathfrak x_1+\mathfrak x_4)&+ \varepsilon_{-4}^2(\mathfrak x_2+\mathfrak x_5) \end{array} \right.

Видим, что получившиеся формулы имеют снова структуру ДПФ — только, в отличие от исходного — не 6-го, а 3_{}-го порядка: величина

\varepsilon_{-2} = \cos \frac{4\pi}{6} - \mathbf i \, \sin \frac{4\pi}{6} = \cos \frac{2\pi}{3} - \mathbf i \, \sin \frac{2\pi}{3}

является корнем 3_{}-й степени из единицы, а входной последовательностью является \mathfrak x_0+\mathfrak x_3, \mathfrak x_1+\mathfrak x_4, \mathfrak x_2+\mathfrak x_5.

Для \mathfrak X_k c нечетными индексами, имеем:

\left\{ \begin{array}{lrrrr} \mathfrak X_1 &=\mathfrak x_0+\varepsilon_{-1}\mathfrak x_1+\varepsilon_{-1}^2\mathfrak x_2 +\varepsilon_{-1}^3\mathfrak x_3+\varepsilon_{-1}^4\mathfrak x_4+\varepsilon_{-1}^5\mathfrak x_5 &=(\mathfrak x_0-\mathfrak x_3)&+\varepsilon_{-1}(\mathfrak x_1-\mathfrak x_4)&+\varepsilon_{-1}^2 (\mathfrak x_2-\mathfrak x_5) \\ \mathfrak X_3&=\mathfrak x_0+\varepsilon_{-3}\mathfrak x_1+\varepsilon_{-3}^2\mathfrak x_2 +\varepsilon_{-3}^3\mathfrak x_3+\varepsilon_{-3}^4\mathfrak x_4+\varepsilon_{-3}^5\mathfrak x_5 & =(\mathfrak x_0-\mathfrak x_3)&+\varepsilon_{-3}(\mathfrak x_1-\mathfrak x_4)&+\varepsilon_{-3}^2 (\mathfrak x_2-\mathfrak x_5) \\ \mathfrak X_5&=\mathfrak x_0+\varepsilon_{-5}\mathfrak x_1+\varepsilon_{-5}^2\mathfrak x_2 +\varepsilon_{-5}^3\mathfrak x_3+\varepsilon_{-5}^4\mathfrak x_4+\varepsilon_{-5}^5\mathfrak x_5 & =(\mathfrak x_0-\mathfrak x_3)&+\varepsilon_{-5}(\mathfrak x_1-\mathfrak x_4)&+ \varepsilon_{-5}^2(\mathfrak x_2-\mathfrak x_5) \end{array} \right.

Здесь труднее увидеть ДПФ, но всё-таки можно: если за входную последовательность взять \mathfrak x_0-\mathfrak x_3,\varepsilon_{-1}(\mathfrak x_1-\mathfrak x_4),\varepsilon_{-1}^2(\mathfrak x_2-\mathfrak x_5), то снова получим ДПФ третьего порядка! Таким образом, перегруппировка формул позволяет свести вычисление ДПФ от шестого порядка к третьему.

Идея очевидно распространяется на случай произвольного четного числа N_{}: на основе входной последовательности \mathfrak x_0,\mathfrak x_1,\dots, \mathfrak x_{N-1} формируются две новые последовательности из N/2 элементов:

\begin{array}{lr} \mathfrak y_j= & \mathfrak x_j+\mathfrak x_{N/2+j} \\ \mathfrak z_j= & \varepsilon_{-j}(\mathfrak x_j-\mathfrak x_{N/2+j}) \end{array} \quad npu \quad j\in \{1,2,\dots, N/2\} ;\quad \varepsilon_{-j}=\cos \frac{2\pi j}{N}- \mathbf i \sin \frac{2\pi j}{N}=e^{-\mathbf i (2\pi j/N)} .

К каждой из них применяется ДПФ порядка N/2, из выходных последовательностей формируется \mathfrak X_0,\mathfrak X_1,\dots, \mathfrak X_{N-1}.

В отличие от предыдущего метода, основанного на измельчении входной последовательности \{\mathfrak x_j\}_{j=0}^{N-1}, в настоящем методе предлагается измельчать выходную последовательность \{\mathfrak X_k\}_{k=0}^{N-1}; он называется бинарным прореживанием по частоте5).

§

Идею метода может прояснить и аналогия с задачей алгебраической (степенной) интерполяции ☞ ЗДЕСЬ.

§

В обоснование терминологии «прореживание по времени», «прореживание по частоте» следует обратиться к исходной задаче тригонометрической интерполяции. Напомним, что для этой задачи индексы у \mathfrak x_j отвечают за порядковый номер узла интерполяции — т.е., в случае, когда основная переменная x_{} интерпретируется как время, — за отсчеты времени; индексы у \mathfrak X_k нумеруют кратные углы косинусов и синусов в тригонометрическом полиноме, т.е. частоты колебаний.


Оба рассмотренных метода относятся к классу методов вычисления ДПФ, известных под названием быстрого преобразования Фурье6).

Интерполяция с кратными узлами

Задача. Построить тригонометрический полином порядка n_{}, имеющий в узлах интерполяции x_{1},\dots,x_n заданные значения и заданные значения своей производной:

g_n(x_j)=y_j, g_n^{\prime}(x_j)=y_j^{\prime} \quad npu \quad j\in\{1,\dots,n\} \ .

Вспомним, что тригонометрический полином порядка n_{} в общем случае имеет 2n+1 коэффициент. Поставленная задача формирует условия на эти коэффициенты в виде линейных уравнений. Поскольку число 2n_{} этих уравнений меньше числа определяемых коэффициентов, то задача, скорее всего (см. ☞ ЗДЕСЬ ), будет иметь бесконечное множество решений. Так оно и оказывается.

Источники

Оппенгейм А., Шафер Р. Цифровая обработка сигналов. М.Техносфера. 2009.

Полиа Г., Сеге Г. Задачи и теоремы из анализа. М.Наука. 1978. Часть 2, отдел 6, сс. 85-90.

Турецкий А.Х. Теория интерполирования в задачах. Минск. Вышэйшая школа. 1968.

1) Discrete Fourier Transform (DFT).
2) Twiddle factor.
3) Radix-2 decimation-in-time (DIT), что дословно переводися как «бинарная децимация по времени»; не очень удачное, на мой взгляд, название, поскольку латинское слово decimatio означает казнь каждого десятого, а вовсе не второго
4) На самом деле, одну, но мы делаем грубые оценки сверху.
5) Radix-2 decimation-in-frequency (DIF)
6) Fast Fourier Transform (FFT).

2017/03/06 09:02 редактировал au