Преобразование Фурье в цифровой обработке сигналов. Часть 1: Дискретное преобразование Фурье (ДПФ)

В этой статье мы начинаем погружение в одну из фундаментальных тем цифровой обработки сигналов (ЦОС) — дискретное преобразование Фурье (ДПФ). Именно ДПФ служит математической основой для понимания более сложных методов спектрального анализа и является отправной точкой для изучения всех остальных видов преобразования Фурье в ЦОС. Будет рассмотрено ДПФ действительных сигналов.

Материал построен так, чтобы объединить теорию, наглядные графики и практический код на Python.

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

Краткое содержание статьи. Часть 1

В ней поэтапно разбираются фундаментальные понятия, необходимые для понимания ДПФ:

Аналоговые и дискретные сигналы — краткое введение, объясняющее разницу между двумя типами сигналов. В статье используются только дискретные сигналы, поэтому эта тема кратко рассмотрена.
Корреляция и ортогональность функций — ключевые понятия, лежащие в основе ДПФ. Именно эти свойства позволяют разложить сложный сигнал на отдельные частотные составляющие.
Зависимость корреляции от сдвига фазы — на примерах показано, как изменение фазы влияет на результат скалярного произведения и почему это важно для анализа сигналов.
Формула ДПФ — финальный этап, где все изученные понятия сводятся к классической формуле дискретного преобразования Фурье.

Аналоговые и дискретные сигналы

Аналоговый (непрерывный) сигнал — это способ представления информации, при котором физическая величина (напряжение, звуковое давление, температура) изменяется непрерывно во времени. В любой момент времени сигнал имеет какое-то значение.

Простой пример: температура воздуха в течение дня. Она меняется постепенно: +15,0°, +15,1°, +15,2°… Вы можете измерить её в любой момент — значение всегда будет определено. Это естественная, «природная» форма сигнала.

Ключевое свойство: аналоговый сигнал непрерывен во времени, поэтому его также называют непрерывным сигналом.

Дискретный сигнал — это представление того же самого, но только в отдельные, заранее выбранные моменты времени.

Простой пример: вы измеряете температуру не постоянно, а только раз в час. У вас будет ряд чисел: +15,2° в 12:00, +15,5° в 13:00 и т. д. Что происходило между этими измерениями — остаётся неизвестным.

Как получается дискретный сигнал: аналоговый сигнал (например, непрерывная звуковая волна) пропускают через специальную схему — аналого-цифровой преобразователь (АЦП). АЦП с определённой частотой измеряет мгновенное значение сигнала и записывает его.

Ключевое свойство: дискретный сигнал дискретен во времени. Он известен только в виде последовательности значений: Между отсчётами о сигнале нет никакой информации.

Иллюстрация аналогового и дискретного сигналов на рис. 1.

Верхняя часть рисунка — классический аналоговый сигнал (плавная синусоида). Он непрерывен.
Нижняя часть рисунка — результат дискретизации той же синусоиды. Красные точки — это отсчёты (моменты, в которые сигнал был измерен). Пунктирная линия напоминает об исходном аналоговом сигнале, но в реальном дискретном сигнале есть только точки.

Частота дискретизации — это частота, с которой производятся измерения сигнала. Каждая точка — один отсчёт.

Чем выше частота дискретизации, тем больше отсчётов приходится на секунду времени и тем точнее дискретный сигнал описывает исходный аналоговый. Частота дискретизации сигнала во всех примерах статьи — 100 Гц. То есть 100 отсчётов за 1 с.

Примечание: здесь освещён только самый минимум, необходимый для дальнейшего понимания статьи. Глубокое погружение в теорию дискретизации выходит за её рамки.

Корреляция и ортогональность функций в цифровой обработке сигналов

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

Раньше для меня было настоящей магией то, как в эфире, переполненном самыми разными сигналами, можно чётко выделить отдельные сигналы и понять, какой у них состав. Но как только я разобрался с теорией, поэкспериментировал с программами и построил графики, эта магия исчезла, уступив место чёткому пониманию.

Ортогональность функций в цифровой обработке сигналов: дискретный подход

В контексте ЦОС понятие «ортогональность» является частным случаем корреляции. Два сигнала считаются ортогональными когда их взаимная корреляция равна . Проще говоря, это означает, что сигналы «не похожи» друг на друга.

Прежде чем мы углубимся в практические примеры и формулы, напомню базовый математический инструмент, без которого дальше будет сложно обойтись: что такое знак суммы $\Sigma$ («сигма») и как он используется в наших формулах.

Знак суммы

Общее понятие

Знак суммы (обозначается греческой буквой Σ — «сигма») — это математический символ, который используется для краткой записи суммы большого количества слагаемых, подчиняющихся определённому правилу.

Формула суммы выглядит так:

$\sum_{k=n}^{m} a_k\tag{1.1}$

В развёрнутом виде это читается как «сумма (a_k) для всех (k) от (n) до (m)».

Расшифровка обозначений:

Σ — знак суммы (оператор суммирования);
(a_k) — общий член последовательности, который мы суммируем;
(k) — индекс (счётчик), который меняет своё значение;
(n) — начальное значение индекса;
(m) — конечное значение индекса.

Правило вычисления

Чтобы вычислить сумму, нужно последовательно подставить в выражение (a_k) все значения (k) от (n) до (m) и сложить получившиеся результаты.

Это эквивалентно длинной записи:

$\sum_{k=n}^{m} a_k = a_n + a_{n+1} + a_{n+2} + \dots + a_{m-1} + a_m \quad (1.2)$

Наглядный пример с числами

Условие: Дана последовательность чисел: (a_1 = 5), (a_2 = 7), (a_3 = 10), (a_4 = 13). Нужно найти сумму всех этих чисел.

Решение через знак суммы:

$S = \sum_{k=1}^{4} a_k$

Развёрнутая запись:

Подставляем значения индекса (k) от 1 до 4 в выражение (a_k):

При (k = 1): (a_1 = 5)
При (k = 2): (a_2 = 7)
При (k = 3): (a_3 = 10)
При (k = 4): (a_4 = 13)

Складываем их:

Итог:

$\sum_{k=1}^{4} a_k = 35$

Дискретный подход к корреляции и ортогональности

Дискретное скалярное произведение

В ЦОС мы работаем не с непрерывными функциями, а с последовательностями отсчётов. Поэтому все операции, включая проверку ортогональности, выполняются в дискретной форме.

Формула скалярного произведения

Пусть даны две последовательности отсчётов длины :

$x[n], \quad y[n], \quad n = 0, 1, 2, \ldots, N-1$

Дискретное скалярное произведение определяется как:

$\left\langle x, y \right\rangle = \sum_{n=0}^{N-1} x[n] \cdot y[n] \quad (1.3)$

Интерпретация формулы на примере (рис. 2)

Рассмотрим вычисление скалярного произведения на примере двух косинусоидальных сигналов (верхний график на рис. 2):

Поэлементное умножение: для каждого отсчёта перемножаем значения и . Результаты показаны на графике «Поэлементное произведение» (рис. 2).
Суммирование: все полученные произведения последовательно суммируются от до . Процесс накопления суммы показан на графике «Накопленная сумма произведений» (рис. 2).
Результат: итоговое значение скалярного произведения — на рис. 2 это значение при .

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

Критерий ортогональности

Две дискретные последовательности и называются ортогональными, если их скалярное произведение строго равно нулю:

$\langle x, y \rangle = \sum_{n=0}^{N-1} x[n] \cdot y[n] = 0$

Однако на практике проверяется условие:

$|\langle x, y \rangle| < \varepsilon$

где $\varepsilon$ — малый порог (например, $10^{-10}$ ).

Зачем нужен порог

Идеальный ноль ( $\langle x, y \rangle = 0$ ) на практике часто недостижим по трём основным причинам:

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

Поэтому на практике сигналы считаются ортогональными, если их скалярное произведение практически равно нулю, то есть меньше выбранного порога $\varepsilon$ .

После введения математического определения дискретного скалярного произведения и критерия ортогональности целесообразно рассмотреть практическую иллюстрацию данных понятий. Графическая интерпретация позволяет наглядно продемонстрировать механизм обнуления суммы произведений отсчётов для ортогональных сигналов и сохранения значимой величины скалярного произведения для коррелированных сигналов.

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

Генерация тестовых сигналов

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

$x[n] = A \cdot \cos\left(2\pi \cdot \frac{k \cdot n}{N} + \varphi\right) \quad (1.4)$

$x[n] = A \cdot \sin\left(2\pi \cdot \frac{k \cdot n}{N} + \varphi\right) \quad (1.5)$

Параметры сигнала:

Символ	Значение
	Амплитуда сигнала
	Частотный индекс (число периодов на интервале )
	Длина последовательности (количество отсчётов)
$\varphi$	Начальная фаза
	Индекс отсчёта, $n = 0, \ldots, N-1$

Варьируя параметры и $\varphi$ , демонстрируют случаи ортогональных и коррелированных сигналов, а также влияние фазового сдвига на величину скалярного произведения.

В примерах статьи используются следующие параметры:

частота дискретизации: Гц (100 отсчётов в секунду);
длина последовательности: отсчётов.

Длительность наблюдаемого сигнала составляет:

$T = \frac{N}{F_s} = \frac{100}{100} = 1 \text{ с}$

В формулах (1.4) и (1.5) индекс показывает, сколько полных периодов синусоиды или косинусоиды укладывается на интервале . Соответственно, частота генерируемого сигнала определяется как:

$f = \frac{k}{T} = \frac{k \cdot F_s}{N}$

Для примеров из статьи () формула упрощается:

$f = k \cdot \frac{100}{100} = k \quad \text{(Гц)}$

То есть означает 1 период в секунду (1 Гц), — 5 периодов в секунду (5 Гц) и т. д.

Важный вывод: если перед нами лежит график, например, из 100 отсчётов, но неизвестна частота дискретизации , то определить реальную частоту сигнала невозможно. Те же 100 отсчётов могли быть получены за 1 секунду (тогда Гц), а могли и за 1 час — тогда частота сигнала в 3600 раз меньше. Без знания мы видим лишь «форму» сигнала, но не его истинную частоту.

Далее приведены примеры применения формул (1.4) и (1.5).

Визуализация ортогональных функций

Пример 1: ортогональные гармонические сигналы с различными частотными индексами

Для иллюстрации дискретного скалярного произведения и критерия ортогональности рассматривается пара гармонических сигналов, сгенерированных согласно формуле 1.4.

Параметры сигналов:

Сигнал	Амплитуда	Частотный индекс	Фаза $\varphi$	Формула
(исследуемый)				$x[n] = 2\cos\left(2\pi \cdot \frac{1 \cdot n}{N}\right)$
(базисная функция)				$y[n] = 1\cos\left(2\pi \cdot \frac{2 \cdot n}{N}\right)$

Длина последовательности отсчётов, дискретное время $n = 0, 1, \dots, N-1$ . Оба сигнала содержат целое число периодов на интервале наблюдения: и . При целых различных частотных индексах $k_1 \neq k_2$ гармонические сигналы ортогональны на интервале длины , что гарантирует нулевое значение скалярного произведения.

Структура визуализации (рис. 2)

Для анализа скалярного произведения приводится набор из трёх графиков:

Исходные сигналы и — дискретные отсчёты, отображаемые в виде дискретного графика (stem-plot, синим и красным цветами).
Произведение сигналов $p[n] = x[n] \cdot y[n]$ — последовательность значений, подлежащих суммированию.
Скалярное произведение $S[m] = \sum_{n=0}^{m} x[n] \cdot y[n]$ — кривая, отражающая динамику формирования скалярного произведения; итоговое значение $S[N-1] = \langle x, y \rangle$ .

Первый график: исходные сигналы

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

Второй график: поэлементное произведение

Для каждого момента времени мы перемножаем значения сигналов: $p[n] = x[n] \cdot y[n]$ . Результат показан фиолетовыми столбиками. Когда оба сигнала имеют одинаковый знак (оба положительные или оба отрицательные), произведение положительно. Когда знаки разные — произведение отрицательно.

Можно заметить, как положительные и отрицательные столбики чередуются. График симметричен, и визуально в сумме положительные и отрицательные площади одинаковы. Это первый намёк на то, что при суммировании всех произведений положительные и отрицательные вклады могут скомпенсировать друг друга.

Третий график: скалярное произведение

Зелёная линия показывает, как меняется сумма произведений по мере продвижения по времени:

$S = \sum_{n=0}^{N-1} x[n] \cdot y[n]$

Это «бегущая» сумма: в каждый новый момент мы добавляем очередное произведение к предыдущей сумме. Кривая начинается с нуля (при сумма равна первому произведению), затем колеблется, то поднимаясь выше нуля, то опускаясь ниже.

Оранжевая пунктирная горизонтальная линия отмечает итоговое значение скалярного произведения — сумму всех произведений, то есть .

Ключевой момент: зелёная кривая после начальных колебаний к концу интервала возвращается почти к нулю. Это значит, что положительные и отрицательные вклады практически полностью уравновесили друг друга. Оранжевая линия оказывается очень близко к нулю (в нашем примере значение скалярного произведения имеет порядок $10^{-15}$ , что является машинным нулём).

Что даёт нам внимательное сравнение графиков

Если присмотреться к тому, как устроены все три графика, можно увидеть прямую связь между ними. Первый график показывает исходные сигналы — их форму, амплитуду, частоту и фазу. Именно эти характеристики определяют, каким будет второй график — поэлементное произведение. В каждом отсчёте мы просто перемножаем значения синего и красного сигналов: если оба большие и положительные — произведение даёт положительный столбик; если один положительный, а другой отрицательный — столбик уходит в отрицательную область. Таким образом, первый график напрямую диктует, где фиолетовые столбики окажутся выше нуля, а где — ниже.

В свою очередь, второй график полностью определяет поведение третьего. Накопленная сумма (зелёная кривая) в каждый следующий момент времени добавляет к себе очередное значение произведения. Когда фиолетовые столбики положительные — зелёная кривая возрастает; когда отрицательные — она оубывает. Если положительные и отрицательные столбики примерно равны по площади, кривая колеблется около нуля и в конце возвращается к нему.

Обратите внимание на симметрию, которая здесь проявляется: положительные и отрицательные фиолетовые столбики на втором графике чередуются, и их площади визуально почти одинаковы. Именно эта симметрия — результат того, что исходные сигналы имеют целое число периодов и разные частоты — приводит к тому, что зелёная кривая после серии возрастаний и убываний возвращается к нулю.

Сопоставляя все три графика, можно в деталях проследить, как форма исходных сигналов через промежуточные произведения формирует итоговую накопленную сумму (скалярное произведение). Это не просто упражнение — именно такое сравнение даёт ключ к интуитивному пониманию того, как работает преобразование Фурье.

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

Вывод

Мы наглядно убедились, что скалярное произведение двух сигналов с частотными индексами 1 и 2 (косинусов) действительно близко к нулю. Математически это записывается как:

$\langle x, y \rangle \approx \sum_{n=0}^{N-1} x[n] \cdot y[n] \approx 0 \quad (1.6)$

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

В терминах цифровой обработки сигналов это означает, что исследуемый сигнал (сигнал 1) и базисный (эталонный) сигнал (сигнал 2) ортогональны. На практике состав исследуемого сигнала нам заранее неизвестен. Чтобы определить, присутствует ли в нём некоторая частота (например, 2 Гц), мы сравниваем его с эталонным сигналом (базисной функцией) этой частоты — вычисляем их скалярное произведение.

Если результат близок к нулю (сигналы ортогональны), это означает, что исследуемый сигнал не содержит составляющей с частотой эталонного сигнала. В нашем примере исследуемый сигнал (частота 1) не имеет ничего общего с частотой 2 — она в нём отсутствует.

Именно так ортогональность позволяет «просвечивать» сигнал и определять его частотный состав: каждый эталонный сигнал (базисная функция) выделяет только «свою» частоту, не реагируя на другие. Если бы в исследуемом сигнале присутствовала частота 2, скалярное произведение дало бы заметную ненулевую величину, указывающую на наличие этой компоненты.

Пример 2

Ортогональность синуса и косинуса с одинаковой частотой

Рассматривается случай ортогональности гармонических сигналов с идентичными частотными индексами, но с фазовым сдвигом $\pi/2$ (90°), то есть между синусом и косинусом.

Параметры сигналов

Сигнал	Тип	Амплитуда	Частотный индекс	Фаза $\varphi$	Формула
(исследуемый)	$\sin$				$x[n] = 1.0 \cdot \sin\left(2\pi \cdot \frac{1 \cdot n}{N}\right)$
(базисная функция)	$\cos$				$y[n] = 1.0 \cdot \cos\left(2\pi \cdot \frac{1 \cdot n}{N}\right)$

где отсчётов, $n = 0, 1, \ldots, N-1$ .

Вывод

Важный момент: в отличие от предыдущего примера, где ортогональность достигалась за счёт разных частот, здесь оба сигнала имеют одинаковую частоту ( Гц). Ортогональность обеспечивается исключительно фазовым сдвигом на 90° между синусом и косинусом.

Фазовый сдвиг и ортогональность: ортогональность между синусом и косинусом наблюдается именно при фазовом сдвиге 90° ( $\pi/2$ рад) и 270° ( $3\pi/2$ рад). При других значениях фазового сдвига (например, 45° или 200°) сигналы не будут ортогональны — их скалярное произведение примет ненулевое значение.

Ещё одно важное замечание: в данном примере проверка на ортогональность производится на интервале, содержащем целое число периодов. Если сигналы имеют нецелое количество периодов на исследуемом отрезке, то их скалярное произведение может оказаться отличным от нуля, и условие ортогональности не выполнится.

Пример 3

Неортогональность синусоид с нецелым числом периодов

Рассматривается случай проверки ортогональности двух гармонических сигналов одинакового типа ( $\sin$ ) с разными частотными индексами ( $k_1 \neq k_2$ ) и нулевыми начальными фазами.

Параметры сигналов

Сигнал	Тип	Амплитуда	Частотный индекс	Фаза $\varphi$	Формула
(базисная функция)	$\sin$				$x[n] = 1 \cdot \sin\left(2\pi \cdot \frac{1 \cdot n}{N}\right)$
(исследуемый)	$\sin$				$y[n] = 1 \cdot \sin\left(2\pi \cdot \frac{1.25 \cdot n}{N}\right)$

где отсчётов, $n = 0, 1, \ldots, N-1$ .

Вывод

Почему сигналы не ортогональны при разных частотах?

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

Первый график (исходные сигналы): первый сигнал совершает ровно 1 период, второй — 1,25 периода. Эта «лишняя» четверть периода нарушает симметрию.
Второй график (произведение сигналов): положительные и отрицательные значения не компенсируют друг друга полностью, поэтому суммарная площадь не уравновешена.
Третий график (скалярное произведение): кривая не возвращается к нулю к концу интервала. Итоговое значение — это и есть $\langle x, y \rangle \neq 0$ , что является прямым доказательством отсутствия ортогональности.

Эффект «утечки» простыми словами: когда измерение обрывается на неполном периоде, «хвост» сигнала не находит компенсирующей пары. Энергия как бы «утекает» за границы окна наблюдения, искажая результат. Аналогичный эффект наблюдается при любых нецелых частотных индексах: , , и т. д.

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

Условия ортогональности

Целое число периодов на интервале . Сигнал должен полностью «помещаться» в окно наблюдения — только тогда положительные и отрицательные площади произведения полностью компенсируют друг друга, давая в сумме ноль.
Соотношение частот или фаз. Либо частотные индексы являются целыми числами ( $k_1, k_2 \in \mathbb{Z}$ ), либо, при одинаковой частоте, фазовый сдвиг составляет $90^\circ$ или $270^\circ$ (синус и косинус).
Амплитуда не влияет на факт ортогональности. Значения и могут быть любыми — они лишь масштабируют итоговый результат, но не меняют факт равенства скалярного произведения нулю при выполнении первых двух условий.

Детальный анализ ситуаций с нецелочисленными частотными индексами будет представлен в разделе «Свойства преобразования Фурье» одной из последующих статей, если данный цикл публикаций получит развитие. Забегая вперёд: именно описанная выше причина лежит в основе явления, известного в цифровой обработке сигналов как утечка спектра.

Пример 4

Корреляция сигналов с одинаковой частотой и противофазой (сдвиг 180°)

Рассматривается случай максимальной (отрицательной) корреляции: два сигнала имеют идентичные частоты, но противоположные фазы. В отличие от примера 2, где сдвиг 90° обеспечивал ортогональность, здесь сдвиг 180° приводит к полному отсутствию компенсации — сигналы «усиливают» друг друга в произведении, сохраняя один знак.

Параметры сигналов

Сигнал	Тип	Амплитуда	Частотный индекс	Фаза $\varphi$	Формула
(исследуемый)	$\sin$				$x[n] = 1 \cdot \sin\left(2\pi \cdot \frac{1 \cdot n}{N}\right)$
(базисный)	$\sin$			$\pi$	$y[n] = 1 \cdot \sin\left(2\pi \cdot \frac{1 \cdot n}{N} + \pi\right) = -x[n]$

где отсчётов, $n = 0, 1, \ldots, N-1$ .

Вывод

Почему нет ортогональности?
При фазовом сдвиге 180° ( $\pi$ рад) сигналы становятся зеркальными: . Их поэлементное произведение $x[n] \cdot y[n] = -\sin^2(\dots)$ всегда $\le 0$ . Все значения лежат ниже нуля (или равны ему), поэтому положительные и отрицательные вклады не компенсируются, а складываются, сохраняя один знак.

Что показывают три графика:

График 1 (исходные сигналы): сигналы симметричны относительно оси времени — пики одного совпадают с впадинами другого.
График 2 (произведение сигналов): все столбцы направлены вниз, площадь под кривой не уравновешена — компенсации нет.
График 3 (скалярное произведение): кривая монотонно уходит от нуля, финальное значение $\langle x, y \rangle \approx -50$ — максимально возможное по модулю для данных амплитуд, что указывает на сильную отрицательную корреляцию.

Прикладная интерпретация (задача обнаружения):
Представим, что — неизвестный входной сигнал. Мы последовательно сравниваем его с набором базисных синусоид с частотными индексами $k = 0, 1, 2, 3, \ldots$ :

Для всех $k \neq 1$ скалярное произведение близко к нулю (сигналы ортогональны) $\rightarrow$ этих частот в сигнале нет.
Только при получаем значимый отклик $\langle x, y \rangle \neq 0$ $\rightarrow$ частота обнаружена.

Пример 5

Корреляция сигналов с одинаковой частотой и фазой

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

Параметры сигналов

Сигнал	Тип	Амплитуда	Частотный индекс	Фаза $\varphi$	Формула
(исследуемый)	$\sin$				$x[n] = 1.5 \cdot \sin\left(2\pi \cdot \frac{1 \cdot n}{N}\right)$
(базисный)	$\sin$				$y[n] = 1.0 \cdot \sin\left(2\pi \cdot \frac{1 \cdot n}{N}\right)$

где отсчётов, $n = 0, 1, \ldots, N-1$ .

Вывод

Почему корреляция максимальна?
При одинаковых частотах и нулевом фазовом сдвиге сигналы полностью синхронизированы: их пики совпадают. Поэлементное произведение $x[n] \cdot y[n] = 1.5 \cdot \sin^2(\dots)$ всегда $\ge 0$ — все значения лежат выше нуля (или равны ему), поэтому положительные вклады не компенсируются, а накапливаются.

Что показывают три графика:

График 1 (исходные сигналы): обе синусоиды синхронны, различаясь только амплитудой.
График 2 (произведение сигналов): все столбцы направлены вверх, площадь под кривой выше нуля — компенсации нет.
График 3 (скалярное произведение): кривая монотонно растёт, финальное значение $\langle x, y \rangle \approx 75$ — максимально возможное для данных амплитуд, что указывает на выраженную корреляцию.

Зависимость корреляции от фазы сигнала

В предыдущих примерах оба сигнала имели нулевую начальную фазу ( $\varphi = 0^\circ$ ). Но что изменится, если добавить фазовый сдвиг? Например, косинус, сдвинутый на $45^\circ$ , оказывается «посередине» между чистым синусом и чистым косинусом. В этом случае корреляции сигнала с базисными функциями $\sin$ и $\cos$ будут одинаковыми.

Почему это важно на практике? Фаза реального сигнала заранее неизвестна. Если сравнивать входной сигнал только с одним базисом (например, только с $\sin$ ), возможны следующие ситуации. Когда фаза сигнала близка к $0^\circ$ относительно базисного синуса, мы получаем максимальную корреляцию, и сигнал считается «найденным» как в примере 5. Но если фаза сигнала составляет примерно $90^\circ$ относительно базисного синуса, корреляция оказывается близка к нулю, и сигнал «пропускается» как в примере 2. При промежуточной фазе корреляция занижена, из-за чего амплитуда оценивается неверно.

По этой причине необходимо вычислять значения корреляции с обеими базисными функциями — $\sin$ и $\cos$ , что в совокупности даст полную информацию о сигнале независимо от его фазового сдвига. Любой гармонический сигнал со сдвигом фазы можно представить как линейную комбинацию синуса и косинуса той же частоты:

$A \cdot \cos(\omega t + \varphi) = B \cdot \cos(\omega t) + C \cdot \sin(\omega t) \quad (1.7)$

Следующий пример более детально рассмотрит случай косинуса, сдвинутого на $45^\circ$ .

Пример 6

В этом примере мы переходим к двойной корреляции: одновременно вычисляем скалярное произведение исследуемого сигнала с обеими ортогональными базисными функциями — $\cos$ и $\sin$ одной и той же частоты (рис. 7). Это позволяет восстановить полную информацию об амплитуде и фазе сигнала независимо от его начального сдвига.

Параметры сигналов

Сигнал	Тип	Фаза $\varphi$	Формула
(исследуемый)	$\cos$	$45^\circ$	$x[n] = 1.0 \cdot \cos\left(2\pi \cdot \frac{1 \cdot n}{N} + 45^\circ\right)$
$y_{\cos}[n]$ (базис $\cos$ )	$\cos$	$0^\circ$	$y_{\cos}[n] = 1.0 \cdot \cos\left(2\pi \cdot \frac{1 \cdot n}{N}\right)$
$y_{\sin}[n]$ (базис $\sin$ )	$\sin$	$0^\circ$	$y_{\sin}[n] = 1.0 \cdot \sin\left(2\pi \cdot \frac{1 \cdot n}{N}\right)$

где отсчётов, $n = 0, 1, \ldots, N-1$ .

От скалярного произведения к комплексной плоскости

В предыдущих разделах мы рассмотрели, как с помощью корреляции (скалярного произведения) можно определить степень сходства сигнала с базисными функциями — синусом и косинусом. Мы убедились, что синус и косинус одной частоты ортогональны, поскольку их скалярное произведение равно нулю (см. пример 2). Это свойство открывает путь к очень удобному и наглядному представлению информации о сигнале — с помощью комплексной плоскости.

Пусть у нас есть исследуемый сигнал и две базисные функции одной частоты: $\cos\left(2\pi k \frac{n}{N}\right)$ и $\sin\left(2\pi k \frac{n}{N}\right)$ . Вычислим два числа:

$\operatorname{Re} X[k]$ — корреляция (скалярное произведение) сигнала с косинусом:

$\operatorname{Re} X[k] = \sum_{n=0}^{N-1} x[n] \cdot \cos\left(2\pi k \frac{n}{N}\right) \quad (1.8)$

$\operatorname{Im} X[k]$ — корреляция (скалярное произведение) сигнала с синусом (со знаком «минус» в соответствии со стандартным инженерным соглашением для ДПФ):

$\operatorname{Im} X[k] = -\sum_{n=0}^{N-1} x[n] \cdot \sin\left(2\pi k \frac{n}{N}\right) \quad (1.9)$

Знак минус в формуле (1.9) соответствует стандартному инженерному соглашению при определении дискретного преобразования Фурье. На рис. 7 графики 2.2 и 3.2 сформированы с учётом минуса. Без учёта знака графики были бы зеркально симметричны относительно горизонтальной оси. Иными словами, на рисунке 2.2 значения в пиках составили бы не около 0,8, а −0,8, а на рисунке 3.2 значение скалярного произведения равнялось бы −35,36. Смысл знака минус будет раскрыт ниже.

Расшифровка обозначений

— значения анализируемого сигнала, $n = 0, 1, \dots, N-1$ ;
— общее количество отсчётов (длина выборки);
— номер гармоники (частота в периодах на интервал );
— индекс текущего отсчёта;
$\operatorname{Re} X[k]$ — проекция на косинус (действительная часть);
$\operatorname{Im} X[k]$ — проекция на синус (мнимая часть).

Ось, на которой откладывается $\operatorname{Re} X[k]$ , ассоциируется с косинусной базисной функцией. Ось, на которой откладывается $\operatorname{Im} X[k]$ , ассоциируется с синусной базисной функцией.

Поскольку сами базисные функции ( $\cos$ и $\sin$ одной частоты) ортогональны, соответствующие им оси также ортогональны. Это даёт нам возможность представить результат корреляции в виде точки на комплексной плоскости.

На рис. 7 (график 4.1) показано представление корреляции сигнала с косинусом и синусом как комплексного числа: горизонтальная ось — действительная часть ( $\operatorname{Re}$ ), вертикальная — мнимая часть ( $\operatorname{Im}$ ). Оси перпендикулярны, как и базисные функции. (Для справки по комплексным числам: Комплексные числа для чайников)

Такое представление является фундаментом комплексного преобразования Фурье и позволяет компактно и наглядно описывать гармонические составляющие сигнала. Для большей наглядности приведён отдельный рис. 8.

Что дают и ?

Зная две проекции, мы можем восстановить полную информацию о гармонике сигнала на данной частоте:

Амплитуда гармоники:

$\operatorname{Mag}[k] = \sqrt{(\operatorname{Re} X[k])^2 + (\operatorname{Im} X[k])^2} \quad (1.10)$

Это длина радиус-вектора от начала координат до точки $(\operatorname{Re} X[k], \operatorname{Im} X[k])$ .

Начальная фаза гармоники (угол относительно оси косинусов):

$\varphi[k] = \arctan\left(\frac{\operatorname{Im} X[k]}{\operatorname{Re} X[k]}\right) \quad (1.11)$

Угол отсчитывается от положительного направления оси $\operatorname{Re}$ против часовой стрелки.

Таким образом, пара чисел $(\operatorname{Re} X[k], \operatorname{Im} X[k])$ полностью описывает амплитуду и фазу синусоидальной компоненты сигнала.

Ось $\operatorname{Re}$ соответствует корреляции с косинусом, ось $\operatorname{Im}$ — с синусом. Каждая точка (или вектор) однозначно определяет амплитуду и фазу гармонического сигнала.

Действительная часть $\operatorname{Re} X[k]$ равна скалярному произведению сигнала с косинусом, а мнимая часть $\operatorname{Im} X[k]$ — со знаком «минус» — с синусом (знак зависит от принятого соглашения). Таким образом, мы перешли к мощному аппарату комплексных чисел, который лежит в основе современного спектрального анализа.

Механическая аналогия ортогональности

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

Представьте себе тяжёлый деревянный брусок, лежащий на идеально гладкой плоскости. Если вы начнёте толкать его строго вертикально вниз или поднимать строго вверх, он будет двигаться только по вертикали. Горизонтальные силы на него не действуют, поэтому он не поедет влево или вправо. И наоборот: если вы толкнёте брусок строго горизонтально, он не взлетит в воздух. Две оси движения — вертикальная и горизонтальная — независимы друг от друга. Сила, приложенная вдоль одной оси, никак не влияет на движение по другой. В физике такие оси называют ортогональными.

Как это связано с обработкой сигналов?

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

Здесь происходит то же самое, что и с бруском. Если мы сравниваем наш сигнал с «вертикальным» базисом (например, с косинусом), то это скалярное произведение покажет только ту «долю» косинуса, которая есть в сигнале. И это никак не повлияет на результат сравнения с «горизонтальным» базисом (синусом).

Сигнал, который является «чистым» косинусом, будет иметь большую корреляцию с базисным косинусом и нулевую корреляцию с базисным синусом. Они друг другу не мешают. И наоборот, сигнал-синус отзовётся только на свой базис.

Так же, как подъём бруска не вынуждает перемещаться в сторону, наличие в сигнале косинусоидальной составляющей не создаёт «ложной» синусоидальной компоненты. Это и есть ортогональность функций. Косинус и синус одной и той же частоты — это две независимые «оси», две независимые «силы», которые можно измерять по отдельности.

Частотный и фазовый спектры

В примере 6 также представлены два дополнительных графика (см. рис. 7, графики 5.1 и 5.2): частотный спектр и фазовый спектр.

На частотном спектре отображаются амплитуды $\operatorname{Mag}[k]$ гармонических составляющих сигнала. В рассматриваемом примере исследуемый сигнал представляет собой косинус частотой 1 Гц и амплитудой 1. По горизонтальной оси графика можно определить, какой частоте соответствует каждый пик: в данном случае наблюдается единственный выраженный пик на отметке 1 Гц, что наглядно подтверждает наличие в сигнале лишь одной гармонической компоненты.

Прежде чем перейти к деталям, важно осветить два понятия ЦОС: временная область и частотная область.

Временная область — это представление сигнала в зависимости от времени. На графиках частей 1.1 и 1.2 по горизонтальной оси откладывались отсчёты — дискретные моменты времени, в которые измерялся сигнал. Здесь мы видим, как изменяется амплитуда сигнала во времени.
Частотная область — это представление сигнала в зависимости от частоты. На графиках 5.1 и 5.2 по горизонтальной оси откладываются уже не отсчёты, а бины — дискретные частотные интервалы, каждый из которых соответствует определённой частоте в герцах. Здесь мы видим, какие частоты и с какой амплитудой присутствуют в сигнале.

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

В нашем конкретном примере номер бина численно совпадает со значением частоты: бин 1 соответствует частоте 1 Гц, бин 2 — частоте 2 Гц и так далее. Однако это частный случай, обусловленный выбором параметров дискретизации. В общем случае масштаб частотной оси может быть иным: например, бин 1 может соответствовать частоте 1 кГц, бин 2 — 2 кГц и т. д. Соотношение между номером бина и его частотой определяется частотой дискретизации и количеством отсчётов согласно формуле:

$f_k = k \cdot \frac{f_s}{N}$

где — номер бина, — соответствующая ему частота.

На фазовом спектре показан сдвиг фазы $\varphi[k]$ для каждой частотной составляющей. В нашем примере фаза косинуса сдвинута на $45^\circ$ , поэтому на графике для частоты 1 Гц будет отображено значение $\varphi = 45^\circ$ (или $\pi/4$ рад). Для всех остальных частот, где амплитуда близка к нулю, фаза не имеет физического смысла и обычно не отображается или принимается равной нулю.

Эти два графика становятся особенно полезными при анализе сложных сигналов, состоящих из набора простых гармонических компонент. Например, если сигнал представляет собой сумму косинусов с частотами 1 Гц, 3 Гц и 5 Гц и разными амплитудами и фазами, то:

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

Таким образом, пара спектров — частотный и фазовый — даёт полное описание частотного состава сигнала: что в нём есть (частоты), насколько сильно это выражено (амплитуды) и с каким сдвигом во времени начинается каждая компонента (фазы). Это и есть основа спектрального анализа.

Расчёт амплитуды по формуле (1.10)

В примере исследуемый сигнал — косинус с амплитудой , частотным индексом и фазовым сдвигом $\varphi = 45^\circ$ . После вычисления корреляций с базисными функциями получаем:

$\operatorname{Re} X[1] \approx 35.36, \quad \operatorname{Im} X[1] \approx 35.36$ $\operatorname{Mag}[1] = \sqrt{35.36^2 + 35.36^2} \approx 50$

Полученное значение не соответствует физической амплитуде исходного сигнала (), отображённой на графиках.

Применение нормировки

Частотный индекс	Компонента	Коэффициент	Почему
	Постоянная составляющая (среднее значение)	$\frac{1}{N}$	Не имеет симметричной пары
$k = 1 \dots \frac{N}{2}-1$	Основные частоты	$\frac{2}{N}$	Энергия «собирается» из пары $\pm k$
$k = \frac{N}{2}$	Максимальная частота (теорема Котельникова)	$\frac{1}{N}$	Не имеет симметричной пары

табл. 1

Постоянная составляющая (DC) ДПФ — нулевая частотная компонента, равная сумме отсчётов сигнала. При наличии, например, постоянного напряжения на графиках она отображается горизонтальной прямой, а не гармоническим сигналом.

Для приведения спектра к физическим единицам используется коэффициент из таблицы нормировки (табл. 1). Поскольку относится к основным частотам ( $1 \leq k \leq N/2-1$ ), применяем множитель :

$A[1] = \operatorname{Mag}[1] \cdot \frac{2}{N} = 50 \cdot \frac{2}{100} = 1$

Результат: нормированная амплитуда точно совпадает с амплитудой исходного сигнала.

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

Наглядный пример: сдвиг фазы

В примере 6 рассматривался косинус, сдвинутый на . В анимации 1 показан косинус со сдвигом фазы, плавно изменяющимся от до .

На первом графике анимации видно, как постепенно изменяется положение сигнала. На третьем графике показана корреляция этого сигнала с базисными функциями (косинусом и синусом) — то есть результаты скалярных произведений. На втором графике отображаются фаза и амплитуда сигнала. Это то же самое, что и в примере 6, но не расчёт одного конкретного сдвига фазы, а отслеживание сдвига от до . Для наглядности значения $\operatorname{Re} X$ и $\operatorname{Im} X$ нормированы, поэтому и корреляция, и $\operatorname{Mag}$ лежат в пределах от −1 до 1.

Из анимации видно: при увеличении сдвига фазы вектор на комплексной плоскости вращается против часовой стрелки. Если в формуле (1.9) убрать знак «минус», то при увеличении сдвига фазы вектор вращался бы по часовой стрелке. В математике и ЦОС стандартом принято вращение против часовой стрелки. Забегая вперёд, знак важен также в обратном преобразовании Фурье.

На анимации 1 показаны три графика — их можно объединить в один трёхмерный график, что и сделано в анимации 2. На этом 3D-графике видны корреляции действительной ( $\operatorname{Re} X$ ) и мнимой ( $\operatorname{Im} X$ ) составляющих, а вид сбоку представляет собой комплексную плоскость. Пример 6, анимация 1 и анимация 2 — это, по сути, три разных способа показать одно и то же явление. Отличие анимации 2 в том, что изображено два периода сигнала, а не один: в 3D с двумя периодами получается более наглядно.

Вывод

Рассмотренные примеры и анимации наглядно демонстрируют ключевой принцип частотного анализа гармонических сигналов. При сдвиге фазы исследуемого сигнала его корреляция с каждой из базисных функций — косинусом (действительная часть $\operatorname{Re} X$ ) и синусом (мнимая часть $\operatorname{Im} X$ ) — изменяется по отдельности: одна составляющая возрастает, другая убывает, следуя гармоническому закону. Однако в совокупности эти две проекции полностью и однозначно описывают сигнал.

Модуль комплексного числа, согласно формуле (1.10),

$\operatorname{Mag} = \sqrt{(\operatorname{Re} X)^2 + (\operatorname{Im} X)^2}$

остаётся инвариантным относительно сдвига фазы — он всегда равен амплитуде сигнала, независимо от того, насколько сигнал сдвинут во времени. Сам же сдвиг фазы однозначно определяется аргументом комплексного числа, согласно формуле (1.11):

$\varphi = \arctan\left(\frac{\operatorname{Im} X}{\operatorname{Re} X}\right)$

Таким образом, именно совместное рассмотрение скалярных произведений с двумя ортогональными базисными функциями (косинусом и синусом) позволяет разложить гармонический сигнал на два независимых параметра — амплитуду и фазу. Это и составляет математическую основу преобразования Фурье: любой гармонический компонент произвольной амплитуды и произвольного сдвига фазы может быть полностью описан парой действительных чисел $(\operatorname{Re} X, \operatorname{Im} X)$ , или, что эквивалентно, парой $(\operatorname{Mag}, \varphi)$ .

Поиск неизвестного сигнала

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

Но на практике ситуация обратная: параметры сигнала неизвестны, и наша задача — определить их по имеющимся отсчётам. Как это сделать?

Принцип последовательного перебора частот

Решение основано на том же подходе, что и в предыдущих примерах: мы вычисляем корреляцию исследуемого сигнала с набором базисных функций. Разница лишь в том, что теперь мы не знаем, с какой именно частотой нужно сравнивать сигнал.

Поэтому действуем системно:

Задаём набор частотных индексов $k = 0, 1, 2, 3, \dots$
Для каждого значения формируем пару базисных функций: $\cos\left(2\pi \cdot k \cdot \frac{n}{N}\right)$ и $\sin\left(2\pi \cdot k \cdot \frac{n}{N}\right)$ .
Вычисляем корреляцию исследуемого сигнала с каждой из этих функций по формулам:

$\operatorname{Re} X[k] = \sum_{n=0}^{N-1} x[n] \cdot \cos\left(2\pi \cdot k \cdot \frac{n}{N}\right) \quad (1.12)$ $\operatorname{Im} X[k] = -\sum_{n=0}^{N-1} x[n] \cdot \sin\left(2\pi \cdot k \cdot \frac{n}{N}\right) \quad (1.13)$

Анализируем результат: если обе проекции $\operatorname{Re} X[k]$ и $\operatorname{Im} X[k]$ близки к нулю — частота в сигнале отсутствует; если хотя бы одна проекция значима — частота присутствует.

Пример: сигнал из примера 6

Возьмём сигнал, который мы уже анализировали: гармоническое колебание частотой 1 Гц. При последовательном переборе частот $k = 0, 1, 2, 3, \dots$ получим:

При : корреляция с базисными функциями даст ненулевые значения $\operatorname{Re} X(1)$ и $\operatorname{Im} X(1)$ $\rightarrow$ частота 1 Гц обнаружена.
При $k = 0, 2, 3, \dots$ : корреляция будет близка к нулю $\rightarrow$ эти частоты отсутствуют в сигнале.

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

Принцип остаётся прежним: мы так же последовательно перебираем частоты и вычисляем корреляцию. Разница в том, что теперь ненулевой отклик будет наблюдаться не для одной, а для нескольких частот — ровно тех, которые входят в состав сложного сигнала.

Бесконечно увеличивать частоты базисных функций не имеет смысла. Существует фундаментальный предел: максимальная частота, которую можно корректно проанализировать, не может превышать половину частоты дискретизации.

Это ограничение описывается теоремой Котельникова:

$f_{\max} \leq \frac{f_d}{2}$

где:

$f_{\max}$ — максимальная частота в сигнале, которую можно определить;
— частота дискретизации (количество отсчётов в секунду).

В наших примерах частота дискретизации составляет 100 Гц (100 отсчетов исследуемого сигнала). Значит, можно определить сигналы частоты которых лежат в пределах от 0 до 50 Гц. Подробный разбор этой темы выходит за рамки первой части статьи и заслуживает отдельного рассмотрения.

Благодаря ортогональности базисных функций, каждая частота анализируется независимо от остальных: наличие одной составляющей не искажает результат поиска другой. Это ключевое свойство, позволяющее «разделять» сложный сигнал на простые компоненты.

Обобщение: формула для нахождения спектра частот

Обобщая вышесказанное, можно получить формулу дискретного преобразования Фурье (ДПФ). ДПФ в тригонометрической форме для последовательности длины выглядит следующим образом:

$X[k] = \sum_{n=0}^{N-1} x[n] \left( \cos\left(\frac{2\pi k n}{N}\right) - j \sin\left(\frac{2\pi k n}{N}\right) \right) \quad (1.14)$

Вынесем знак «минус» из правой части формулы (1.13) и перенесем его в итоговое выражение перед мнимой единицей j Получим:

$X[k] = \operatorname{Re} X[k] - j \cdot \operatorname{Im} X[k] \quad (1.14.1)$

где — комплексный спектр, — мнимая единица. Действительная и мнимая части могут быть записаны отдельно в виде формул (1.12) и (1.13).

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

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

Расшифровка обозначений:

— комплексный спектр сигнала на дискретной частоте с индексом . Это выход ДПФ.
— индекс частотной компоненты ( $k = 0, 1, \dots, N-1$ ). Пропорционален частоте: $\omega_k = 2\pi k / N$ (в радианах на отсчёт).
$\sum_{n=0}^{N-1}$ — сумма по всем временным отсчётам от до .
— исходный дискретный сигнал (вещественный или комплексный) в момент времени с индексом .
— индекс временного отсчёта ( $n = 0, 1, \dots, N-1$ ).
— длина сигнала (количество отсчётов); также равно количеству частотных компонент.
$\cos\left(\frac{2\pi k n}{N}\right)$ — базисная косинусная функция.
$\sin\left(\frac{2\pi k n}{N}\right)$ — базисная синусная функция.
— мнимая единица (). В некоторых областях используют .
Знак «» перед $j \sin(\dots)$ — стандартное соглашение для прямого ДПФ (соответствует комплексной экспоненте $e^{-j\theta}$ ).

Переход к экспоненциальной форме

Формулу (1.14) можно сделать короче и элегантнее, если вспомнить формулу Эйлера:

$e^{-j\theta} = \cos\theta - j \cdot \sin\theta \quad (1.15)$

Это одно из самых красивых равенств в математике. Оно связывает комплексную экспоненту с тригонометрическими функциями. Если подставить в это выражение $\theta = 2\pi \cdot k \cdot \frac{n}{N}$ , то наша громоздкая формула превратится в одну компактную строчку:

$X[k] = \sum_{n=0}^{N-1} x[n] \cdot e^{-j \cdot 2\pi \cdot k \cdot \frac{n}{N}} \quad (1.16)$

Это и есть экспоненциальная форма ДПФ. Она не содержит ни синусов, ни косинусов по отдельности — вся информация о частоте и фазе «упакована» в комплексную экспоненту $e^{-j\theta}$ .

Зачем нужна экспоненциальная форма?

Переход к экспоненте — это не просто «математический трюк». Это даёт важные преимущества:

Краткость и единообразие. Одна формула заменяет две (для синуса и косинуса). Все вычисления становятся компактнее, особенно при выводе теории.
Удобство анализа. Экспоненциальная форма тесно связана с преобразованием Фурье в непрерывном мире. Это делает математические выкладки (например, доказательство свойств ДПФ, таких как линейность или сдвиг) значительно проще и изящнее.

В следующем примере мы сгенерируем сигнал из четырёх гармоник с разными частотами, амплитудами и фазами и покажем, как метод последовательной корреляции позволяет точно восстановить все параметры каждой составляющей.

Пример 7

Анализ сложного сигнала (сумма нескольких гармоник)

До сих пор мы рассматривали сигналы, состоящие из одной гармонической составляющей. Параметры такого сигнала (частота, амплитуда, фаза) были нам известны заранее, так как мы сами их задавали при генерации. Однако реальные сигналы почти всегда содержат несколько частотных составляющих одновременно.

Как получить такой сложный сигнал? Математически это делается путём поэлементного сложения нескольких гармонических сигналов.

Алгоритм сложения:

Берём нулевые отсчёты всех компонент: $x_1[0], x_2[0], x_3[0], \dots$ и складываем их $\rightarrow$ получаем $x_{\text{сум}}[0]$ .
Берём первые отсчёты: $x_1[1], x_2[1], x_3[1], \dots$ и складываем $\rightarrow$ получаем $x_{\text{сум}}[1]$ .
Продолжаем этот процесс для всех $n = 0, 1, \dots, N-1$ .

Параметры сигналов

Тип	Фаза $\varphi$	Формула
DC (пост.)	$0^\circ$
$\cos$	$250^\circ$	$x_2[n] = 2.0 \cdot \cos\left(2\pi \cdot \frac{5 \cdot n}{N} + 250^\circ\right)$
$\cos$	$135^\circ$	$x_3[n] = 3.0 \cdot \cos\left(2\pi \cdot \frac{10 \cdot n}{N} + 135^\circ\right)$
$\sin$	$135^\circ$	$x_4[n] = 1.5 \cdot \sin\left(2\pi \cdot \frac{15 \cdot n}{N} + 135^\circ\right)$

где отсчётов, $n = 0, 1, \dots, N-1$ .

Визуализация и интерпретация (рис. 9)

На рис. 9 представлены графики отдельных компонент (1–4) и результирующий сигнал (5), полученный их суммированием. Ниже показаны итоговые частотный и фазовый спектры.

Почему здесь нет детальных графиков поэлементного умножения и накопленной суммы, как в предыдущих примерах? Дело в том, что в примерах 1–5 мы фокусировались на одной конкретной частоте, чтобы детально показать механизм скалярного произведения. В данном же примере мы применяем алгоритм ДПФ, анализируя спектр частот (в нашем случае от до Гц). Отображать 50 наборов промежуточных графиков было бы избыточно и затруднило бы восприятие общей картины.

Результирующий сигнал проходит те же этапы математических расчётов, что и в примере 5, но с одним ключевым отличием: базисные функции последовательно «сканируют» сигнал, меняя частотный индекс от до . Именно это последовательное вычисление корреляции и описывают формулы прямого ДПФ (1.14) или (1.16).

Вывод: линейность и независимость частот

Для алгоритма ДПФ не имеет значения, анализируем мы гармоники по отдельности или их сумму. Благодаря корреляции и ортогональности базисных функций, каждая частота извлекается независимо от остальных. Корреляция с базисом «вытащит» нужную составляющую, точно определив её амплитуду и фазу, словно других компонентов в сигнале в этот момент вообще не существует.

Именно это свойство делает преобразование Фурье таким мощным инструментом: оно позволяет анализировать сложные сигналы, раскладывая их на простые, понятные нам составляющие.

Важное уточнение: такая идеальная независимость гарантирована только при целых частотных индексах. Если бы хотя бы одна из четырёх гармоник имела дробную частоту, её энергия начала бы «размазываться» по соседним частотам, искажая результаты измерения. Эту проблему мы уже затрагивали в примере 3, и она заслуживает отдельного внимания (явление «утечки спектра» может быть разобрано в рамках одной из следующих частей статьи).

Форматы представления фазы: от к

В математике и программировании фазу часто измеряют в пределах полного круга: от $0^\circ$ до $360^\circ$ (или от до $2\pi$ ). Однако в цифровой обработке сигналов (ЦОС) стандартом де-факто является диапазон $-180^\circ\dots180^\circ$ (или $-\pi\dots\pi$ ).

Для удобства восприятия и перехода между этими системами координат воспользуемся таблицей 2:

Полный круг ( $0\dots360^\circ$ )	$-180^\circ\dots180^\circ$	Радианы	Физический смысл в ЦОС
$0^\circ$	$0^\circ$		Совпадение фаз (базовый косинус)
$90^\circ$	$90^\circ$	$\pi/2$	Опережение на периода
$180^\circ$	$\pm 180^\circ$	$\pm\pi$	Полная противофаза
$270^\circ$	$-90^\circ$	$-\pi/2$	Отставание на периода

табл. 2

Математический алгоритм перевода

Чтобы перевести угол из диапазона $0\dots360^\circ$ в «инженерный» диапазон $-180^\circ\dots180^\circ$ , используется простое кусочное правило:

$\varphi_{\text{ЦОС}} =\begin{cases}\varphi, & \text{если } 0^\circ \le \varphi \le 180^\circ \\\varphi - 360^\circ, & \text{если } 180^\circ < \varphi < 360^\circ\end{cases}$

Разберём на примерах:

$\varphi = 45^\circ \implies \mathbf{45^\circ}$ (остаётся без изменений, так как $\le 180^\circ$ ).
$\varphi = 180^\circ \implies \mathbf{\pm 180^\circ}$ (граница диапазона, оба варианта физически эквивалентны).
$\varphi = 250^\circ \implies 250^\circ - 360^\circ = \mathbf{-110^\circ}$ .
$\varphi = 359^\circ \implies 359^\circ - 360^\circ = \mathbf{-1^\circ}$ .

💡 Важное замечание для программистов: Если вы пишете код на Python (например, используя numpy.arctan2 или cmath.phase для вычисления фазы), библиотеки автоматически возвращают результат именно в диапазоне $-\pi\dots\pi$ (или $-180^\circ\dots180^\circ$ ).

Зачем нужен диапазон ?

Этот формат показывает кратчайшее отклонение от нулевой фазы. Главное практическое преимущество — устранение искусственных вертикальных скачков на графиках фазового спектра. Если бы мы использовали диапазон $0\dots360^\circ$ , то значения $359^\circ$ и $1^\circ$ оказались бы на противоположных концах графика, хотя физически это почти одно и то же колебание (разница всего в $2^\circ$ ). Диапазон $-\pi\dots\pi$ делает анализ сдвигов интуитивно понятным и непрерывным.

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

Результаты преобразования Фурье для нашего сложного сигнала представлены на рис. 9. В итоговой таблице сведены восстановленные характеристики (частота, амплитуда и фаза) для каждой значимой составляющей.

Давайте подробно разберём каждую строку и сопоставим её с исходными параметрами, которые мы задавали при генерации.

Сигнал 1: постоянная составляющая (DC, )

Задано: , $\varphi = 0^\circ$
Найдено ДПФ: , $\varphi = 0^\circ$
Вывод: все параметры полностью совпадают. Постоянная составляющая (смещение графика вверх) корректно распознана алгоритмом.

Сигнал 2: косинус ()

Задано: , $\varphi = 250^\circ$
Найдено ДПФ: , $\varphi = -110^\circ$
Вывод: амплитуда и частота совпадают. Фаза $250^\circ$ автоматически переведена алгоритмом в стандарт ЦОС. Применяем математическое правило перевода (так как $250^\circ > 180^\circ$ ): $250^\circ - 360^\circ = \mathbf{-110^\circ}$ . Значение абсолютно корректно.

Сигнал 3: косинус ()

Задано: , $\varphi = 135^\circ$
Найдено ДПФ: , $\varphi = 135^\circ$
Вывод: идеальное совпадение. Исходная фаза уже лежала в допустимом диапазоне ( $0^\circ\dots180^\circ$ ), поэтому преобразований не потребовалось.

Сигнал 4: синус и «неожиданный» сдвиг фазы ()

Задано: , $\varphi = 135^\circ$ (тип: синус)
Найдено ДПФ: , $\varphi = \mathbf{45^\circ}$
Вывод: амплитуда и частота совпадают, но фаза отличается от заданной на $90^\circ$ . В чём причина?

Секрет сигнала 4

Дело в том, что классический алгоритм ДПФ «смотрит» на мир через призму косинуса. Фаза в ДПФ всегда отсчитывается строго от оси $\operatorname{Re}$ комплексной плоскости.

Поскольку наш 4-й сигнал задан как синус, алгоритм математически сводит его к эквивалентному косинусу, используя фундаментальное тригонометрическое тождество:

$\sin(\alpha) = \cos(\alpha - 90^\circ) \quad (1.17)$

Применим это правило к нашему сигналу:

$x_4[n] = 1.5 \cdot \sin\left(\dots + 135^\circ\right) \implies 1.5 \cdot \cos\left(\dots + 135^\circ - 90^\circ\right)$ $135^\circ - 90^\circ = \mathbf{45^\circ}$

Таким образом, фаза $45^\circ$ , полученная в результате ДПФ — это фаза эквивалентного косинуса. Алгоритм не ошибся, он просто честно перевёл синусоидальный сигнал в косинусоидальный базис, как того требует математика преобразования Фурье.

Краткий вывод по интерпретации фазы

ДПФ «видит мир» строго через косинус (базис $\operatorname{Re}$ ). Алгоритм ДПФ всегда сравнивает анализируемый сигнал с базисным косинусом (осью $\operatorname{Re}$ комплексной плоскости). Если исходный сигнал задан как синус, ДПФ автоматически сведёт его к эквивалентному косинусу.
- Формула: $\sin(\alpha) = \cos(\alpha - 90^\circ)$
- Следствие: фаза синусоиды в спектре всегда будет смещена ровно на $-90^\circ$ относительно её исходной «синусной» фазы.
Стандартный диапазон: от $-180^\circ$ до $+180^\circ$ . Все современные математические библиотеки (например, функция arctan2 в Python) возвращают фазу в диапазоне $-\pi \dots \pi$ (или $-180^\circ \dots 180^\circ$ ). Это сделано для того, чтобы показывать кратчайшее отклонение от нуля и избегать разрывов на графиках.
Перевод ( $0\dots360^\circ$ ). Перевод выполняется в уме за одно действие по простому условию:
- Если угол $\varphi \le 180^\circ$ — оставляем его без изменений.
- Если угол $\varphi > 180^\circ$ — вычитаем из него $360^\circ$ .

С формулами из примеров на Python можно ознакомиться в следующем разделе. Разобран 7-й пример, все его формулы и алгоритмы на Python.

Перевод формул (1.12), (1.13), (1.14) в код на Python

import math# ============================================================================# ДПФ ПРЯМОЕ В ТРИГОНОМЕТРИЧЕСКОЙ ФОРМЕ# ============================================================================def dft_smith(x):    """    Прямое вычисление ДПФ (Дискретного Преобразования Фурье).            📖 Связь со статьёй: Формулы (1.12), (1.13), (1.14)        X[k] = Σ_{n=0}^{N-1} x[n] · cos(2πkn/N)  -  j · Σ_{n=0}^{N-1} x[n] · sin(2πkn/N)             └─────────────┬───────────────────┘   └────────────┬───────────────────────┘                   Действительная часть                  Мнимая часть    📐 ФИЗИЧЕСКИЙ СМЫСЛ:        ДПФ последовательно сравнивает сигнал x[n] с двумя эталонными         тригонометрическими гармониками частоты k:        • Косинусная волна cos(2πkn/N)  → "фазовый ноль" (in-phase)        • Синусная волна sin(2πkn/N)    → "фазовый сдвиг -90°" (quadrature)        Результат сравнения (скалярное произведение) даёт два числа:        re_X[k] = Σ x[n]·cos(θ)  → насколько сигнал похож на косинус        im_X[k] = -Σ x[n]·sin(θ) → насколько сигнал похож на синус (со знаком минус)    📦 ВОЗВРАЩАЕТ:        (re_X, im_X) — списки действительной и мнимой частей спектра.    """        # ========================================================================    # 1: ОПРЕДЕЛЕНИЕ ДЛИНЫ СИГНАЛА    # ========================================================================    # len(x) — возвращает общее количество отсчётов (N) в анализируемом сигнале.    N_val = len(x)        # ========================================================================    # 2: ИНИЦИАЛИЗАЦИЯ МАССИВОВ ДЛЯ НАКОПЛЕНИЯ ПРОЕКЦИЙ    # ========================================================================    # Создаём два списка длиной N_val, заполненных нулями.    # re_X[k] будет накапливать проекцию сигнала на косинусный базис частоты k.    # im_X[k] будет накапливать проекцию на синусный базис (с учётом знака ДПФ).    re_X = [0.0] * N_val    im_X = [0.0] * N_val        # ========================================================================    # 3: ВНЕШНИЙ ЦИКЛ — ПЕРЕБОР ЧАСТОТНЫХ ГАРМОНИК [k]    # ========================================================================    # 📖 Связь со статьёй: Принцип последовательного перебора частот    # Цикл перебирает все возможные частотные индексы от 0 до N-1:    # • k = 0   → постоянная составляющая (среднее значение сигнала, DC)    # • k = 1   → основная частота (один полный период на длине N)    # • k = N/2 → максимальная частота (по теореме Котельникова)    for k in range(N_val):                # Обнуляем накопители перед анализом новой гармоники,         # чтобы суммы не "перетекали" с предыдущей частоты.        re_sum = 0.0        im_sum = 0.0                # ====================================================================        # 4: ВНУТРЕННИЙ ЦИКЛ — ПЕРЕБОР ВРЕМЕННЫХ ОТСЧЁТОВ (n)        # ====================================================================        # Используем n_idx, чтобы не затенять внешние переменные и явно         # указать, что это индекс текущего временного отсчёта.        for n_idx in range(N_val):                        # ----------------------------------------------------------------            # 5: ВЫЧИСЛЕНИЕ МГНОВЕННОГО УГЛА ЭТАЛОННОЙ ГАРМОНИКИ            # ----------------------------------------------------------------            # 📖 Связь со статьёй: Аргумент тригонометрических функций            # Формула: θ = 2π · k · n / N            # Это мгновенный фазовый угол базисной функции для текущего отсчёта.            angle = 2.0 * math.pi * k * n_idx / N_val                        # ----------------------------------------------------------------            # 6: РАЗЛОЖЕНИЕ ПО ТРИГОНОМЕТРИЧЕСКИМ БАЗИСАМ            # ----------------------------------------------------------------            # 📖 Связь со статьёй: Формулы (1.14) и (1.16)            # ДПФ вычисляет два скалярных произведения одновременно.                        # 6.1. ПРОЕКЦИЯ НА КОСИНУС (Действительная часть):            # re_sum += x[n] · cos(θ)            # Если сигнал x[n] совпадает по форме с косинусом частоты k,            # эта сумма будет максимальной. Если сигналы ортогональны → ~0.            re_sum += x[n_idx] * math.cos(angle)                        # 6.2. ПРОЕКЦИЯ НА СИНУС (Мнимая часть) ⚠️ КЛЮЧЕВОЙ ЗНАК "МИНУС":            # im_sum -= x[n] · sin(θ)              # Если убрать этот минус, то при изменении фазы сигнала вектор             # на комплексной плоскости вращался бы по часовой стрелке.             # В математике и ЦОС стандартом принято вращение против часовой стрелки.            #            # 📐 ТРИГОНОМЕТРИЧЕСКАЯ ИНТЕРПРЕТАЦИЯ:            # • Если x[n] = cos(θ)  → re_sum > 0, im_sum = 0  → φ = 0°            # • Если x[n] = sin(θ)  → re_sum = 0, im_sum < 0  → φ = -90°            # • Если x[n] = -sin(θ) → re_sum = 0, im_sum > 0  → φ = +90°            # Без этого минуса знаки фаз для синусоид инвертировались бы.            im_sum -= x[n_idx] * math.sin(angle)                # ====================================================================        # 7: СОХРАНЕНИЕ ИТОГОВЫХ ПРОЕКЦИЙ ДЛЯ k-й ГАРМОНИКИ        # ====================================================================        # После прохождения всех отсчётов n, накопленные суммы записываются         # в соответствующие ячейки выходных массивов.        re_X[k] = re_sum        im_X[k] = im_sum            # ========================================================================    # 8: ВОЗВРАТ РЕЗУЛЬТАТА    # ========================================================================    # 📖 Связь со статьёй: Формулы (1.10) и (1.11)    # Возвращаем спектр в виде двух параллельных списков.    # Далее в коде (в функции get_dft_spectrum) они преобразуются в:    # • Амплитуду:   mag   = sqrt(re² + im²)    # • Фазу:        phase = degrees(atan2(im, re))    return re_X, im_X

Перевод формул (1.10), (1.11) в код на Python

import numpy as npimport math# ============================================================================# СПЕКТР АМПЛИТУД И ФАЗ (ДЛЯ ГРАФИКОВ)# ============================================================================def get_dft_spectrum(re_X, im_X, return_physical_amp=False, N_signal=None,                      phase_threshold_ratio=0.01, unwrap_phase=False,                      return_half_spectrum=True):    """        Амплитуда гармоники k:  A[k] = √(Re[k]² + Im[k]²)        Фаза гармоники k:       φ[k] = arctan2(Im[k], Re[k])                Физическая нормировка (для действительного сигнала):        • A[0] (постоянная составляющая) и A[N/2] (частота Найквиста) делятся на N.        • Все остальные гармоники (1 ≤ k < N/2) делятся на N/2 (или умножаются на 2/N).    """    N = len(re_X)        # Преобразуем входные списки в массивы numpy для векторных вычислений    re = np.array(re_X, dtype=float)    im = np.array(im_X, dtype=float)        # ========================================================================    # 1. АМПЛИТУДНЫЙ СПЕКТР    # ========================================================================    # 📖 Связь со статьёй:: Формула (1.10)    # Mag[k] = √(Re X[k]² + Im X[k]²)    # np.sqrt() — извлекает квадратный корень из каждого элемента массива    amplitudes = np.sqrt(re**2 + im**2)        # ========================================================================    # 2. НОРМИРОВКА НА ФИЗИЧЕСКУЮ АМПЛИТУДУ    # ========================================================================    # 📖 Связь со статьёй: Таблица 1 (Таб. 1) и раздел "Расчёт амплитуды..."    #    # "Сырые" значения ДПФ необходимо масштабировать, чтобы получить реальные     # физические амплитуды исходного сигнала.        if return_physical_amp:        if N_signal is None:            N_signal = N                    # for k in range(N): — это цикл, который перебирает все частоты от 0 до N-1.        for k in range(N):                        # 2.1: Проверка специальных частот            # Условие: k == 0 or (N % 2 == 0 and k == N // 2)            if k == 0 or (N % 2 == 0 and k == N // 2):                # /= N_signal — это сокращённая запись: amplitudes[k] = amplitudes[k] / N_signal                # Делим текущую амплитуду на N (количество отсчётов).                amplitudes[k] /= N_signal                            # 2.2: Все остальные частоты            # Для частот от 1 до N/2-1  применяется коэффициент 2/N.            else:                # Делим на (N/2), что математически эквивалентно умножению на 2/N.                amplitudes[k] /= (N_signal / 2)        # ========================================================================    # 3. ФАЗОВЫЙ СПЕКТР С МАСКИРОВКОЙ ШУМА    # ========================================================================    # np.zeros(N): Создаёт новый список (массив) из N элементов,     # заполненный нулями. Это "пустая заготовка" под будущие значения фаз.    phases = np.zeros(N)        # np.max(amplitudes): Находит самое большое число в списке амплитуд.    # Мы умножаем его на 0.01 (1%), чтобы получить порог (amp_threshold).    amp_threshold = phase_threshold_ratio * np.max(amplitudes)        for k in range(N):        if amplitudes[k] < amp_threshold:            # np.nan (Not a Number): Специальная метка "Нет данных" или "Неопределено".            # Мы ставим её вместо фазы, чтобы на графике не рисовался случайный "мусор"             # там, где амплитуда сигнала практически равна нулю.            phases[k] = np.nan          else:            # np.arctan2(im[k], re[k]): "Умный" арктангенс.             # В отличие от обычного arctan(y/x), он корректно определяет угол             # во всех 4-х четвертях круга и не ломается при делении на ноль.            # np.degrees(...): Переводит результат из радиан (стандарт математики)             # в привычные нам градусы (стандарт человека).            phases[k] = np.degrees(np.arctan2(im[k], re[k]))        # ========================================================================    # 4. РАЗВЁРТКА ФАЗЫ (UNWRAP)    # ========================================================================    # 📖 Связь со статьёй: Таблица 2 (Таб. 2).    # Устраняет искусственные скачки фазы при переходе через границу ±180°.        if unwrap_phase:        # 4.1: np.radians(phases)        # Переводит список фаз из градусов обратно в радианы.         # Функция "развёртки" (unwrap) понимает только радианы.        phases_rad = np.radians(phases)                # 4.2: np.isnan(phases_rad)        # "Детектор дырок". Проверяет каждый элемент списка и возвращает         # True (если это np.nan) или False (если это нормальное число).                # 4.3: np.where(условие, значение_если_ДА, значение_если_НЕТ)        # Это аналог конструкции "ЕСЛИ ... ТО ... ИНАЧЕ ...", применённый ко всему списку сразу.        # Читаем так: "ЕСЛИ элемент является np.isnan (дыркой), ТО замени его на 0,         # ИНАЧЕ оставь исходное значение (phases_rad)".        # Зачем? Функция unwrap не умеет работать с "дырками" (np.nan) и выдаст ошибку.         # Мы временно "замазываем" их нулями, чтобы алгоритм отработал.        phases_rad_no_nan = np.where(np.isnan(phases_rad), 0, phases_rad)                # 4.4: np.unwrap(phases_rad_no_nan)        # Главная функция "развёртки". Она идёт по списку слева направо и смотрит         # на разницу между соседними числами. Если она видит резкий скачок         # (например, с +179° сразу на -179°, что в радианах выглядит как скачок через границу π),         # она понимает: "Это не настоящий скачок, просто счётчик переполнился!".         # И автоматически добавляет или вычитает 360° (2π), чтобы линия на графике стала плавной.        # Аналогия: одометр в машине, который вместо сброса с 999 на 000 продолжает считать: 1000, 1001...        phases_unwrapped = np.unwrap(phases_rad_no_nan)                # 4.5: Возврат "дырок" и перевод в градусы        # Снова используем np.where, но с обратной логикой:        # "ЕСЛИ в исходном списке (phases_rad) тут была дырка (np.isnan),         # ТО верни дырку (np.nan) на место, ИНАЧЕ возьми сглаженное значение из phases_unwrapped".        # Внешняя функция np.degrees(...) переводит итоговый сглаженный список обратно в градусы.        phases = np.degrees(np.where(np.isnan(phases_rad), np.nan, phases_unwrapped))        # ========================================================================    # 5. МАССИВ ИНДЕКСОВ ЧАСТОТ    # ========================================================================    # np.arange(N) — создаёт массив целых чисел [0, 1, 2, ..., N-1]    freq_indices = np.arange(N)        # ========================================================================    # 6. ВОЗВРАТ ТОЛЬКО НИЖНЕЙ ПОЛОВИНЫ СПЕКТРА    # ========================================================================        if return_half_spectrum:        half_N = N // 2 + 1        # [:half_N] — срез массива (берёт первые half_N элементов)        return amplitudes[:half_N], phases[:half_N], freq_indices[:half_N]        return amplitudes, phases, freq_indices

Перевод формул (1.4), (1.5) в код на Python

import numpy as np# ============================================================================# ГЕНЕРАЦИЯ ДИСКРЕТНЫХ ГАРМОНИЧЕСКИХ СИГНАЛОВ# ============================================================================def generate_signal(amp, freq, phase_deg, signal_type, n):    """    Генерирует дискретный гармонический сигнал (синус или косинус) для тестирования ДПФ.        📘 МАТЕМАТИЧЕСКАЯ ОСНОВА:        📖 Связь со статьёй: Формулы (1.4) и (1.5)        x[n] = A · cos(2π · (k/N) · n + φ)   или   x[n] = A · sin(2π · (k/N) · n + φ)                Где:        • A (amp) — амплитуда сигнала.        • k (freq) — частотный индекс (количество полных периодов на длине N).        • N — общее количество отсчётов.        • n — индекс текущего отсчёта (дискретное время).        • φ (phase_deg) — начальная фаза в градусах.    📐 ФИЗИЧЕСКИЙ СМЫСЛ:        Функция создаёт "идеальную" эталонную гармонику, которая укладывается         целым числом периодов в интервал наблюдения N. Это критически важно         для демонстрации ортогональности базисных функций в ДПФ и отсутствия         эффекта "утечки спектра" (spectral leakage).    📦 ПАРАМЕТРЫ:        • amp (float) — амплитуда сигнала (A).        • freq (float) — частотный индекс [k], число периодов на длине N.        • phase_deg (float) — начальная фаза в градусах (φ).        • signal_type (str) — тип базовой функции: 'sin' или 'cos'.        • n (array_like) — массив дискретного времени [0, 1, 2, ..., N-1].    📦 ВОЗВРАЩАЕТ:        (signal, name) — кортеж из:        • signal (np.ndarray) — массив значений сигнала длины N.        • name (str) — человекочитаемая строка с формулой для подписи на графике.    """        # ========================================================================    # 1: ОПРЕДЕЛЕНИЕ ДЛИНЫ СИГНАЛА    # ========================================================================    # len(n) возвращает общее количество отсчётов (N) в анализируемом интервале.    N = len(n)    # ========================================================================    # 2: ВЫБОР ТИПА СИГНАЛА И УНИФИКАЦИЯ ЧЕРЕЗ СИНУС    # ========================================================================    # 📖 Связь со статьёй: Раздел "Генерация тестовых сигналов"    # Мы используем математическую хитрость для унификации кода:     # косинус можно представить как синус, сдвинутый на 90 градусов.    # Формула: cos(α) = sin(α + 90°)        if signal_type == 'sin':        # Для синуса фаза применяется напрямую        total_phase_deg = phase_deg        name = f"{amp}·sin(2π·{freq}·n/N"    else:  # signal_type == 'cos'        # Для косинуса добавляем 90 градусов к исходной фазе        total_phase_deg = phase_deg + 90        name = f"{amp}·cos(2π·{freq}·n/N"    # ========================================================================    # 3: ПЕРЕВОД ФАЗЫ ИЗ ГРАДУСОВ В РАДИАНЫ    # ========================================================================    # Тригонометрические функции в Python работают с радианами.    # Формула перевода: радианы = градусы × (π / 180)    # 📌 Мы используем np.pi — это стандартная, максимально точная константа     # числа Пи, встроенная в библиотеку NumPy.    phase_rad = total_phase_deg * np.pi / 180.0    # ========================================================================    # 4: РАСЧЁТ АРГУМЕНТА ТРИГОНОМЕТРИЧЕСКОЙ ФУНКЦИИ    # ========================================================================    # 📖 Связь со статьёй: Аргумент гармонического колебания θ[n]    # Формула: θ[n] = 2π · (freq / N) · n + phase_rad    #     # Где:    # • 2π · (freq / N) — угловая частота, приходящаяся на один отсчёт.    # • n — текущий номер отсчёта (дискретное время).    # • phase_rad — постоянный фазовый сдвиг.    #    # Благодаря векторизации NumPy, эта операция мгновенно вычисляет     # мгновенный фазовый угол для каждого элемента массива n.    argument = 2 * np.pi * (freq / N) * n + phase_rad    # ========================================================================    # 5: ГЕНЕРАЦИЯ МАССИВА ОТСЧЁТОВ СИГНАЛА    # ========================================================================    # Инициализируем массив нулями нужной длины.    signal = np.zeros(N)        # Используем явный цикл для максимальной наглядности процесса генерации     # (в реальных высоконагруженных задачах здесь была бы векторная операция:     # signal = amp * np.sin(argument)).    for i in range(N):        # Унифицированная генерация: мы ВСЕГДА используем функцию np.sin().        # Если на Шаге 2 был выбран 'cos', нужный сдвиг фазы уже учтён         # в переменной total_phase_deg, поэтому формула остаётся единой.        signal[i] = amp * np.sin(argument[i])    # ========================================================================    # 6: ФОРМИРОВАНИЕ ТЕКСТОВОЙ ПОДПИСИ ДЛЯ ГРАФИКА    # ========================================================================    # Дополняем строку name информацией о фазовом сдвиге, если он не равен нулю.    if phase_deg != 0:        # Добавляем знак "+" для положительных значений для красоты формулы        name += f"+{phase_deg}°" if phase_deg > 0 else f"{phase_deg}°"        # Закрываем скобку математического выражения    name += ")"      # ========================================================================    # 7: ВОЗВРАТ РЕЗУЛЬТАТА    # ========================================================================    # Возвращаем кортеж: сам массив данных и его строковое математическое описание.    return signal, name

Пример 7 — применение вышеописанных функций с генерацией графиков

import numpy as npimport matplotlib.pyplot as pltimport os# ============================================================================# 1: ГЕНЕРАЦИЯ ТЕСТОВЫХ СИГНАЛОВ# ============================================================================# 📖 Связь со статьёй: Пример 7, формулы (1.4), (1.5).# Формируем сложный сигнал как сумму 4 гармоник.N = 100  # Длина последовательности (количество отсчётов)n = np.arange(N)  # np.arange: создание массива дискретного времени [0, 1, ..., N-1]SIGNAL_PARAMS = [    {"amp": 1.5, "freq": 0,   "phase": 0,   "type": "cos"},  # DC (постоянная составляющая)    {"amp": 2.0, "freq": 5,   "phase": 250, "type": "cos"},      {"amp": 3.0, "freq": 10,  "phase": 135, "type": "cos"},      {"amp": 1.5, "freq": 15,  "phase": 135, "type": "sin"}   # Синус (важно для проверки фазы)]# generate_signal(): 📖 Определена ранее в статье. # Генерирует "идеальную" эталонную гармонику, укладывающуюся целым числом периодов # в интервал N. Это критически важно для предотвращения утечки спектра и # демонстрации строгой ортогональности базисных функций ДПФ.signals_raw = [generate_signal(p["amp"], p["freq"], p["phase"], p["type"], n)                for p in SIGNAL_PARAMS]signal_arrays = [s[0] for s in signals_raw]  # Извлекаем только массивы значений# np.zeros(N): инициализация нулевого массива, чтобы sum() не модифицировал первый элемент.# sum(): поэлементное сложение массивов для формирования составного сигнала x_sum[n].sig_sum = sum(signal_arrays, np.zeros(N, dtype=float))# ============================================================================# 2: ВЫЧИСЛЕНИЕ ДПФ И ПОЛУЧЕНИЕ СПЕКТРОВ# ============================================================================# 📖 Связь со статьёй: Формулы (1.12), (1.13), (1.14).# dft_smith(): 📖 Определена ранее в статье. # Выполняет прямое вычисление ДПФ в тригонометрической форме (O(N²)). # Последовательно вычисляет скалярные произведения сигнала с базисными # косинусами и синусами, возвращая "сырые" проекции re_X и im_X.re_X, im_X = dft_smith(sig_sum)# get_dft_spectrum(): 📖 Определена ранее в статье. # Выполняет постобработку "сырых" проекций:# • return_physical_amp=True: нормировка амплитуд на N (для k=0, N/2) #   и N/2 (для остальных частот) согласно Таблице 1 статьи.# • phase_threshold_ratio=0.01: маскировка фазового шума (замена на np.nan), #   если амплитуда гармоники менее 1% от максимума.# • return_half_spectrum=True: возврат только первой половины спектра (0...N/2), #   так как вторая половина для вещественного сигнала является зеркальной.amplitudes, phases, freq_indices = get_dft_spectrum(    re_X, im_X, return_physical_amp=True, N_signal=N,    phase_threshold_ratio=0.01, return_half_spectrum=True)# Фильтрация значимых амплитуд для чистоты графиковnonzero = amplitudes > 1e-6freq_nz = freq_indices[nonzero]amp_nz = amplitudes[nonzero]# ============================================================================# !!! ПРИМЕЧАНИЕ: ДАЛЕЕ ИДЁТ ВИЗУАЛИЗАЦИЯ, СОЗДАННАЯ С ПОМОЩЬЮ ИИ !!!# ============================================================================# Блок построения графиков (matplotlib) сгенерирован ИИ.# Математика, вызовы функций и параметры сигналов — авторские.# ============================================================================# ============================================================================# 3: НАСТРОЙКИ СОХРАНЕНИЯ# ============================================================================IMAGES_FOLDER = "Images_part_1"DPI = 300os.makedirs(IMAGES_FOLDER, exist_ok=True)# ============================================================================#  4: ПОДГОТОВКА МАСШТАБА ОСЕЙ# ============================================================================# Общий диапазон по оси Y для сопоставимости всех 5 временных графиков.all_signals = signal_arrays + [sig_sum]max_abs = max(np.max(np.abs(s)) for s in all_signals)ylim = (-max_abs * 1.05, max_abs * 1.05)horiz_levels = [-max_abs, -max_abs/2, 0, max_abs/2, max_abs]x_ticks = np.arange(0, N+1, 10)# ============================================================================# 5: ПОСТРОЕНИЕ ФИГУРЫ (7 ГРАФИКОВ + ТАБЛИЦА)# ============================================================================# 📖 Связь со статьёй: Рис. 9. Компоновка 7×1, стиль и отступы — от ИИ.plt.style.use('seaborn-v0_8-darkgrid' if 'seaborn-v0_8-darkgrid' in plt.style.available else 'default')fig = plt.figure(figsize=(14, 24))# Увеличен нижний отступ, чтобы таблица не перекрывала нижний графикplt.subplots_adjust(hspace=0.45, top=0.96, bottom=0.18)  # ИЗМЕНЕНО: bottom 0.12 → 0.18# --- 5 временных графиков ---for idx in range(5):    ax = fig.add_subplot(7, 1, idx+1)    if idx < 4:        markerline, stemlines, baseline = ax.stem(n, signal_arrays[idx], linefmt=f'C{idx}-', markerfmt='o', basefmt='none')        plt.setp(markerline, markersize=4.0, markeredgewidth=0.5)        plt.setp(stemlines, linewidth=1.2, alpha=0.7)        ax.set_title(f"{idx+1}. Сигнал {idx+1}", fontsize=13, fontweight='bold')    else:        markerline, stemlines, baseline = ax.stem(n, sig_sum, linefmt='-', markerfmt='o', basefmt='none')        plt.setp(stemlines, color='steelblue', linewidth=1.5, alpha=0.8)        plt.setp(markerline, markersize=4.0, markeredgewidth=0.5, color='steelblue')        ax.set_title("5. Результирующий сигнал (сумма)", fontsize=13, fontweight='bold')    ax.set_ylabel("Амплитуда", fontsize=12)    ax.set_xlabel("Отсчёты (n)", fontsize=12)    ax.set_ylim(ylim)    for xv in x_ticks: ax.axvline(x=xv, color='gray', linestyle='--', linewidth=0.8, alpha=0.5)    for yh in horiz_levels: ax.axhline(y=yh, color='gray', linestyle='--', linewidth=0.8, alpha=0.5)    ax.grid(True, alpha=0.5, linestyle='-', linewidth=0.6, which='major')    ax.grid(True, alpha=0.3, linestyle=':', linewidth=0.5, which='minor')    ax.set_xticks(x_ticks)    ax.tick_params(axis='both', labelsize=10)# --- Амплитудный спектр (график 6) ---# 📖 Связь со статьёй: Формула (1.10).ax_amp = fig.add_subplot(7, 1, 6)markerline, stemlines, baseline = ax_amp.stem(freq_nz, amp_nz, linefmt='C0-', markerfmt='o', basefmt='none')plt.setp(markerline, markersize=7.0, markeredgewidth=1)plt.setp(stemlines, linewidth=2, alpha=0.8)ax_amp.set_title("6. Амплитудный спектр", fontweight='bold', fontsize=14)ax_amp.set_ylabel('Амплитуда', fontsize=12)ax_amp.set_xlabel('Частота (Гц)', fontsize=12)if len(freq_nz) > 0: ax_amp.set_xlim(-0.5, max(freq_nz) + 0.5)ax_amp.set_xticks(np.arange(0, N//2 + 1, 5))ax_amp.grid(True, alpha=0.5, linestyle='-', linewidth=0.6, which='major')ax_amp.grid(True, alpha=0.3, linestyle=':', linewidth=0.5, which='minor')ax_amp.axhline(0, color='gray', linewidth=0.8, alpha=0.6)ax_amp.tick_params(axis='both', labelsize=10)# --- Фазовый спектр (график 7) ---# 📖 Связь со статьёй: Формула (1.11). Фаза синуса вычисляется относительно косинуса.ax_phase = fig.add_subplot(7, 1, 7)mask = ~np.isnan(phases)freq_phase = freq_indices[mask]phase_vals = phases[mask]markerline, stemlines, baseline = ax_phase.stem(freq_phase, phase_vals, linefmt='C1-', markerfmt='o', basefmt='none')plt.setp(markerline, markersize=7.0, markeredgewidth=1)plt.setp(stemlines, linewidth=2, alpha=0.8)ax_phase.set_title("7. Фазовый спектр", fontweight='bold', fontsize=14)ax_phase.set_ylabel('Фаза (градусы)', fontsize=12)ax_phase.set_xlabel('Частота (Гц)', fontsize=12)if len(freq_phase) > 0: ax_phase.set_xlim(-0.5, max(freq_phase) + 0.5)ax_phase.set_ylim(-200, 200)ax_phase.set_xticks(np.arange(0, N//2 + 1, 5))ax_phase.grid(True, alpha=0.5, linestyle='-', linewidth=0.6, which='major')ax_phase.grid(True, alpha=0.3, linestyle=':', linewidth=0.5, which='minor')ax_phase.axhline(0, color='gray', linewidth=0.8, alpha=0.6)ax_phase.tick_params(axis='both', labelsize=10)# ============================================================================# ДОБАВЛЕНА ТАБЛИЦА С ДАННЫМИ АМПЛИТУД И ФАЗ (ГРАФИКИ 6 и 7)# ============================================================================# Таблица размещена ниже, чтобы не перекрывать фазовый спектрax_table = fig.add_axes([0.1, 0.02, 0.8, 0.10])  # [left, bottom, width, height]ax_table.axis('off')# Подготовка данных: частота, амплитуда, фаза (для всех значимых частот)table_data = [['Частота (Гц)', 'Амплитуда', 'Фаза (°)']]for f, a, p in zip(freq_nz, amp_nz, phases[nonzero]):    phase_str = f"{p:.1f}" if not np.isnan(p) else "N/A"    table_data.append([f"{f}", f"{a:.3f}", phase_str])table = ax_table.table(cellText=table_data, loc='center', cellLoc='center')table.auto_set_font_size(False)table.set_fontsize(11)                     # Увеличен шрифт# Делаем весь текст жирнымfor (row, col), cell in table.get_celld().items():    cell.set_text_props(fontweight='bold')table.scale(1, 1.5)ax_table.set_title("Результат ДПФ", fontsize=11, fontweight='bold')# ============================================================================# 6: СОХРАНЕНИЕ# ============================================================================filepath = os.path.join(IMAGES_FOLDER, "orthogonality_7.png")fig.savefig(filepath, dpi=DPI, bbox_inches='tight', facecolor='white')print(f"\n✅ Изображение сохранено: {filepath}")# ============================================================================# 7: ДИАГНОСТИКА# ============================================================================# 📖 Связь со статьёй: Пример 7. Финальная проверка восстановления A и φ.print("\n" + "="*50)print("ПРОВЕРКА КОРРЕКТНОСТИ АМПЛИТУД И ФАЗ")print("="*50)for idx, p in enumerate(SIGNAL_PARAMS):    k = p["freq"]    if k < len(amplitudes):        A_calc = amplitudes[k]        phi_calc = phases[k] if not np.isnan(phases[k]) else "N/A"        print(f"Сигнал {idx+1}: Задано A={p['amp']:.1f}, φ={p['phase']:>3}° | "              f"Найдено A={A_calc:.3f}, φ={phi_calc}")    else:        print(f"Частота {k} выходит за пределы спектра (k ≤ {N//2})")print("="*50)plt.show()

На этом первая часть статьи завершается. Однако остаётся ещё много важных тем, связанных с преобразованием Фурье, которые заслуживают отдельного внимания: теорема Котельникова, утечка спектра, свёртка, быстрое преобразование Фурье (БПФ), обратное ДПФ, оконные функции и многие другие.

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

ссылка на оригинал статьи https://habr.com/ru/articles/1045118/