Вспомогательная страница к разделу ((:interpolation ИНТЕРПОЛЯЦИЯ)). Cущественно используются материалы разделов ((:complex_num КОМПЛЕКСНЫЕ ЧИСЛА)) и ((:polynomial ПОЛИНОМ ОДНОЙ ПЕРЕМЕННОЙ)). ==Тригонометрическая интерполяция== ~~TOC~~ Этот раздел посвящен периодическим функциям. Они рассматриваются для случая периода равного $ 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 \ . $$ Для общего случая доказано ☞ ((:complex_num:vspom2 ЗДЕСЬ)). !!Т!! **Теорема.** //Произвольный тригонометрический полином// $ 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 \ , $$ Тогда, на основании свойств операций над комплексными числами (в том числе ((:complex_num#формула_муавра формулы Муавра)) ), имеем $$ 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 $ (см. ☞((:complex_num#экспоненциальное_представление_комплексного_числа ЗДЕСЬ)) ). ♦ !!=>!! Если все коэффициенты $ 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_{} $ (т. е. ((:modular#сравнения сравнимые по модулю)) $ 2\pi_{} $). Определение **кратного корня** тригонометрического полинома вводится ((:polynomial#основная_теорема_высшей_алгебры аналогично)) случаю полинома алгебраического. Следующая теорема является аналогом ((:polynomial#основная_теорема_высшей_алгебры основной теоремы высшей алгебры)) для степенных полиномов. !!Т!! **Теорема.** //Тригонометрический полином// $ n_{} $//-го порядка, у которого старшие коэффициенты// $ a_{n} $ и $ b_n $ //удовлетворяют условиям// $ a_n+ \mathbf i b_n \ne 0, a_n- \mathbf i b_n \ne 0 $,// имеет в точности// $ 2n $ //комплексных корней с учетом их кратностей.// **Доказательство** ☞ ((:interpolation:dft:vspom1 ЗДЕСЬ)). ♦ !!=>!! Если числа $ \{\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} \quad \mbox{ при } \quad 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) \ . $$ **Доказательство** обоих следствий ☞ ((:interpolation:dft:vspom4 ЗДЕСЬ)). ♦ ===Тригонометрическая интерполяция== **Задача.** Построить тригонометрический полином порядка $ 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) $ удовлетворяет интерполяционной таблице очевидно: имеем полный аналог интерполяционного __степенного__ полинома в ((:interpolation#интерполяционый_полином_в_форме_лагранжа форме Лагранжа)). Образно говоря, каждое слагаемое $ 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 $$ $$ \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 \ . $$ В результате, мы избавляемся от половинного аргумента под синусами. Дальнейшие преобразования аналогичны. Окончательная формула ☞ ((:interpolation:dft:vspom2 ЗДЕСЬ)). Единственность этого полинома при условии ограничения на порядок следует из теоремы предыдущего пункта: если $ 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 $. Согласно теореме, эта система должна иметь решение относительно столбца коэффициентов, причем решение единственное. Из этого факта, в частности, следует, что определитель матрицы этой системы должнен быть отличен от нуля (см. ☞ ((:algebra2/linearsystems/cramert формулы Крамера)) ). Так оно и есть. !!Т!! **Теорема.** //Определитель матрицы системы равен// $$ 2^{2n^2} \prod_{0\le k < j \le 2n} \sin \frac{x_k-x_j}{2} \ . $$ Очень похоже на выражение для ((:algebra2:vander определителя Вандермонда)) (к которому, собственно, и сводится вычисление). Теперь обсудим возможные упрощения представления тригонометрического полинома из теоремы, которые получатся за счет удачного выбора узлов интерполяции. !!Т!! **Теорема.** //Если в качестве узлов интерполяции взять равноотстоящие// $$ 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 \} . $$ **Доказательство** ☞ ((:interpolation:dft:vspom3 ЗДЕСЬ)). Формулы последней теоремы остаются справедливыми и в случае "сдвига" системы узлов интерполяции на величину, кратную $ \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 $. !!Т!! **Теорема.** //Это преобразование линейно.// **Доказательство** этого факта разбивается на проверку двух свойств ((:mapping:operator#оператор линейности преобразования)), которые мы переведем "на язык интерполяции". Умножение всех значений $ 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} $$ Оба утверждения следуют непосредственно из доказательства теоремы ((#тригонометрическая_интерполяция1 о существовании тригонометрического полинома)): $$ 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}\} $. ♦ Теперь наша задача --- наиболее удобным способом формализовать это преобразование и выписать его ((:mapping:operator#матрица_оператора матрицу)). В общем случае, решение этой задачи громоздко и мы будем заниматься случаем равноотстоящих узлов $$ 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 \} . $$ Обратим внимание на симметрию вхождения косинусов и синусов в две последние формулы. Вспоминая ((:complex_num#тригонометрическая_форма_комплексного_числа тригонометрическую форму комплексного числа)), хочется объединить эти формулы в одну, с использованием мнимой единицы $ \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) = $$ Теперь еще одно приятное воспоминание о замечательной ((:complex_num#формула_муавра формуле Муавра)) позволит нам привести выражение к виду $$ =\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} $$ и в правой части формулы стоит объект, известный как ((:complex_num#корни_из_единицы корень из единицы степени)) $ 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) $$ и видим, что в правой части возникла ((:algebra2:vander матрица Вандермонда)) так характерная для задачи ((:interpolation алгебраической интерполяции)). Пришло время давать самое главное определение. ---- Дискретное преобразование Фурье Пусть $ N_{} $ --- произвольное натуральное число. Говорят, что cистема (последовательность) чисел $ \{ \mathfrak X_0,\dots, \mathfrak X_{N-1} \} $ получается из системы (последовательности) чисел $ \{ \mathfrak x_0,\dots, \mathfrak x_{N-1} \} $ в результате **дискретного преобразования Фурье** (**ДПФ**)[[Discrete Fourier Transform (DFT).]], если $$ \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}) \ ; $$ оно, очевидно, является ((:mapping:operator линейным преобразованием пространства)) $ \mathbb C^N $. Естественно считать наборы чисел $ (\mathfrak x_0,\dots, \mathfrak x_{N-1}) $ и $ (\mathfrak X_0,\dots, \mathfrak X_{N-1}) $ векторами этого пространства; но в ((:signal ТЕОРИИ СИГНАЛОВ)) принято говорить о последовательности $ \{ \mathfrak x_j \}_{j=0}^{N-1} $ как о **входной последовательности**; а о последовательности $ \{ \mathfrak X_k \}_{k=0}^{N-1} $ --- как о **выходной последовательности**. Еще раз подчеркнем смысл чисел $ \{\varepsilon_{-j}\}_{j=0}^{N-1} $: это --- система всех различных корней $ N_{} $-й степени из $ 1_{} $ ; только в отличие от ((:complex_num#корни_из_единицы стандартного определения)), они пронумерованы в обратном порядке (на единичной окружности комплексной плоскости движемся по часовой стрелке). Поскольку имеет место равенство $$ \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\} \ . $$ Наконец, использование ((:complex_num#экспоненциальное_представление_комплексного_числа экспоненциального представления)) корня из единицы $ \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 $, а корни четвертой степени из единицы ((:complex_num#корни_из_единицы хорошо известны)): $ \{\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}) \ . $$ ((:algebra2#обращение_матрицы Обращение)) матрицы интерполяционной системы оказывается крайне простым! ---- Опять-таки, отвлекаясь на время от задачи интерполяции, введем для такой удивительной матрицы специальное название. Для произвольного натурального $ 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} \ . $$ !!§!! ''Непосредственное доказательство этого свойства, а также вывод других свойств матрицы ДПФ'' ☞ ((:algebra2:fourier ЗДЕСЬ)). ---- Итак, хороший выбор узлов интерполяции и удачное переобозначение коэффициентов существенно упростили жизнь аналитическое решение задачи. В общем случае произвольных узлов для нахождения коэффициентов полинома пришлось бы решать систему линейных уравнений с матрицей, состоящей из синусов и косинусов от значений узлов (см. ☞ ((#тригонометрическая_интерполяция1 ЗДЕСЬ)) ). ===Быстрое преобразование Фурье== Дискретное преобразование Фурье, рассмотренное в предыдущем пункте, представляет собой только лишь альтернативную форму записи решения задачи тригонометрической интерполяции с равноотстоящими узлами, приведенного в конце ☞ ((#тригонометрическая_интерполяция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} $ называется **дополнительным множителем**[[Twiddle factor.]]. Если число $ N/2 $ само является четным, то процедура повторяется. Если число $ N_{} $ представляет собой степень двойки: $ N=2^{\eta} $, то рекурсивный алгоритм сведет вычисление **ДПФ** порядка $ N_{} $ к вычислению **ДПФ** порядка $ 2_{} $. Этот алгоритм вычисления **ДПФ** называется **бинарным прореживанием по времени**[[Radix-2 decimation-in-time (DIT), что дословно переводися как "бинарная децимация по времени"; не очень удачное, на мой взгляд, название, поскольку латинское слово //decimatio// означает казнь каждого //десятого//, а вовсе не //второго//...]]. !!П!! **Пример.** Для случая $ N_{}=8 $ схема вычислений будет следующей: {{ interpolation:dft:dft81.jpg |}} ♦ Подсчитаем теперь стоимость бинарного прореживания в элементарных операциях. Если считать, что **ДПФ** порядка $ 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 $ шагов когда мы добираемся до **ДПФ** вторых порядков, каждое из которых стоит две операции[[На самом деле, одну, но мы делаем грубые оценки сверху.]]. В результате общая стоимость вычислений $ 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} $; он называется **бинарным прореживанием по частоте**[[Radix-2 decimation-in-frequency (DIF)]]. !!§!! Идею метода может прояснить и аналогия с задачей алгебраической (степенной) интерполяции ☞ ((:interpolation:symm ЗДЕСЬ)). В обоснование терминологии //"прореживание по времени"//, //"прореживание по частоте"// следует обратиться к исходной задаче тригонометрической интерполяции. Напомним, что для этой задачи индексы у $ \mathfrak x_j $ отвечают за порядковый номер узла интерполяции --- т.е., в случае, когда основная переменная $ x_{} $ интерпретируется как время, --- за отсчеты времени; индексы у $ \mathfrak X_k $ нумеруют кратные углы косинусов и синусов в тригонометрическом полиноме, т.е. частоты колебаний. ---- Оба рассмотренных метода относятся к классу методов вычисления **ДПФ**, известных под названием **быстрого преобразования Фурье**[[Fast Fourier Transform (FFT).]]. ===Интерполяция с кратными узлами== **Задача.** Построить тригонометрический полином порядка $ 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_{} $ этих уравнений меньше числа определяемых коэффициентов, то задача, скорее всего (см. ☞ ((:algebra2:linearsystems#теорема_кронекера-капелли ЗДЕСЬ)) ), будет иметь бесконечное множество решений. Так оно и оказывается. ==Источники== **Оппенгейм А., Шафер Р.** //Цифровая обработка сигналов.// М.Техносфера. 2009. **Полиа Г., Сеге Г.** //Задачи и теоремы из анализа.// М.Наука. 1978. Часть 2, отдел 6, сс. 85-90. **Турецкий А.Х.** //Теория интерполирования в задачах.// Минск. Вышэйшая школа. 1968.