Если у вас возникли сложности с курсовой, контрольной, дипломной, рефератом, отчетом по практике, научно-исследовательской и любой другой работой - мы готовы помочь.
Предоплата всего
Подписываем
Глава 8 - Дискретное преобразование Фурье 22
Дискретное преобразование Фурье
Анализ Фурье есть семейство математических приемов, которые базируются на разложение сигналов на синусоиды. Дискретное преобразование Фурье (ДПФ) член семейства, который использует цифровые сигналы. Это первая из четырех глав, посвященных реальному ДПФ, версии дискретного преобразования Фурье, использующему реальные (вещественные) числа для представления входных и выходных сигналов. Комплексное ДПФ, более совершенный способ, использует комплексные числа, которые будут рассмотрены в главе 31. В этой главе мы будем рассматривать математические выражения и алгоритмы разложения Фурье, которые являются сердцем ДПФ.
Семейство преобразований Фурье
Фурье анализ назван в честь Джона Баптиста Джозефа Фурье (1768-1830), французского математика и физика. Хотя многие внесли вклад в эту область, Фурье является героем благодаря его математическим открытиям и практической интуиции в технике. Фурье интересовался распространением тепла и представил статью в 1807 году во Французскую Академию об использовании синусоид для описания температурного распределения. Статья содержала спорное утверждение, что любой непрерывный периодический сигнал может быть представлен как сумма соответственно выбранных синусоидальных волн. Два знаменитых математика просмотрели эту статью Джозеф Луис Лагранж (1736-1813) и Пьеро Симон де Лаплас (1749-1827).
В то время как Лаплас и другие вотировали публикацию статьи, Лагранж резко возражал. Около 50 лет Лагранж настаивал на том, что такой подход не может быть использован для представления сигналов с углами, т.е. изломанных подъемов, таких как квадратные волны. Французская Академия покорилась авторитету Лагранжа и не опубликовала работу Фурье. И только после смерти Лагранжа, статья была, в конце концов, опубликована через 15 лет. Фурье занимался и другими вещами: политикой, экспедицией в Египет с Наполеоном и пытался избежать гильотины после Французской революции (буквально!).
Кто был прав? Ответ двойственен. Лагранж был прав в своем утверждении, что сумма синусоид не сможет сформировать сигнал с углом. Однако она может подойти к нему так близко, что разница между суммой синусоид и угловым сигналом будет иметь нулевую энергию. В этом смысле Фурье был прав, хотя наука 18-го столетия мало знала о концепции энергии. Этот феномен теперь получил наименование эффект Гиббса (Gibbs Effect), который мы обсудим в Главе 11.
Рисунок 8-1 иллюстрирует, как сигнал может быть разложен на синусоидальные и косинусоидальные волны. Рисунок (а) показывает пример сигнала длиной в 16 точек с номерами отсчетов от 0 до 15. Рисунок (b) показывает разложение Фурье этого сигнала на 9 косинусов и 9 синусов каждый с разной частотой и амплитудой. Хотя это не очень очевидно, но эти 18 синусоид при сложении дают исходный сигнал (а). Необходимо заметить, что возражение сделанное Лагранжем применимо только к непрерывным сигналам. Для дискретных сигналов это разложение математически точно. Нет разницы между сигналом (а) и суммой сигналов (b) также как нет разницы между 7 и 3+4.
РИСУНОК 8-1а
Почему используются синусоиды, а не квадратные или треугольные волны? Вспомните, имеется бесконечное число способов разложения сигнала. Цель разложения заключается в том, чтобы получить что-то более легкое для работы, чем оригинальный сигнал. Например, импульсное разложение позволяет сигнал исследовать по каждой точке во времени, используя мощный аппарат свертки. Сигналы в виде синуса и косинуса проще, чем оригинальный сигнал потому, что они имеют свойства, которых нет у него - синусоидальную точность воспроизведения (sinusoidal fidelity). Как обсуждалось в главе 5, синусоида на входе системы гарантированно воспроизводит синусоиду на выходе. Только амплитуда и фаза сигнала может измениться, частота и форма волны должны сохраниться точно такими же. Синусоида единственная форма сигнала, которая обладает этим полезным свойством. Хотя разложение на квадраты и треугольники возможно, но нет причин для их использования.
Основной тип преобразования Фурье может быть представлен 4-мя категориями по 4-м типам сигналов, которые могут встретиться. Сигнал может быть непрерывным или дискретным и может быть периодическим или апериодическим. Комбинации этих 2-х характеристик дает 4 категории, описанные ниже и показанные на рисунке 8-2.
РИСУНОК 8-1b
Пример разложения Фурье. 16-точечный сигнал (с противоположной страницы) разлагается на 9 косинусных и 9 синусных волн. Частота каждой синусоиды фиксирована, только амплитуды изменяются в зависимости от формы сигнала.
Апериодический, непрерывный сигнал. Сюда включаются, например, экспоненциально затухающий сигнал и Гауссова кривая. Эти сигналы простираются от минус бесконечности до плюс бесконечности без какого-либо периодического повторения. Преобразование Фурье для этих сигналов называется просто преобразование Фурье.
Периодический непрерывный сигнал. Примеры – синусоидальные волны, квадратные волны и любой вид сигнала, который повторяется с регулярным рисунком от минус бесконечности до плюс бесконечности. Эта версия Фурье преобразования называется ряды Фурье (Fourier Series).
Апериодический дискретный сигнал. Эти сигналы определяются дискретными точками от минуса бесконечности до плюс бесконечности и не повторяются каким-либо периодическим образом. Этот тип преобразования Фурье называется дискретным временным преобразованием Фурье (Discrete Time Fourier Transform).
Периодический дискретный сигнал. Это дискретный сигнал, который повторяется каким-либо периодическим образом от минус бесконечности до плюс бесконечности. Этот класс преобразования Фурье иногда называется дискретные ряды Фурье, но чаще называется дискретное преобразование Фурье (Discrete Fourier Transform).
Вы можете подумать, что наименования этих типов преобразования Фурье случайны и плохо продуманы. Вы правы. Имена даны случайно более 200 лет назад. Здесь ничего вы не можете сделать, но только запомнить и пользоваться ими.
Все сигналы этих 4-х классов простирают от минуса бесконечности до плюс бесконечности. Ничего себе, скажете вы. Что, если вы имеете в вашем компьютере только конечное число отсчетов, скажем сигнал из 1024 точек. Ведь нет же версии для преобразования Фурье сигнала конечной длины. Нет. Синусы и косинусы определены, как простирающиеся от минуса бесконечности до плюс бесконечности. Вы не можете использовать группу сигналов бесконечной длины для синтезирования чего-либо конечного. Решение этой дилеммы заключается в том, чтобы конечные данные выглядели как бесконечные данные. Это осуществляется представлением, что сигнал имеет бесконечное число отсчетов справа и слева от точек с реальными данными. Если все эти «воображаемые» отсчеты имеют величину, равную нулю, сигнал выглядит дискретным и апериодическим и может применяться дискретное временное преобразование Фурье. Как альтернатива, воображаемые отсчеты могут быть представлены как копии реальных 1024 точек. В этом случае, сигнал выглядит дискретным и периодическим с периодом 1024 отсчета. Это позволяет использовать дискретное преобразование Фурье.
Оказывается, для синтеза апериодического сигнала требуется бесконечное число синусоид. Это делает невозможным вычисление дискретного временного преобразования Фурье в компьютерных алгоритмах. Как исключение, единственный тип преобразования Фурье, который может быть использован в ЦОС, есть ДПФ. Другими словами, цифровой компьютер может работать только с информацией, которая дискретна и конечна. Когда вы боретесь с теоретическими вопросами, пытаясь преодолеть домашней работой проблемы, и обдумываете математические выкладки, вы можете найти для себя полезным первые 3 члена из семейства преобразований Фурье. Но, когда вы сядете за свой компьютер, вы сможете использовать только ДПФ. Мы будем коротко рассматривать другие преобразования Фурье в последующих главах, но теперь сконцентрируйтесь на понимании дискретного преобразования Фурье.
Посмотрите пример разложения ДПФ на рисунке 8-1. 16-ти точечный сигнал раскладывается на 18 синусоид, каждая состоит из 16 точек. В более формальных терминах, 16-ти точечный сигнал, показанный в (а), может рассматриваться, как сигнал с бесконечным периодом. Подобно этому, каждая из 18-ти синусоид, показанная в (b), представляет 16-ти точечный сегмент синусоиды бесконечной длины. Имеет ли значение, будем ли мы рассматривать это разложение, как 16-ти точечный сигнал, синтезированный из 16-ти точечных синусоид, или как сигнал с бесконечно длинным периодом, синтезированный из бесконечно длинных синусоид? Обычно нет, но иногда да. В следующих главах мы встретимся со свойствами ДПФ, которые кажутся сбивающими с толку, если сигнал рассматривать как бесконечный, но становятся очевидными, когда предполагается периодическая природа сигнала. Главное понять, что эта периодичность вызвана порядком использования математического инструмента, т.е. ДПФ. И не играет роли, откуда берется сигнал или как он получается.
РИСУНОК 8-2
Иллюстрация четырех преобразований Фурье. Сигнал может быть непрерывным и дискретным; периодическим и апериодическим. Всего четыре комбинации, каждая имеет собственную версию преобразования Фурье. Имена неудачны, просто запомните их.
Каждое из четырех преобразований Фурье может быть подразделено на реальную и комплексную версии. Реальная версия самая простая, использует обычные числа и алгебру для синтеза и разложения. Например, рисунок 8-1 есть пример реального ДПФ. Комплексная версия безмерно более сложная, требует использования комплексных чисел. Это числа типа 3 + 4j, где j равна корню из минус единицы (инженеры используют j, математики i). Математика комплексных чисел может быстро ошеломлять даже тех, кто специализируется в ЦОС. В действительности, главная цель этой книги представить основы ЦОС без использования комплексных чисел, делая материал понятным широкому кругу научных работников и инженеров. Комплексное преобразование Фурье есть царство тех, кто специализируется в ЦОС и охотно погружается по самую шею в болото математики. Если у вас есть к этому наклонность, главы 30 – 33 перенесут вас туда.
Математический термин преобразование широко используется в методах ЦОС таких как: преобразование Фурье, преобразование Лапласа, Z-преобразование, преобразование Гильберта, дискретно-косинусное преобразование и т.д. Так что ж такое преобразование? Отвечая на этот вопрос, помните, что это функция. Функция есть алгоритм или процедура, которая изменяет одну величину в другую. Например, y = 2x + 1 есть функция. Вы берете некоторую величину х, вставляете ее в формулу, и на выходе получаете у. Функция так же может изменять несколько величин в одну величину, как: y = 2a + 3b + 4c, где a, b, c преобразуются в у.
Преобразование есть прямое расширение этого, позволяющее как входной величине, так и выходной иметь многочисленные (multiple) значения. Предположим, вы имеет сигнал, состоящий из 100 отсчетов. Если вы придумали некоторую формулу, алгоритм или процедуру для изменения этих 100 отсчетов в другие 100 отсчетов, вы имеете собственное преобразование. Если вы думаете, что оно достаточно полезно, то вы имеете полное право присвоить ему свое имя и изложить своим коллегам его достоинства. (Особенно, если вы знаменитый французский математик 18-го столетия). Преобразования не лимитируются каким-либо специфическим типом или числом данных. Например, вы можете иметь 100 отсчетов дискретных данных для входного сигнала и 200 отсчетов дискретных данных для выходного. Подобно этому, вы можете иметь непрерывный сигнал на входе и непрерывный сигнал на выходе. Смешенные сигналы также допускаются, дискретный на входе и непрерывный на выходе, и наоборот. Короче, преобразование есть любой стационарный процесс, который изменяет один кусок данных в другой кусок данных. Давайте посмотрим, как это осуществляется в разделе с заголовком: дискретное преобразование Фурье.
Система обозначений и формат реального ДПФ
Как показано на рисунке 8-3, дискретное преобразование Фурье изменяет N точек входного сигнала в два выходных сигнала из N/2 + 1 точек. Входной сигнал это сигнал, который подвергается разложению, а два выходных сигнала содержат амплитуды синусов и косинусов (способ определения масштаба мы обсудим немного позже). Входной сигнал предназначен находиться во временной области. Это вызвано тем, что большинство типичных сигналов, встречающихся в ЦОС, состоят из отсчетов, полученных через регулярные интервалы времени. Конечно, любой вид отсчетов может быть использован в ЦОС независимо от того, как он получен. Когда вы в анализе Фурье встречаете термин «временная область», это может действительно означать, что отсчеты получены во времени, или это может означать любой дискретный сигнал, подвергающийся разложению. Термин частотная область используется для описания амплитуд синусов и косинусов (включая специальное масштабирование, которое мы обещаем объяснить).
РИСУНОК 8-3
Терминология ДПФ. Во временной области сигнал x[ ] состоит из N точек с номерами от 0 до N-1. В частотной области ДПФ создает два сигнала: реальная часть, ReX[ ], и мнимая часть, ImX[ ]. Каждый из этих сигналов имеет N/2 + 1 точек с номерами от 0 до N/2. Прямое ДПФ преобразует временную область в частотную, инверсное ДПФ преобразует частотную область во временную. (Здесь описывается реальное ДПФ; комплексное ДПФ, описанное в главе 31, изменяет N комплексных точек другой набор N комплексных точек).
Частотная область содержит точно такую же информацию, как и временная область, но в другой форме. Если вам известна одна область, мы можете вычислить другую. Если получен сигнал во временной области, процесс вычисления сигнала в частотной области называется разложение, анализ, прямое ДПФ или просто ДПФ. Если известна частотная область, вычисление сигнала во временной области называется синтезом или инверсным ДПФ. Как синтез, так и анализ могут быть представлены в виде формул и в виде компьютерных алгоритмов.
Число отсчетов во временной области обычно обозначаются переменой N. В то время как N может быть любой положительной величиной, выбирается оно обычно кратным двум, т.е. 128, 256, 512, 1024 и т.д. Имеется две причины для этого. Первая, при хранении цифровых данных используется бинарная адресация, что делает естественную длину сигнала кратной двум. Вторая, большинство эффективных алгоритмов для вычисления ДПФ, как, например, быстрое преобразование Фурье (БПФ), используют величину N кратную двойке. Обычно N выбирают между 32 и 4096. В большинстве случаев отсчеты нумеруют от 0 до N-1, а не от 1 до N.
Стандартные обозначения в ЦОС используют строчные буквы для обозначения сигналов во временной области, такие как x[ ], y[ ], z[ ]. Соответственные заглавные буквы используются для обозначения сигналов в частотной области X[ ], Y[ ], Z[ ]. Для иллюстрации, возьмем N точек сигнала во временной области, который содержится в x[ ]. В частотной области этот сигнал будет называться X[ ] и состоять из двух частей, каждая из которых представляет массив из N/2 + 1 отсчетов. Эти части называются реальная часть X[ ], записывается как Re X[ ], и мнимая часть X[ ], записывается как Im X[ ]. Величины в Re Х[ ] есть амплитуды косинусов, величины в Im Y[ ] есть амплитуды синусов (пока не заботьтесь о масштабе).Так же как во временной области отсчеты идут от x[0] до x[N-1], в частотной области сигналы идут от Re X[0] до Re X[N/2] и от Im X[0] до Im X[N/2]. Изучите внимательно эти обозначения, они важны для понимания выражений в ЦОС. К сожалению, некоторые компьютерные языки не отличают строчные буквы от заглавных, что приводит к разным обозначениям у разных программистов. В программах этой книги используется обозначение XX[ ] для массива сигнала во временной области и REX[ ], IMX[ ] для массивов в частотной области.
Наименования реальная часть и мнимая часть берут начало от комплексного ДПФ, где они используются для обозначения отличия между реальными и мнимыми числами. Ничего такого сложного не требуется для реального ДПФ. Вплоть до главы 31 просто считайте, что «реальная часть» обозначает амплитуду косинусов, а «мнимая часть» - амплитуду синусов. Не позволяйте этим наименованиям сбивать вас с толку, здесь везде используются обычные числа.
Также не смущайтесь термином длина по отношению к сигналу в частотной области. В литературе по ЦОС часто можно увидеть предложение, которое выглядит так: «ДПФ изменяет N точек сигнала во временной области в N точек сигнала в частотной области». Это относится к комплексному ДПФ, где каждая «точка» есть комплексное число, состоящее из реальной и мнимой части. Теперь сфокусируемся на реальном ДПФ, скоро подойдет трудная математика.
Независимая переменная в частотной области
Рисунок 8-4 показывает пример ДПФ для N = 128. Сигнал временной области состоит из массива x[0].....x[127]. Сигнал частотной области содержит два массива ReX[0]...ReX[64] и ImX[0]...ImX[64]. Заметим, что 128 точек временной области соответствуют 65 точкам в каждом массиве частотной области, с номерами частот от 0 до 64. То есть, N точек временной области соответствуют N/2 + 1 точек частотной области (не N/2 точек). Забывание этой дополнительной точки есть обычная ошибка при программировании ДПФ.
Горизонтальная ось в частотной области может быть обозначена четырьмя различными способами, все они общеприняты в ЦОС. При первом способе горизонтальная ось размечена номерами отсчетов от 0 до 64 (N/2). В этом случае индексы в частотной области есть целые числа, например, ReX[k] и ImX[k], где k меняется от 0 до N/2 с шагом единица. Программисты любят этот способ, потому что он соответствует коду, который они используют для доступа к массиву. Это обозначение используется на рисунке 8-4b.
РИСУНОК 8-4
Пример ДПФ. ДПФ преобразует сигнал временной области, x[ ], в сигналы частотной области, ReX[ ] и ImX[ ]. Горизонтальная ось в частотной области может быть отмечена тремя способами: (1) как массив индексов от 0 до N/2; (2) как часть частоты дискретизации от 0 до 0,5; (3) как натуральная частота от 0 до . В этом примере (b) использует первый метод, (с) – второй.
При втором методе, который показан в (с), горизонтальная ось размечена в долях частоты дискретизации (fraction of the sampling rate). Это означает, что величина вдоль горизонтальной оси всегда изменяется от 0 до 0,5; так как дискретные данные могут содержать частоты только от постоянного тока до 1\2 частоты дискретизации. Индекс, используемый для этого обозначения, есть f. Реальная и мнимая части записываются как ReX[f] и ImX[f], где f принимает значения N/2 + 1 величин, равномерно расположенных между 0 и 0,5. Преобразование первого обозначения, k, во второе, f, заключается в делении горизонтальной оси на N. То есть, f = k/N. Большинство графиков в этой книге используют этот второй метод, фокусируя внимание на том, что дискретный сигнал может содержать частоты только от 0 до 0,5 частоты дискретизации.
Третий метод – это просто второй, за исключением того, что горизонтальная ось умножена на 2. Индекс, используемый для разметки - . При этом обозначении реальная и мнимая части записываются как: ReX[] и ImX[], где принимает значения N/2 + 1 величин, равномерно расположенных между 0 и . Параметр называется натуральная частота и имеет размерность в радианах. Основная идея заключается в том, что в цикле 2 радиан. Математики любят этот способ обозначений, потому что он позволяет короче записать формулы. Например, рассмотрим запись косинуса для каждого из этих обозначений. Используя k: c[n] = cos(2kn/N), используя f: c[n] = cos(2fn), используя : c[n] = cos(n).
Четвертый способ заключается в разметке горизонтальной оси в величинах аналогичных частотам, используемых в практических приложениях. Например, если в системе применяется частота дискретизации 10 кГц (т.е. 10 000 отсчетов в секунду), то на графике частота будет меняться от 0 до 5 кГц. Этот способ хорош для представления частотных данных в обозначениях, используемых на практике. Недостаток заключается в том, что он привязан к частному уровню дискретизации, и поэтому не применим для развития общих алгоритмов ЦОС, таких как создание цифровых фильтров.
Все эти четыре способа обозначения используются в ЦОС, и вам необходимо хорошо разобраться с преобразованием их друг в друга. Это относится как к графикам, так и к математическим формулам. Чтобы определить какой способ обозначения используется, посмотрите на независимую величину и уровень ее величины. Вы можете увидеть одно из четырех обозначений: k (или какой-либо другой целочисленный индекс), меняющийся от 0 до N/2; f, меняющееся от 0 до 0,5; или частоту в герцах, меняющуюся от постоянного тока (DC) до одной второй частоты дискретизации.
Базовые функции ДПФ
Синусы и косинусы, используемые в ДПФ, обычно называются базовыми функциями. Другими словами, выход ДПФ есть набор чисел, которые представляют амплитуды. Базовые функции есть набор синусов и косинусов с единичной амплитудой. Если вы определили каждую амплитуду (в частотной области) соответствующего синуса или косинуса (базовые функции), то в результате вы получите масштабированные синусы и косинусы, сумма которых сможет сформировать сигнал во временной области.
Базовые функции ДПФ определяются следующими выражениями:
ck[i] = cos(2ki/N)
sk[i] = sin(2ki/N)
где: ck[i] – косинусная волна с амплитудой, содержащейся в ReX[k], и sk[i] – синусная волна с амплитудой, содержащейся в ImX[k]. Например, рисунок 8-5 показывает некоторые из 17 синусов и 17 косинусов, используемых в N = 32-точечном ДПФ. Поскольку эти синусоиды складываются для формирования входного сигнала, они должны быть той же самой длины, как и входной сигнал. В этом случае каждая синусоида имеет 32 точки с номерами от 0 до 31. Параметр k определяет частоту каждой синусоиды. Практически, c1[ ] это косинусная волна, которая совершает один полный цикл колебаний на N точках, c5[ ] – косинусная волна с пятью полными колебаниями на N точках и т.д. Это важная концепция в понимании базовых функций; частотный параметр, k, равен числу полных колебаний (циклов), приходящихся на N точек сигнала.
РИСУНОК 8-5
Базовые функции ДПФ. 32-точечное ДПФ имеет 17 дискретных косинусных волн и 17 дискретных синусных волн в качестве базовых функций. Восемь из них показаны на рисунке. Эти сигналы дискретны, непрерывная линия на рисунках только помогает глазу отслеживать формы волны.
Давайте детально рассмотрим некоторые из этих базовых функций. Рисунок (а) показывает косинусную волну c0[ ]. Эта косинусная волна имеет нулевую частоту, т.е. это постоянная единичной величины. Это означает, что ReX[0] содержит величину, усредненную по всем точкам сигнала во временной области. В электронике можно было бы сказать, что ReX[0] содержит постоянный ток (DC offset). Синусная волна нулевой частоты s0[ ], показанная в (b), состоит из одних нулей. Поскольку она не воздействует на формирование сигнала во временной области, то величина ImX[0] не относится к делу и всегда устанавливается в 0. Подробнее об этом позже.
Рисунки (с) и (d) показывают с2 [ ]и s2[ ] синусоиды, которые совершают 2 полных колебания на отрезке из N точек. Им соответствуют ReX[2] и ImX[2]. Подобным образом, (e) и (f) показывают c10[ ] и s10[ ] синусоиды, которые совершают 10 полных колебаний на N точках. Этим синусоидам соответствуют амплитуды ReX[10] и ImX[10]. Проблема заключается в том, что отсчеты в (e) и (f) не отражают наглядно характер синусных и косинусных волн. Если бы не было непрерывной кривой на графике, вам было бы очень трудно определить вид сигнала. Это может доставить вам маленькое неудобство, но не беспокойтесь об этом. С математической точки зрения эти отсчеты формируют дискретные синусоиды, хотя ваш глаз и не отслеживает этот рисунок.
Самые высокие частоты в базовых функциях показаны в (g) и (h). Это cN/2 и sN/2, или для нашего примера c16[ ] и s16[ ]. Дискретная косинусная волна чередуется между 1 и -1, что может быть интерпретировать как отсчеты по пикам непрерывной синусоиды. В противоположность этому дискретная синусная волна состоит из всех нулей, в результате попадания отсчетов на нулевые пересечения. Это делает величину ImX[N/2] такой же, как и ImX[0], всегда равной нулю и не влияющей на синтез сигнала во временной области.
Здесь есть загадка: Если имеется N отсчетов на входе ДПФ и N+2 отсчетов на выходе, то откуда берется дополнительная информация? Ответ: два выходных отсчета не содержат информацию, позволяя другим N быть полностью независимыми. Как Вы, возможно, догадались, точки, которые не несут информацию, есть ImX[0] и ImX[N/2] отсчеты, которые всегда имеют нулевую величину.
Синтез, вычисление инверсного (обратного) ДПФ
Обобщая все, что будет сказано дальше, мы можем написать уравнение синтеза:
Иными словами, любой сигнал из N точек x[i] может быть создан суммированием N/2+1 косинусных волн и N/2 + 1 синусных волн. Амплитуды косинусных и синусных волн содержатся в массивах Re[k] и Im[k], соответственно. Уравнение синтеза перемножает эти амплитуды на базовые функции, создавая набор масштабированных синусных и косинусных волн. Сложение этих масштабированных синусных и косинусных волн формирует сигнал временной области x[i].
В выражении 8-2 массивы называются Im[k] и Re[k], вместо Im X [k] и Re X [k]. Это вызвано тем, что амплитуды, необходимые для синтеза, (названные в этом обсуждении Im[k] и Re[k]) слегка отличаются от сигналов частотной области (обозначены Im X [k] и Re X [k]). Это влияние фактора масштабирования, о котором мы упоминали ранее. Хотя это преобразование всего лишь нормализация, оно часто приводит к ошибкам при программировании. Будьте внимательны! В формульном выражении преобразование выглядит следующим образом:
За исключением двух специальных случаев:
Предположим, вам представлена частотная область, и предлагается синтезировать сигнал временной области. Для начала вы должны найти амплитуды синусов и косинусов. Другими словами, взяв Im X [k] и Re X [k] вы должны найти Im[k] и Re[k]. Формула 8-3 показывает, как это сделать математически. В компьютерной программе необходимо выполнить три действия. Первое, разделить все величины в частотной области на N/2. Второе, изменить знак перед всеми мнимыми величинами. Третье, разделить первую и последнюю величину в реальной части, Re X [0] и Re X [N/2], на 2. Это обеспечит величину амплитуд, необходимую для синтеза по формуле 8-2. Взятые вместе формулы 8-2 и 8-3 определяют инверсное ДПФ.
Целиком инверсное ДПФ в компьютерной программе показано в таблице 8-1. Имеется два способа программирования синтеза (по формуле 8-2), и оба они показаны. В первом способе, каждая из масштабированных синусоид вычисляется один раз и суммируется в аккумуляторе, процесс заканчивается получением сигнала временной области. При втором способе, каждый отсчет сигнала временной области вычисляется один раз, как сумма всех соответствующих отсчетов косинусных и синусных волн. Оба метода дают тот же самый результат. Разница между этими двумя программами очень не большая: внутренняя и внешняя петли циклов меняются местами.
100 ИНВЕРСНОЕ ДИСКРЕТНОЕ ПРЕОРАЗОВАНИЕ ФУРЬЕ
110 'Сигнал временной области содержится в XX[ ], он вычисляется из сигнала
120 'частотной области, который содержится в REX[ ] и IMX[ ].
130 '
140 DIM XX[511] 'XX[ ] содержится сигнал временной области
150 DIM REX[256] 'REX[ ] сдержится реальная часть сигнала частотной области
160 DIM IMX[256] 'IMX[ ] содержится мнимая часть сигнала частотной области
170 '
180 PI = 3.14159265 Установка константы, PI
190 N% = 512 'N% число точек в XX[ ]
200 '
210 GOSUB XXXX 'Воображаемая подпрограммы для загрузки REX[ ] и IMX[ ]
220 '
230
240 ' 'Нахождение амплитуд косинусных и синусных волн по форм. 8.3
250 FOR K% = 0 TO 256
260 REX[K%] = REX[K%] / (N%/2)
270 IMX[K%] = -IMX[K%] / (N%/2)
280 NEXT K%
290 '
300 REX[0] = REX[0] / N%
310 REX[256] = REX[256] / N%
320 '
330 '
340 FOR I% = 0 TO 511 'обнуление XX[ ] перед аккумуляцией
350 XX[I%] = 0
360 NEXT I%
370 '
380 ' Формула 8-2 СИНТЕЗ #1. Цикл через каждую частоту
390 ' по целой длине косинусных и синусных волн
400 ' и хранение их в аккумуляторе XX[ ]
410 '
420 FOR K% = 0 TO 256 'K% цикл через каждый отсчет в REX[ ] и IMX[ ]
430 FOR I% = 0 TO 511 'I% цикл через каждый отсчет в XX[ ]
440 '
450 XX[I%] = XX[I%] + REX[K%] * COS(2*PI*K%*I%/N%)
460 XX[I%] = XX[I%] + IMX[K%] * SIN(2*PI*K%*I%/N%)
470 '
480 NEXT I%
490 NEXT K%
500 '
510 END
Альтернативный код в линиях от 380 до 510
380 ' Формула 8.2 СИНТЕЗ #2. Цикл через каждый отсчет
390 ' временной области, суммирование соответствующих
400 ' отсчетов из каждой косинусной и синусной волны
410 '
420 FOR I% = 0 TO 511 'I% цикл через каждый отсчет в XX[ ]
430 FOR K% = 0 TO 256 'K% цикл через каждый отсчет в REX[ ] и IMX[ ]
440 '
450 XX[I%] = XX[I%] + REX[K%] * COS(2*PI*K%*I%/N%)
460 XX[I%] = XX[I%] + IMX[K%] * SIN(2*PI*K%*I%/N%)
470 '
480 NEXT K%
490 NEXT I%
500 '
510 END
ТАБЛИЦА 8-1
РИСУНОК 8-6
Пример инверсного ДПФ. Рисунок (а) показывает пример сигнала временной области, импульс амплитудой 32. Рисунок (b) показывает реальную часть частотной области этого сигнала, постоянную величину равную 32. Мнимая часть не показана, так как она вся состоит из нулей. Рисунок (с) показывает амплитуду косинусных волн, необходимых для реконструкции (а) по формуле 8-2. Величина (с) находится из (b) по формуле 8-3.
Рисунок 8-6 иллюстрирует инверсное ДПФ, а так же разницу между частотной областью и амплитудами, необходимыми для синтеза. Рисунок 8-6а показывает сигнал, который синтезируется, - импульс с амплитудой равной 32 в нулевом отсчете. Рисунок 8-6b показывает частотную область, соответствующую этому сигналу. Мнимая часть (не показана) состоит из одних нулей. Как обсуждается в следующей главе, это важная пара ДПФ: импульсу во временной области соответствует постоянная величина в частотной области. Теперь важный момент: (b) есть ДПФ от (а) и (а) есть инверсное ДПФ от (b).
Формулы 8-3 используются для преобразования сигнала частотной области, (b), в амплитуды косинусных волн, (с). Как показано, амплитуды всех косинусных волн имеет амплитуду равную двум, за исключением отсчетов 0 и 16, имеющих единичную величину амплитуды. Амплитуды синусных волн не показаны в этом примере, потому что они все имеют нулевую амплитуду, и поэтому роли не играют. Затем используется формула синтеза, 8-2, для преобразования амплитуды косинусных волн, (b), в сигнал временной области, (а).
Здесь описывается, чем частотная область отличается от синусоидальных амплитуд, но не объясняется почему. Отличие возникает из-за того, что частотная область определена как спектральная плотность. Рисунок 8-7 поясняет это. На рисунке показан пример реальной части частотной области для 32 точечного сигнала. Как вы и ожидали, отсчеты идут от 0 до 16, представляя 17 частот равномерно расположенных между 0 и ½ частоты дискретизации. Спектральная плотность описывает как много сигнала (амплитуд) представлено на единицу полосы частот. Для преобразования синусоидальных амплитуд в спектральную плотность разделите каждую амплитуду на ширину полосы частот, приходящуюся на каждую амплитуду. Это вызывает следующий вопрос: Как мы определим ширину полосы для каждой дискретной частоты в частотной области?
Как показано на рисунке, ширина полосы может быть определена путем нанесения разделяющих линий между отсчетами. Например, отсчет номер 5 оказывается в полосе между 4,5 и 5,5; отсчет номер 6 оказывается в полосе между 5,5 и 6,5 и т.д. Представляя эту полосу как часть общей полосы (т.е. N/2), мы получим ширину полосы для каждого отсчета равной 2/N. Исключение составят отсчеты на концах, каждый будет иметь половинку от этой полосы, 1/N. Это согласуется с коэффициентом масштабирования, 2/N, меду синусоидальными амплитудами и частотной областью, и с добавочным коэффициентом, два, необходимым для первого и последнего отсчета.
Почему необходим знак минус перед мнимой частью? Только лишь для согласования реального ДПФ с его старшим братом, комплексным ДПФ. В главе 29 мы покажем, что это необходимо для работы комплексного ДПФ с математической точки зрения. Когда вопрос касается только реального ДПФ, многие авторы не учитывают этот минус. Случается, многие авторы даже не включают коэффициент масштабирования 2/N. Будьте готовы встретить эти ошибки в некоторых статьях. Они отмечены здесь по двум важным причинам. Большинство эффективных способов вычисления ДПФ заключается в использовании алгоритмов быстрого преобразования Фурье (БПФ), которые представлены в главе 12. БПФ создает частотную область согласно формулам 8-2 и 8-3. Если вы упустили этот фактор нормализации, ваша программа, содержащая БПФ, не будет работать так, как ожидалось.
РИСУНОК 8-7
Ширина отсчетов частотной области. Каждый отсчет в частотной области может рассматриваться как находящийся в полосе частот шириной 2/N по отношению к общей полосе частот. Исключение составляют первый и последний отсчеты, которые имеют половинную ширину, 1/N.
Анализ, вычисление ДПФ
ДПФ может быть вычислено тремя совершенно разными способами. Первый, проблема решается с помощью системы уравнений. Этот метод полезен для понимания ДПФ, но слишком непроизводителен для практического использования. Второй метод основан на идее из последней главы, корреляции. Он базируется на обнаружении известных волн в данном сигнале. Третий метод, называемый быстрое преобразование Фурье (БПФ), есть остроумный алгоритм, который заменяет N-точечное ДПФ на N одноточечных ДПФ. БПФ обычно в сотни раз быстрее, чем другие методы. Первые два метода обсуждаются здесь, БПФ – в главе 12. Важно помнить, что все три метода дают один и тот же результат. Который метод использовать? На практике, корреляционный метод предпочтителен, если ДПФ имеет меньше чем 32 точки, в других случаях надо использовать БПФ.
ДПФ с помощью системы уравнений
Это ДПФ вычисляется следующим способом. Вы имеете N величин из временной области, и вам предлагается вычислить N величин в частотной области (игнорируя две величины в частотной области, которые, как вы знаете, должны быть равны нулю). В соответствии с правилами алгебры, для нахождения N неизвестных вы должны составить N независимых линейных уравнений. Для этого, берутся первые отсчеты из каждой синусоиды и складываются вместе. Сумма должна быть равна первому отсчету сигнала временной области, так получается первое уравнение. Аналогичным образом, могут быть составлены уравнения для оставшихся точек сигнала временной области, в результате получаются необходимые N уравнений. Решение может быть найдено использованием основных методов для решения подобных уравнений, таких как подстановка Гаусса. К несчастью этот метод требует громадного числа вычислений, и, фактически, никогда не используется в ЦОС. Однако он важен по другой причине, он показывает, почему возможно разложение сигнала на синусоиды, как много синусоид требуется, и что базовые функции должны быть линейно независимые (более подробно об этом позже).
ДПФ с помощью корреляции
Давайте перейдем к лучшему способу, стандартному методу вычисления ДПФ. Покажем этот метод на примере. Предположим, мы пытаемся вычислить ДПФ для 64-точечного сигнала. Это означает, что нам необходимо вычислить 33 точки в реальной части и 33 точки в мнимой части частотной области. В этом примере мы покажем вычисление только одного отсчета, ImX[3], то есть амплитуду синусной волны, которая делает три полных колебания между точками 0 и 63. Все другие величины в частотной области вычисляются подобным образом.
Рисунок 8-8 иллюстрирует использование корреляции для вычисления ImX[3]. Рисунки (а) и (b) показывают пример двух сигналов временной области, названных x1[ ] и x2[ ], соответственно. Первый сигнал, x1[ ], состоит только из одной синусоиды, делающей три колебания меду точками 0 и 63. В противоположность этому, x2[ ] состоит из нескольких синусных и косинусных волн, ни одна из которых не делает три колебания между точками 0 и 63. Эти два сигнала иллюстрируют, что должен делать алгоритм для вычисления ImX[3]. Для x1[ ] алгоритм должен получить величину равную 32, амплитуду синусной волны представленной в сигнале (с учетом коэффициента масштабирования согласно формуле 8-3). Для x2[ ] должна получиться нулевая величина, свидетельствующая о том, что данная частная синусоида не присутствует в сигнале.
РИСУНОК 8-8
Два примера сигналов (а) и (b), анализируемых на наличие в них базовых функций, показанных в (с) и (d). Рисунки (е) и (f) показывают результат перемножения анализируемых сигналов на базовые функции. Рисунок (е) имеет среднее значение 0,5, что свидетельствует о наличии базовой функции с амплитудой 1 в анализируемом сигнале. (f) имеет нулевое среднее значение, что свидетельствует об отсутствии базовой функции в x2[ ].
Концепция корреляции была представлена в главе 7. Как вы помните, для обнаружения сигнала известной формы в исследуемом сигнале необходимо перемножить эти два сигнала и сложить все точки полученного сигнала. Число, полученное в результате этой процедуры, является мерой похожести этих двух сигналов. Рисунок 8-8 иллюстрирует этот подход. Рисунки (с) и (d) показывают сигнал, который мы сравниваем на похожесть, - синусоиду с тремя колебаниями между отсчетами 0 и 63. Рисунок (е) показывает результат перемножения (а) и (с). Рисунок (f) показывает результат перемножения (b) и (d). Сумма всех точек в (е) равна 32, в то время как сумма всех точек в (f) равна нулю. Таким образом, мы нашли желаемый алгоритм.
Другие отсчеты в частотной области вычисляются точно так же. Эта процедура формализуется в следующем уравнении анализа, математических формулах для нахождения сигналов частотной области по сигналам временной области:
Другими словами, каждый отсчет в частотной области находится перемножением сигнала временной области на рассматриваемые синусные и косинусные волны и суммированием результатов умножения. Если кто-нибудь спросит вас, что вы делаете, уверенно говорите: «Я ищу корреляцию входного сигнала с базовыми функциями». Таблица 8-2 показывает компьютерную программу, которая вычисляет ДПФ этим способом.
Уравнения анализа не требуют специального выражения для первой и последней точки в отличие от уравнения синтеза. Имеется, однако, знак минуса перед мнимой частью в формуле 8-4. Так же как и раньше, этот знак минус делает реальное ДПФ совместимым с комплексным ДПФ, и не всегда учитывается.
Для правильной работы этого корреляционного алгоритма, базовые функции должны обладать интересным свойством: каждая из них должна быть совершенно некоррелирована со всеми другими. Это означает, что если вы перемножите любые две базовые функции, то сумма полученных точек должна быть равна нулю. Базовые функции, обладающие этим свойством, называются ортогональными. Существует много ортогональных базовых функций: квадратные волны, треугольные волны, импульсные и т.д. Сигналы могут быть разложены на эти базовые ортогональные функции, так же как и на синусоиды. Это не означает, что это полезно, но только, что это возможно.
Как было показано в таблице 8-1, инверсное ДПФ имеет два способа выполнения в компьютерной программе. Разница связана с тем, что внутренняя и внешняя петли циклов меняются местами. Это не меняет результат вычислений, но меняет ваш взгляд на выполнение программы. Программа ДПФ в таблице 8-2 может быть изменена таким же образом, переменой мест внутренней и внешней петли циклов в линиях от 310 до 380. Так же как и раньше, результат вычислений не меняется, но меняется ваш взгляд на процесс вычисления. (Эти два разных взгляда на ДПФ и инверсное ДПФ могут быть описаны как «входной алгоритм» и «выходной алгоритм» аналогично свертке).
100 'ДИСКРЕТНОЕ ПРЕОБРАЗОВАНИЕ ФУРЬЕ (ДПФ)
110 'Сигнал частотной области храниться в REX[ ] и IMX[ ], вычисляется из
120 'сигнала временной области, который хранится в XX[ ].
130 '
140 DIM XX[511] 'XX[ ] хранится сигнал временной области
150 DIM REX[256] 'REX[ ] хранится реальная часть сигнала частотной области
160 DIM IMX[256] 'IMX[ ] хранится мнимая часть сигнала частотной области
170 '
180 PI = 3.14159265 'Установка константы, PI
190 N% = 512 'N% число точек в XX[ ]
200 '
210 GOSUB XXXX 'Воображаемая подпрограмма загрузки XX[ ]
220 '
230 '
240 FOR K% = 0 TO 256 'Обнуление REX[ ] & IMX[ ] для аккумуляции
250 REX[K%] = 0
260 IMX[K%] = 0
270 NEXT K%
280 '
290 ' 'Корреляция XX[ ] с косинусными и синусными волнами по формуле 8-4
300 '
310 FOR K% = 0 TO 256 'K% цикл через каждый отсчет в REX[ ] и IMX[ ]
320 FOR I% = 0 TO 511 'I% цикл через каждый отсчет в XX[ ]
330 '
340 REX[K%] = REX[K%] + XX[I%] * COS(2*PI*K%*I%/N%)
350 IMX[K%] = IMX[K%] - XX[I%] * SIN(2*PI*K%*I%/N%)
360 '
370 NEXT I%
380 NEXT K%
390 '
400 END
ТАБЛИЦА 8-2
Программа в таблице 8-2 написана таким образом, что показывает воздействие всех отсчетов временной области на отдельный отсчет в частотной области. То есть, программа вычисляет каждую величину частотной области в отдельности, а не как группу. Если внутренний цикл и внешний цикл поменять местами, то программный цикл через каждый отсчет во временной области будет вычислять вклад этой точки в частотную область. Целиком частотная область находится сложением этих вкладов от отдельных точек временной области. Отсюда возникает следующий вопрос: какова природа влияния отдельного отсчета во временной области на создание частотной области? Ответ заключается в интересном аспекте анализа Фурье, называемом двойственность (duality).
Двойственность
Уравнения синтеза и анализа (формулы 8-2 и 8-4) разительно похожи. Для перемещения из одной области в другую надо перемножить известные величины на базовые функции и сложить результаты. То, что ДПФ и инверсное ДПФ используют тот же самый математический подход, является замечательным фактом; учитывая итоговую разность вычислений, мы получаем две процедуры. Фактически, единственная разница между этими двумя формулами заключается в том, что временная область содержит один сигнал из N точек, а частотная область – два сигнала из N/2 + 1 точек. Как обсуждалось в последних главах, комплексное ДПФ выражает оба сигнала, как во временной области, так и в частотной в виде N точечного сигнала. Это делает эти две области совершенно симметричными, и формулы для перехода между ними фактически идентичными.
Эта симметрия между временной областью и частотной областью называется двойственность, и приводит к интересным свойствам. Например, единственная точка в частотной области соответствует синусоиде во временной области. Благодаря двойственности обратное утверждение также справедливо. Одна точка во временной области соответствует синусоиде в частотной области. Другой пример, свертка во временной области соответствует перемножению в частотной области. Справедливо и обратное утверждение: свертка в частотной области соответствует перемножению во временной области. Эти и другие соотношения двойственности обсуждаются более подробно в главах 10 и 11.
Полярная система координат
Как описано далее, частотная область есть группа амплитуд косинусных и синусных волн (с учетом масштабирования). Это называется прямоугольной системой координат. Альтернативным способом частотная область может быть представлена в полярной форме. В этой системе координат ReX[ ] и ImX[ ] преобразуются в два других массива, называемых амплитуда Х (magnitude), записывается в формулах как Mag X[ ], и фаза Х, записывается как Phase X[ ]. Амплитуда и фаза попарно вычисляются из реальной и мнимой частей. Например, Mag X[0] и Phase X[0] вычисляются с использованием только ReX[0] и ImX[0]. Подобным образом Mag X[14] и Phase X[14] вычисляются с использованием только ReX[14] и ImX[14] и так далее. Для понимания этих преобразований, посмотрим, что происходит, когда вы складываете косинусную волну и синусную волну в частотной области. Результат есть косинусная волна той же самой частоты, но с новой амплитудой и новым фазовым сдвигом. В виде формулы это выглядит так:
A cos(x) + B sin(x) = M cos(x + )
Важно, что информация в этом процессе не теряется; имея одно представление, вы можете вычислить другое. Другими словами, информация содержащаяся в амплитудах А и В, также содержится в величинах М и . Хотя формула содержит косинусы и синусы, она полностью соответствует преобразованию простых векторов. Рисунок 8-9 показывает векторную аналогию, переменные А и В могут рассматриваться в прямоугольной системе координат, в то время как переменные М и есть параметры полярной системы координат.
В полярной системе координат Mag X[ ] является амплитудой косинусной волны (М в формуле 8-5 и рисунка 8-9), Phase X[ ] фаза угла косинусной волны ( в формуле 8-5 и рисунка 8-9). Следующие формулы обеспечивают преобразование прямоугольной системы координат в полярную систему координат и наоборот.
Mag X[k] = ( Re X[k]2 + Im X[k]2 )1/2
Re X[k] = Mag X[k] cos(Phase X[k])
Im X[k] = Mag X[k] sin(Phase X[k])
Прямоугольная система координат и полярная система координат позволяют нам рассматривать ДПФ с двух точек зрения. В полярной системе координат ДПФ раскладывает N-точечный сигнал в N/2 + 1 косинусных волн и N/2 + 1 синусных волн, каждая со своей амплитудой. В полярной системе координат ДПФ раскладывает N-точечный сигнал в N/2 + 1 косинусных волн, каждая из которых определяется амплитудой и фазовым сдвигом. Почему в полярной системе координат используются косинусы, а не синусы? Синусные волны не могут представить постоянную составляющую в сигнале, поскольку синусная волна нулевой частоты состоит из одних нулей (см. рисунок 8-5а,b).
РИСУНОК 8-10
Пример отображения частотной области в прямоугольной и полярной системе координат. Полярная система координат обычно более наглядна и понятна для отображения характеристик сигнала. Прямоугольная система координат почти всегда применяется при компьютерных вычислениях. Обратите особое внимание на то, что первый и последний отсчет по фазе должен быть равен нулю, также как в мнимой части.
Хотя полярная и прямоугольная системы координат представляют одну и туже информацию, имеется много примеров, где одна из них легче в использовании, чем другая. Например, рисунок 8-10 показывает сигналы частотной области в обеих системах координат. Предупреждение: не пытайтесь понять форму реальной и мнимой частей, ваша голова лопнет! Для сравнения, кривая в полярной системе координат наглядна, из нее видно, что присутствуют только частоты ниже 0,25 и фазовый сдвиг почти пропорционален частоте. Это частотная характеристика (частотный отклик) низкочастотного (low-pass) фильтра.
Когда мы должны использовать прямоугольную систему координат, а когда полярную? Прямоугольная система координат лучше всего подходит для вычислений, как в формульном виде, так и в компьютерной программе. Графики почти всегда отображаются в полярной системе координат. Как показывает предыдущий пример это более наглядно для человеческого восприятия, чем реальная и мнимая части. В типичной программе сигнал частотной области хранится в прямоугольной системе координат до момента его вывода для наблюдения, тогда он преобразуется в полярную систему координат.
Почему легче понять частотную область в полярной системе координат? Этот вопрос является сердцем другого вопроса: почему полезно разложение сигнала на синусоиды? Вспомните свойство синусоидальной точности воспроизведения (sinusoidal fidelity) из главы 5: если на вход линейной системы подать синусоиду, то на выходе системы будет тоже синусоида точно такой же частоты. Только амплитуда и фаза могут быть изменены. Полярная система координат непосредственно представляет амплитуду и фазу косинусных волн, на которые раскладывается сигнал. В свою очередь система может быть представлена тем, как она преобразует амплитуду и фазу каждой косинусной волны.
Теперь посмотрим, что случится, если мы будем использовать прямоугольную систему координат для этого сценария. Смесь косинусных и синусных волн на входе линейной системы преобразуется в смесь косинусных и синусных волн на ее выходе. Проблема заключается в том, что косинусная волна на входе системы может преобразоваться и в синусную и в косинусную волну. Точно также и синусная волна может преобразоваться и в синусную и в косинусную волну. Поскольку эта путаница может быть исправлена, общий подход не объясняет, почему мы желаем, прежде всего, использовать синусоиду.
Особенности полярной системы координат
Имеется много особенностей, связанных с использованием полярной системы координат. Ни одна из них не смертельна, но все равно досаждают! Таблица 8-3 показывает компьютерную программу преобразования между прямоугольной системой координат и полярной системой координат, и дает решение некоторых из этих неприятностей.
Особенность первая. Радианы или градусы?
Фазу можно выразить либо в градусах, либо в радианах. Когда она выражается в градусах, то находится между минус 180 и 180. При использовании радиан – между - и , то есть, между –3,141592 и 3,141592. Большинство компьютерных языков требует использования радиан для тригонометрических функций, таких как косинус, синус, арктангенс и т. д. Работа с длинными децимальными числами и трудная интерпретация полученных данных могут сильно раздражать. Например, если вы хотите ввести 90 градусный фазовый сдвиг, вам необходимо добавить в фазу 1,570796. Хотя это и не убьет вас, но это надоедает. Лучшее решение этой проблемы заключается в определении постоянной PI=3,141592 в начале вашей программы. 90 градусов могут быть тогда записаны как PI/2. И градусы, и радианы широко используются в ЦОС, и вам необходимо разобраться с ними.
Особенность 2. Деление на ноль.
Когда преобразуется прямоугольная система координат в полярную, очень часто встречаются частоты с реальной частью равной нулю и мнимой частью не равной нулю. Это значит, что фаза точно равна 90 или -90 градусов. Но попробуйте объяснить это вашему компьютеру! Когда ваша программа пытается вычислить фазу по формуле: PhaseX[k] = arctan(ImX[k]/ReX[k]), появляется ошибка деления на ноль. Даже если ваша программа не остановится, фаза для этой частоты будет неверной. Чтобы избежать эту проблему, перед делением необходимо проверить реальную часть на наличие нулей. Если ноль обнаружен, то необходимо проверить мнимую часть на положительное и отрицательное значения для определения значения фазы /2 или -/2, соответственно. Последнее, деление необходимо обойти. Ничего сложного в этих шагах нет, кроме возможности ошибиться. Альтернативный путь решения этой проблемы показан в линии 250 таблицы 8-3. Если реальная часть равна нулю, замените ее очень маленькой величиной, чтобы осчастливить математический процессор во время деления.
100 'ПРЕОБРАЗОВАНИЕ СИСТЕМ КООРДИНАТ
110 '
120 DIM REX[256] 'REX[ ] содержит реальную часть
130 DIM IMX[256] 'IMX[ ] содержит мнимую часть
140 DIM MAG[256] 'MAG[ ] содержит амплитуду (magnitude)
150 DIM PHASE[256] 'PHASE[ ] содержит фазу
160 '
170 PI = 3.14159265
180 '
190 GOSUB XXXX 'Воображаемая подпрограмма загрузки REX[ ] и IMX[ ]
200 '
210 '
220 ' 'Прямоугольную в полярную, формула 8-6
230 FOR K% = 0 TO 256
240 MAG[K%] = SQR( REX[K%]^2 + IMX[K%]^2 ) 'из формулы 8-6
250 IF REX[K%] = 0 THEN REX[K%] = 1E-20 'обход деления на 0 (особенность 2)
260 PHASE[K%] = ATN( IMX[K%] / REX[K%] ) 'из формулы 8-6
270 ' 'корректировка арктангенса (особенность 3)
280 IF REX[K%] < 0 AND IMX[K%] < 0 THEN PHASE[K%] = PHASE[K%] - PI
290 IF REX[K%] < 0 AND IMX[K%] >= 0 THEN PHASE[K%] = PHASE[K%] + PI
300 NEXT K%
310 '
320 '
330 ' Полярную в прямоугольную, формула 8-7
340 FOR K% = 0 TO 256
350 REX[K%] = MAG[K%] * COS( PHASE[K%] )
360 IMX[K%] = MAG[K%] * SIN( PHASE[K%] )
370 NEXT K%
380 '
390 END
ТАБЛИЦА 8-3
Особенность 3. Некорректный арктангенс.
Представьте отсчет частотной области, у которого ReX[k] = 1 и ImX[k] = 1. Формула 8-6 даст соответствующую полярную величину MagX[k] = 1,414 и ReX[k] = 45. Теперь рассмотрим другой отсчет, у которого ReX[k] = -1 и ImX[k] = -1. Снова формулф8-6 даст MagX[k] = 1,414 и ReX[k] = 45. Но фаза неверная! Она должна быть -135. Это ошибка появляется, когда реальная часть отрицательная. Проблема может быть исправлена тестированием реальной и мнимой части после вычисления фазы. Если реальная и мнимая части отрицательные, то надо вычесть 180 (или радиан) из вычисленной фазы. Если реальная часть отрицательная, а мнимая положительная, то добавить180 (или радиан). Линии 340и 350 в программе таблицы 8-3 показывают, как это делается. Если вы упустили эту проблему, то вычисленные фазы будут находиться между -/2 и /2, вместо - и . Муштруйте свой ум! Если вы видите фазу только между 1,5708, то вы забыли скорректировать неоднозначность при вычислении арктангенса.
Особенность 4. Фаза очень маленькой амплитуды.
Представьте следующую ситуацию. Вы разжевали некоторую задачу ЦОС, и неожиданно обнаружили, что часть фазы выглядит неправильно. Это может быть шум, скачки или плавная ошибка. После нескольких часов разглядывания компьютерного кода вы находите ответ. Соответствующая величина амплитуды слишком мала, и она теряется в шуме округления. Если амплитуда чрезвычайно мала, фаза не имеет значения и может принимать необычную величину. Пример этого показан на рисунке 8-11. Это обычно очевидно, когда амплитуда теряется в шумах; когда величина ее слишком мала, то вы усиливаетесь в подозрениях, что ее величина бессмысленна. Фаза меняется. Когда полярный сигнал загрязнен шумами, величина фазы принимает случайные значения между - и . К несчастью, это часто больше похоже на реальный сигнал, чем на вздор.
РИСУНОК 8-11
Фаза сигнала с маленькой амплитудой. Для частот, у которых амплитуда падает до очень низкой величины, шум округления вызывает «дикие» броски фаз. Не ошибайтесь, предполагая, что этот сигнал имеет значение.
Особенность 5. 2 неоднозначность фазы.
Посмотрите снова на рисунок 8-10d, и обратите внимание на прерывистость данных. Каждый раз точка выглядит так, как будто она собирается нырнуть ниже -3,141592 и резко прыгает к 3,141592. Это есть результат периодической природы синусоиды. Например, фазовый сдвиг есть то же самое, что и фазовый сдвиг + 2, + 4, + 6 и т.д. Любая синусоида не изменяется, когда вы добавляете к фазе число кратное 2. Видимая прерывистость в сигнале есть результат компьютерных алгоритмов, осуществляющих предпочтительный отбор из бесконечного числа эквивалентных возможностей. Всегда выбирается самая маленькая величина, сохраняющую фазу между - и .
Очень часто легче понять фазу, если нет этой прерывистости, даже если это означает, что фаза будет превышать или буде ниже -. Это называется разверткой фазы (unwrapping the phase), и пример этого показан на рисунке 8 12. Как показано в программе таблицы 8-4, 2 многократно добавляется или вычитается из каждой величины фазы. Точная величина фазы определяется алгоритмом, который минимизирует разницу между соседними отсчетами.
РИСУНОК 8-12
Пример фазовой развертки. Верхняя кривая показывает типичный сигнал фазы, полученный при преобразовании прямоугольной системы координат в полярную. Каждая величина сигнала должна находится между - и . Нижняя кривая показывает развернутую фазу.
100 ' РАЗВЕРТКА ФАЗЫ
110 '
120 DIM PHASE[256] 'PHASE[ ] содержит исходную фазу
130 DIM UWPHASE[256] 'UWPHASE[ ] содержит развернутую фазу
140 '
150 PI = 3.14159265
160 '
170 GOSUB XXXX 'Воображаемая подпрограмма загрузки PHASE[ ]
180 '
190 UWPHASE[0] = 0 'Первая точка всех фаз, равная нулю
200 '
210 ' 'Алгоритм развертки
220 FOR K% = 1 TO 256
230 C% = CINT( (UWPHASE[K%-1] - PHASE[K%]) / (2 * PI) )
240 UWPHASE[K%] = PHASE[K%] + C%*2*PI
250 NEXT K%
260 '
270 END
ТАБЛИЦА 8-4
Особенность 6. Амплитуда всегда положительная ( неоднозначность фазы)
Рисунок 8-13 показывает сигнал частотной области в прямоугольной и полярной форме. Реальная часть гладкая и вполне легкая для понимания, мнимая часть целиком состоит из нулей. Сигнал в полярной форме содержит резкие изломы и острые углы. Это вызвано тем, что амплитуда в полярной форме (magnitude) должна быть положительной по определению. Хотя реальная часть падает ниже нуля, амплитуда в полярной форме остается положительной за счет изменения фазы на (или -, что одно и тоже). Это не проблема с точки зрения математики, но нерегулярная кривая может быть более трудной для понимания.
РИСУНОК 8-13
Пример сигнала в прямоугольной и полярной форме. Так как амплитуда всегда должна быть положительной (по определению), амплитуда и фаза могут содержать резкие изломы и острые углы. Рисунок (d) показывает другой нюанс: случайный шум может вызывать быстрые скачки фазы между и -.
Некоторые решения позволяют иметь полярной амплитуде (magnitude) отрицательные величины. В примере на рисунке 8-13, полярная амплитуда может выглядеть точно так же, как и реальная часть при полностью нулевой фазе. В этом не будет ошибки, если это поможет вашему пониманию. Только постарайтесь не называть сигнал с отрицательными величинами полярной амплитудой, чтобы не осквернить формальное определение. В этой книге мы будем использовать слово «развернутая амплитуда» (unwrapped magnitude), чтобы указать, что «амплитуде» (magnitude) разрешается иметь отрицательные величины.
Особенность 7. Выбросы между и -.
Так как и - представляют одну и ту же фазу, шум округления может быть причиной резкого скачка фазы между этими величинами у соседних отсчетов. Как показано на рисунке 8-13d, это может приводить к острым изломам и выбросам на гладкой во всем остальном кривой. Не позволяйте себя путать, фаза в действительности непрерывна.