Дискретное двумерное преобразование фурье c. Описание Image Processing Toolbox

Я полагаю что все в общих чертах знают о существовании такого замечательного математического инструмента как преобразование Фурье. Однако в ВУЗах его почему-то преподают настолько плохо, что понимают как это преобразование работает и как им правильно следует пользоваться сравнительно немного людей. Между тем математика данного преобразования на удивление красива, проста и изящна. Я предлагаю всем желающим узнать немного больше о преобразовании Фурье и близкой ему теме того как аналоговые сигналы удается эффективно превращать для вычислительной обработки в цифровые.

Без использования сложных формул и матлаба я постараюсь ответить на следующие вопросы:

  • FT, DTF, DTFT - в чем отличия и как совершенно разные казалось бы формулы дают столь концептуально похожие результаты?
  • Как правильно интерпретировать результаты быстрого преобразования Фурье (FFT)
  • Что делать если дан сигнал из 179 сэмплов а БПФ требует на вход последовательность по длине равную степени двойки
  • Почему при попытке получить с помощью Фурье спектр синусоиды вместо ожидаемой одиночной “палки” на графике вылезает странная загогулина и что с этим можно сделать
  • Зачем перед АЦП и после ЦАП ставят аналоговые фильтры
  • Можно ли оцифровать АЦП сигнал с частотой выше половины частоты дискретизации (школьный ответ неверен, правильный ответ - можно)
  • Как по цифровой последовательности восстанавливают исходный сигнал

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

Начать надо, наверное, с того что обычное преобразование Фурье - это некая такая штука которая, как можно догадаться из названия, преобразует одни функции в другие, то есть ставит в соответствие каждой функции действительного переменного x(t) её спектр или фурье-образ y(w):

Если приводить аналогии, то примером аналогичного по смыслу преобразования может послужить например дифференцирование, превращающее функцию в её производную. То есть преобразование Фурье - такая же, по сути, операция как и взятие производной, и её часто обозначают схожим образом, рисуя треугольную “шапочку” над функцией. Только в отличие от дифференцирования которое можно определить и для действительных чисел, преобразование Фурье всегда “работает” с более общими комплексными числами. Из-за этого постоянно возникают проблемы с отображением результатов этого преобразования, поскольку комплексные числа определяются не одной, а двумя координатами на оперирующем действительными числами графике. Удобнее всего, как правило, оказывается представить комплексные числа в виде модуля и аргумента и нарисовать их по раздельности как два отдельных графика:

График аргумента комплексного значения часто называют в данном случае “фазовым спектром”, а график модуля - “амплитудным спектром”. Амплитудный спектр как правило представляет намного больший интерес, а потому “фазовую” часть спектра нередко пропускают. В этой статье мы тоже сосредоточимся на “амплитудных” вещах, но забывать про существование пропущенной фазовой части графика не следует. Кроме того, вместо обычного модуля комплексного значения часто рисуют его десятичный логарифм умноженный на 10. В результате получается логарифмический график, значения на котором отображаются в децибелах (дБ).

Обратите внимание что не очень сильно отрицательным числам логарифмического графика (-20 дБ и менее) при этом соответствуют практически нулевые числа на графике “обычном”. Поэтому длинные и широкие “хвосты” разнообразных спектров на таких графиках при отображении в “обычные” координаты как правило практически исчезают. Удобство подобного странного на первый взгляд представления возникает из того что фурье-образы различных функций часто необходимо перемножать между собой. При подобном поточечном умножении комплекснозначных фурье-образов их фазовые спектры складываются, а амплитудные - перемножаются. Первое выполняется легко, а второе - сравнительно сложно. Однако логарифмы амплитуды при перемножении амплитуд складываются, поэтому логарифмические графики амплитуды можно, как и графики фаз, просто поточечно складывать. Кроме того, в практических задачах часто удобнее оперировать не «амплитудой» сигнала, а его «мощностью» (квадратом амплитуды). На логарифмической шкале оба графика (и амплитуды и мощности) выглядят идентично и отличаются только коэффициентом - все значения на графике мощности ровно вдвое больше чем на шкале амплитуд. Соответственно для построения графика распределения мощности по частоте (в децибелах) можно не возводить ничего в квадрат, а посчитать десятичный логарифм и умножить его на 20.

Заскучали? Погодите, еще немного, с занудной частью статьи, объясняющей как интерпретировать графики, мы скоро покончим:). Но перед этим следует понять одну крайне важную вещь: хотя все вышеприведенные графики спектров были нарисованы для некоторых ограниченных диапазонов значений (в частности, положительных чисел), все эти графики на самом деле продолжаются в плюс и минус бесконечность. На графиках просто изображается некоторая “наиболее содержательная” часть графика, которая обычно зеркально отражается для отрицательных значений параметра и зачастую периодически повторяется с некоторым шагом, если рассматривать её в более крупном масштабе.

Определившись с тем, что же рисуется на графиках, давайте вернемся собственно к преобразованию Фурье и его свойствам. Существует несколько разных способов как определить это преобразование, отличающихся небольшими деталями (разными нормировками). Например в наших ВУЗах почему-то часто используют нормировку преобразования Фурье определяющую спектр в терминах угловой частоты (радианов в секунду). Я буду использовать более удобную западную формулировку, определяющую спектр в терминах обычной частоты (герцах). Прямое и обратное преобразование Фурье в этом случае определяются формулами слева, а некоторые свойства этого преобразования которые нам понадобятся - списком из семи пунктов справа:

Первое из этих свойств - линейность. Если мы берем какую-то линейную комбинацию функций, то преобразование Фурье этой комбинации будет такой же линейной комбинацией образов Фурье этих функций. Это свойство позволяет сводить сложные функции и их фурье-образы к более простым. Например, фурье-образ синусоидальной функции с частотой f и амплитудой a является комбинацией из двух дельта-функций расположенных в точках f и -f и с коэффициентом a/2:

Если взять функцию, состоящую из суммы множества синусоид с разными частотами, то согласно свойству линейности, фурье-образ этой функции будет состоять из соответствующего набора дельта-функций. Это позволяет дать наивную, но наглядную интерпретацию спектра по принципу “если в спектре функции частоте f соответствует амплитуда a, то исходную функцию можно представить как сумму синусоид, одной из которых будет синусоида с частотой f и амплитудой 2a”. Строго говоря, эта интерпретация неверна, поскольку дельта-функция и точка на графике - это совершенно разные вещи, но как мы увидим дальше, для дискретных преобразований Фурье она будет не так уж и далека от истины.

Второе свойство преобразования Фурье - это независимость амплитудного спектра от сдвига сигнала по времени. Если мы подвинем функцию влево или вправо по оси x, то поменяется лишь её фазовый спектр.

Третье свойство - растяжение (сжатие) исходной функции по оси времени (x) пропорционально сжимает (растягивает) её фурье-образ по шкале частот (w). В частности, спектр сигнала конечной длительности всегда бесконечно широк и наоборот, спектр конечной ширины всегда соответствует сигналу неограниченной длительности.

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

Шестое свойство говорит о симметрии фурье-образов. В частности, из этого свойства следует что в фурье-образе действительнозначной функции (т.е. любого “реального” сигнала) амплитудный спектр всегда является четной функцией, а фазовый спектр (если его привести к диапазону -pi...pi) - нечетной. Именно по этой причине на графиках спектров практически никогда не рисуют отрицательную часть спектра - для действительнозначных сигналов она не дает никакой новой информации (но, повторюсь, и нулевой при этом не является).

Наконец последнее, седьмое свойство, говорит о том, что преобразование Фурье сохраняет “энергию” сигнала. Оно осмысленно только для сигналов конечной продолжительности, энергия которых конечна, и говорит о том, что спектр подобных сигналов на бесконечности быстро приближается к нулю. Именно в силу этого свойства на графиках спектров как правило изображают только “основную” часть сигнала, несущую в себе львиную долю энергии - остальная часть графика просто стремится к нулю (но, опять же, нулем не является).

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

Гребенка Дирака - это просто периодическая последовательность дельта-функций с единичным коэффициентом, начинающаяся в нуле и идущая с шагом T. Для оцифровки сигналов, T выбирают по возможности малым числом, T<<1. Фурье-образ этой функции - тоже гребенка Дирака, только с гораздо большим шагом 1/T и несколько меньшим коэффициентом (1/T). С математической точки зрения, дискретизация сигнала по времени - это просто поточечное умножение исходного сигнала на гребенку Дирака. Значение 1/T при этом называют частотой дискретизации:

Вместо непрерывной функции после подобного перемножения получается последовательность дельта-импульсов определенной высоты. При этом согласно свойству 5 преобразования Фурье, спектр получившегося дискретного сигнала есть свертка исходного спектра с соответствующей гребенкой Дирака. Несложно понять, что исходя из свойств свертки, спектр исходного сигнала при этом как бы “копируется” бесконечное число раз вдоль оси частот с шагом 1/T, а затем суммируется.

Заметим, что если исходный спектр имел конечную ширину и мы использовали достаточно большую частоту дискретизации, то копии исходного спектра не будут перекрываться, а следовательно и суммироваться друг с другом. Несложно понять что по подобному “свернутому” спектру будет легко восстановить исходный - достаточно будет просто взять компоненту спектра в районе нуля, “обрезав” лишние копии уходящие на бесконечность. Простейший способ это сделать - это домножить спектр на прямоугольную функцию, равную T в диапазоне -1/2T...1/2T и нулю - вне этого диапазона. Подобный Фурье-образ соответствует функции sinc (Tx) и согласно свойству 4, подобное умножение равнозначно свертке исходной последовательности дельта-функций с функцией sinc(Tx)



То есть с помощью преобразования Фурье мы получили способ легко восстановить исходный сигнал из дискретизированного по времени, работающий при условии что мы используем частоту дискретизации, по крайней мере вдвое (из-за наличия в спектре отрицательных частот) превышающую максимальную частоту присутствующую в исходном сигнале. Этот результат широко известен и называется “теорема Котельникова / Шеннона-Найквиста” . Однако, как несложно теперь (понимая доказательство) заметить, этот результат вопреки широко распространенному заблуждению определяет достаточное , но не необходимое условие для восстановления исходного сигнала. Все что нам требуется - это добиться того, чтобы интересующая нас часть спектра после дискретизации сигнала не накладывалась друг на друга и если сигнал достаточно узкополосный (имеет малую “ширину” ненулевой части спектра), то этого результата часто можно добиться и при частоте дискретизации намного ниже чем удвоенная максимальная частота сигнале. Подобная техника называется “undersampling” (субдискретизация, полосовая дискретизация) и довольно широко используется при обработке всевозможных радиосигналов. Например, если мы берем FM-радио действующее в полосе частот от 88 до 108 МГц, то для его оцифровки можно использовать АЦП с частотой всего 43.5 МГц вместо предполагающихся по теореме Котельникова 216 МГц. При этом, правда, понадобится качественный АЦП и хороший фильтр.

Замечу, что “дублирование” высоких частот частотами меньших порядков (алиасинг) - непосредственное свойство дискретизации сигнала, необратимо “портящее” результат. Поэтому если в сигнале в принципе могут присутствовать частоты высокого порядка (то есть практически всегда) перед АЦП ставят аналоговый фильтр, “отсекающий” все лишнее непосредственно в исходном сигнале (так как после дискретизации делать это уже будет поздно). Характеристики этих фильтров, как аналоговых устройств, неидеальны, поэтому некоторая “порча” сигнала при этом все равно происходит, и на практике из этого следует что наибольшие частоты в спектре, как правило, недостоверны. Чтобы уменьшить эту проблему, сигнал нередко сэмплируют с завышенной частотой дискретизации, ставя при этом входной аналоговый фильтр на меньшую полосу пропускания и используя только нижнюю часть теоретически доступного частотного диапазона АЦП.

Еще одно распространенное заблуждение, кстати, - это когда сигнал на выходе ЦАП рисуют “ступеньками”. “Ступеньки” соответствуют свертке дискретизированной последовательности сигналов с прямоугольной функцией ширины T и высоты 1:

Спектр сигнала при таком преобразовании умножается на фурье-образ этой прямоугольной функции, а у подобной прямоугольной функции это снова sinc(w), “растянутый” тем сильнее, чем меньше ширина соответствующего прямоугольника. Спектр дискретизированного сигнала при подобном “ЦАП” поточечно умножается на этот спектр. При этом ненужные высокие частоты с “лишними копиями” спектра обрезаются не полностью, а верхняя часть “полезной” части спектра, напротив, ослабляется.

На практике так, естественно, никто не делает. Существует много разных подходов к построению ЦАП, но даже в наиболее близких по смыслу ЦАП взвешивающего типа прямоугольные импульсы в ЦАП напротив выбираются по возможности короткими (приближающимися к настоящей последовательности дельта-функций) чтобы избежать излишнего подавления полезной части спектра. “Лишние” частоты в получившемся широкополосном сигнале практически всегда гасят, пропуская сигнал через аналоговый фильтр низких частот, так что «цифровых ступенек» нет ни «внутри» преобразователя, ни, тем более, на его выходе.

Однако вернемся обратно к преобразованию Фурье. Описанное выше преобразование Фурье, примененное к заранее дискретизированной последовательности сигналов называется преобразованием Фурье дискретного времени (DTFT). Спектр получаемый подобным преобразованием всегда 1/T-периодичен, поэтому спектр DTFT полностью определяется её значениями на отрезке , n=0,…,N-1 - исходный комплексный сигнал, состоящий из N комплексных чисел. Обозначим X[k], k=0,…N-1 - его комплексный спектр, также состоящий из N комплексных чисел. Тогда справедливы следующие формулы прямого и обратного преобразований Фурье:

Если по этим формулам разложить в спектр действительный сигнал, то первые N/2+1 комплексных коэффициентов спектра будут совпадать со спектром "обычного" действительного ДПФ, представленным в "комплексном" виде, а остальные коэффициенты будут их симметричным отражением относительно половины частоты дискретизации. Для косинусных коэффициентов отражение четное, а для синусных - нечетное.

Двумерное ДПФ

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

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

Здесь N 1 xN 2 - размер исходного сигнала, он же - размер спектра. k 1 и k 2 - это номера базисных функций (номера коэффициентов двумерного ДПФ, при которых эти функции находятся). Поскольку размер спектра равен размеру исходного сигнала, то k 1 = 0,…,N 1 -1; k 2 = 0,…,N 2 -1.

n 1 и n 2 - переменные-аргументы базисных функций. Поскольку область определения базисных функций совпадает с областью определения сигнала, то n 1 = 0,…,N 1 -1; n 2 = 0,…,N 2 -1.

Двумерное ДПФ (в комплексной форме) определяется следующими формулами (здесь x - исходный сигнал, а X - его спектр):

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

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

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

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

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

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

(2)

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

(3)

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

(4)

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

, (5)

Только функции sin и cos c частотами k/T будут взаимно ортогональны, а это является непременным условием преобразования Фурье. Набор первых функций разложения в ряд Фурье приведен на рисунке 1. При этом длительность функций совпадает с длительностью анализа T .


Рисунок 1. Функции разложения в ряд Фурье

Теперь спектр сигнала будет выглядеть так, как это показано на рисунке 2.



Рисунок 2. Спектр функции x (t ) при анализе на ограниченном интервале времени

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

(6)

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

(7)

Подобным образом можно определить формулу прямого преобразования Фурье для цифровых отсчетов сигнала. Учитывая, что вместо непрерывного сигнала используются его цифровые отсчеты, в выражении (6) интеграл заменяется на сумму. В данном случае длительность анализируемого сигнала определяется количеством цифровых отсчетов N . Преобразование Фурье для цифровых отсчетов сигнала называется дискретным преобразованием Фурье и записывается следующим образом:

(8)

Теперь рассмотрим как изменились свойства дискретного преобразования Фурье (ДПФ) по сравнению с прямым преобразованием Фурье на ограниченном интервале времени. Когда мы рассматривали дискретизацию аналогового сигнала, мы выяснили, что спектр входного сигнала должен быть ограничен по частоте. Это требование ограничивает количество дискретных составляющих спектра сигнала. Первоначально может показаться, что мы можем ограничить спектр сигнала частотой f д /2, что соответствует количеству частотных составляющих K = N /2 . Однако это не так. Несмотря на то, что спектр сигнала для действительных отсчетов сигнала для положительных частот и отрицательных частот симметричен относительно 0, отрицательные частоты могут потребоваться для некоторых алгоритмов работы со спектрами, например, для . Еще больше отличие получается при выполнении дискретного преобразования Фурье над комплексными отсчетами входного сигнала. В результате для полного описания спектра цифрового сигнала требуется N частотных отсчетов (k = 0, ..., N/2 ).

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

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

Для расчёта преобразования Фурье обычно используется ускоренная процедура расчёта - т.н. быстрое преобразование Фурье (БПФ). Это позволяет в значительной мере сократить процессорное время на достаточно сложные и ресурсоёмкие математические расчёты.

1 Комплексные числа

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

Код класса для описания комплексных чисел (разворачивается) """ """ Комплексное число. """ Public Class ComplexNumber """ """ Действительная часть комплексного числа. """ Public Real As Double = 0 """ """ Мнимая часть комплексного числа. """ Public Imaginary As Double = 0 Public Sub New() Real = 0 Imaginary = 0 End Sub """ """ Создаёт комплексное число. """ """ Действительная часть комплексного числа. """ Мнимая часть комплексного числа. Public Sub New(ByVal r As Double, Optional ByVal im As Double = 0) Real = r Imaginary = im End Sub Private usCult As New Globalization.CultureInfo("en-US") "используем культуру "en-US" чтобы целая и дробная части разделялись точкой, а не запятой """ """ Возвращает строку, состоящую из действительной и мнимой части, разделённых символом табуляции. """ Public Overrides Function ToString() As String Return (Real.ToString(usCult) & ControlChars.Tab & Imaginary.ToString(usCult)) End Function End Class

2 Прямое дискретное быстрое преобразование Фурье

На вход функции передаётся массив комплексных чисел. Действительная часть которого представляет произвольный дискретный сигнал, с отсчётами через равные промежутки времени. Мнимая часть содержит нули. Число отсчётов в сигнале должно равняться степени двойки. Если ваш сигнал короче, то дополните его нулями до числа, кратного степени 2: 256, 512, 1024 и т.д. Чем длиннее сигнал, тем у рассчитанного спектра будет выше разрешение по частоте.

Код для расчёта прямого быстрого преобразования Фурье на VB.NET (разворачивается) """ """ Рассчитывает спектр сигнала методом быстрого преобразования Фурье. Использовать только (N/2+1) возвращаемых значений (до половины частоты дискретизации). """ """ Сигнал, содержащий количество отсчётов, кратное степени двойки, и состоящий из действительной и мнимой частей. Все мнимые части сигнала заполнены нулями. """ Возвращает массив комплексных чисел спектра. """ Значимы только первые N/2+1, остальные - симметричная часть, соответствующая отрицательным частотам. """ Первое значение спектра - это постоянная составляющая, последнее - соответствует половине частоты дискретизации (частота Найквиста). """ Значения выше половины частоты дискретизации - не использовать. """ Public Shared Function FFT(ByVal signal As ComplexNumber()) As ComplexNumber() Dim order As Integer = signal.Length "порядок ДПФ CheckFftOrder(order) "Проверяем, что порядок равен степени двойки Dim spectrumLen As Integer = order \ 2 Dim j As Integer = spectrumLen "Бит-реверсная сортировка: For i As Integer = 1 To order - 2 If (i < j) Then Dim tmpRe As Double = signal(j).Real Dim tmpIm As Double = signal(j).Imaginary signal(j).Real = signal(i).Real signal(j).Imaginary = signal(i).Imaginary signal(i).Real = tmpRe signal(i).Imaginary = tmpIm End If Dim k As Integer = spectrumLen Do Until (k > j) j -= k k \= 2 Loop j += k Next "Цикл по уровням разложения: For level As Integer = 1 To CInt(Math.Log(order) / Math.Log(2)) Dim lvl As Integer = CInt(2 ^ level) Dim lvl2 As Integer = lvl \ 2 Dim tmp As Double = Math.PI / lvl2 Dim sr As Double = Math.Cos(tmp) Dim si As Double = -Math.Sin(tmp) Dim tr As Double = 0 Dim ur As Double = 1 Dim ui As Double = 0 For jj As Integer = 1 To lvl2 "Цикл по спектрам внутри уровня For i As Integer = (jj - 1) To (order - 1) Step lvl "Цикл по отдельным "бабочкам" Dim ip As Integer = i + lvl2 tr = signal(ip).Real * ur - signal(ip).Imaginary * ui "Операция "бабочка" Dim ti As Double = signal(ip).Real * ui + signal(ip).Imaginary * ur signal(ip).Real = signal(i).Real - tr signal(ip).Imaginary = signal(i).Imaginary - ti signal(i).Real = signal(i).Real + tr signal(i).Imaginary = signal(i).Imaginary + ti Next tr = ur ur = tr * sr - ui * si ui = tr * si + ui * sr Next Next "Заполняем массив комплексных чисел, обработанных БПФ: Dim spectrum(order - 1) As ComplexNumber For i As Integer = 0 To order - 1 With signal(i) spectrum(i) = New ComplexNumber(.Real, .Imaginary) End With Next Return spectrum End Function

3 Обратное дискретное быстрое преобразование Фурье

Обратное дискретное преобразование Фурье (ОДПФ) одним из этапов расчёта включает в себя прямое ДПФ на массиве комплексных чисел, где мнимая часть - это инверсия относительно оси X мнимой части спектра.

Код для расчёта обратного быстрого преобразования Фурье на VB.NET (разворачивается) """ """ Восстанавливает сигнал по его спектру методом обратного быстрого преобразования Фурье. """ """ Спектр сигнала, содержащий количество отсчётов, кратное степени двойки, и состоящий из действительной и мнимой частей. Public Shared Function InverseFFT(ByVal spectrum As ComplexNumber()) As ComplexNumber() Dim order As Integer = spectrum.Length "Порядок обратного ДПФ. CheckFftOrder(order) "Изменение арифметического знака элементов мнимой части: For i As Integer = 0 To spectrum.Length - 1 spectrum(i).Imaginary = -spectrum(i).Imaginary Next "Вычисление прямого БПФ: Dim directFFT As ComplexNumber() = FFT(spectrum) "Деление на order во временной области со сменой арифметического знака мнимой части: Dim signal(directFFT.Length - 1) As ComplexNumber For i As Integer = 0 To directFFT.Length - 1 Dim ReX As Double = directFFT(i).Real / order Dim ImX As Double = -directFFT(i).Imaginary / order signal(i) = New ComplexNumber(ReX, ImX) Next Return signal End Function

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

"""

""" Проверяет, является ли порядок БПФ степенью двойки, и если нет - вызывает исключение. """ """ Порядок БПФ. Private Shared Sub CheckFftOrder(ByVal order As Integer) Dim chk As Double = Math.Abs(Math.Floor(Math.Log(order, 2)) - Math.Log(order, 2)) If (chk > 0.0001) Then Throw New ArgumentException(String.Format("Длина массива ({0}) не кратна степени двойки.", order)) End If End Sub

4 Проверка прямого и обратного преобразования Фурье

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

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

Dim cn(127) As ComplexNumber For i As Integer = 0 To cn.Length - 1 cn(i) = New ComplexNumber(Math.Sin(i * 3 * Math.PI / 180)) Next

Получим вот такой сигнал:

Здесь по оси X - номера отсчётов во временной области, по оси Y - амплитуда. Обратим внимание, что сигнал состоит только из действительных частей, а мнимая часть на всём отрезке равна "0".

Теперь передадим этот сигнал на вход функции FFT(). По полученным в ходе прямого преобразования Фурье массивам комплексных чисел построим два графика - действительной (Re) и мнимой (Im) частей спектра:


Здесь по оси X - отсчёты в частотной области, по оси Y - амплитуда. Чтобы получить реальные значения частоты, необходимо рассчитать их, учитывая, что "0" оси Y соответствует нулевой частоте, максимум оси Y соответствует частоте дискретизации.

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


Как видно, восстановленный сигнал полностью повторяет исходный.

Преобразования Фурье

Многие сигналы удобно анализировать, раскладывая их на синусоиды (гармоники). Тому есть несколько причин. Например, подобным образом работает человеческое ухо. Оно раскладывает звук на отдельные колебания различных частот. Кроме того, можно показать, что синусоиды являются «собственными функциями» линейных систем (т.к. они проходят через линейные системы, не изменяя формы, а могут изменять лишь фазу и амплитуду). Еще одна причина в том, что теорема Котельникова формулируется в терминах спектра сигнала.

Преобразование Фурье (Fourier transform)– это разложение функций на синусоиды (далее косинусные функции мы тоже называем синусоидами, т.к. они отличаются от «настоящих» синусоид только фазой). Существует несколько видов преобразования Фурье.

1. Непериодический непрерывный сигнал можно разложить в интеграл Фурье.

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

3. Непериодический дискретный сигнал можно разложить в интеграл Фурье.

4. Периодический дискретный сигнал можно разложить в конечный ряд Фурье.

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

ДПФ вещественного сигнала

Пусть дискретный сигнал x имеет периодN точек. В этом случае его можно представить в виде конечного ряда (т.е. линейной комбинации) дискретных синусоид:

2π k (n + ϕ k )

x = ∑ C k cos

(ряд Фурье)

k = 0

Эквивалентная запись (каждый косинус раскладываем на синус и косинус, но теперь – без фазы):

2 π kn

2 π kn

x = ∑ A k cos

+ ∑ B k sin

(ряд Фурье)

k = 0

k = 0

Рис. 6 . Базисные функции ряда Фурье для 8-точеченого дискретного сигнала. Слева – косинусы, справа – синусы. Частоты увеличиваются сверху вниз.

Базисные синусоиды имеют кратные частоты. Первый член ряда (k =0) – это константа, называемаяпостоянной составляющей (DC offset ) сигнала. Самая первая синусоида (k =1) имеет такую частоту, что ее период совпадает с периодом самого исходного сигнала. Самая высокочастотная составляющая (k =N /2) имеет такую частоту, что ее период равен двум отсчетам. КоэффициентыA k и

B k называютсяспектром сигнала (spectrum ). Они показывают амплитуды си-

нусоид, из которых состоит сигнал. Шаг по частоте между двумя соседними синусоидами из разложения Фурье называется частотным разрешением спектра.

На рис. 6 показаны синусоиды, по которым происходит разложение дискретного сигнала из 8 точек. Каждая из синусоид состоит из 8 точек, то есть является обычным дискретным сигналом. Непрерывные синусоиды показаны на рисунке для наглядности.

вить исходный сигнал, вычислив сумму ряда Фурье в каждой точке. Разложение сигнала на синусоиды (т.е. получение коэффициентов) называется прямым преобразованием Фурье . Обратный процесс – синтез сигнала по синусоидам – называетсяобратным преобразованием Фурье (inverse Fourier transform ).

Алгоритм обратного преобразования Фурье очевиден (он содержится в формуле ряда Фурье; для проведения синтеза нужно просто подставить в нее коэффициенты). Рассмотрим алгоритм прямого преобразования Фурье, т.е. нахождения коэффициентов A k иB k .

2 π kn

2 π kn

от аргумента n является ор-

Система функций

K = 0,...,

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

Итак, коэффициенты A k иB k вычисляются как скалярные произведения (в не-

прерывном случае – интегралы от произведения функций, в дискретном случае

– суммы от произведения дискретных сигналов):

N − 1

2 π ki , приk = 1,...,

A k=

∑ x cos

−1

N i = 0

N − 1

A k=

∑ x cos2 π ki , приk = 0,

N i = 0

N − 1

2 π ki

NB 0 иB N 2 всегда равны нулю (т.к. соответствующие им «базисные»

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

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

Вычисление преобразований Фурье требует очень большого числа умножений (около N 2 ) и вычислений синусов. Существует способ выполнить эти преобразования значительно быстрее: примерно заN log2 N операций умножения.

Этот способ называется быстрым преобразованием Фурье(БПФ, FFT, fast Fourier transform). Он основан на том, что среди множителей (синусов) есть много повторяющихся значений (в силу периодичности синуса). Алгоритм БПФ группирует слагаемые с одинаковыми множителями, значительно сокращая число умножений. В результате быстродействие БПФ может в сотни раз превосходить быстродействие стандартного алгоритма (в зависимости от N). При этом следует подчеркнуть, что алгоритм БПФ является точным. Он даже точнее стандартного, т.к. сокращая число операций, он приводит к меньшим ошибкам округления.

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

N называется размеромили длиной БПФ(FFT size).

Комплексное ДПФ

До сих пор мы рассматривали ДПФ от действительных сигналов. Обобщим теперь ДПФ на случай комплексных сигналов. Пусть x , n =0,…,N -1 – исходный комплексный сигнал, состоящий изN комплексных чисел. ОбозначимX , k =0,…N -1 – его комплексный спектр, также состоящий изN комплексных чисел. Тогда справедливы следующие формулы прямого и обратного преобразо-

ваний Фурье (здесь j = − 1):

N − 1

X [ k] = ∑ x[ n] e− jnk (2 π N )

n= 0

N − 1

∑ X [ k ] e jnk(2 π N)

N k = 0

Если по этим формулам разложить в спектр действительный сигнал, то первые N /2+1 комплексных коэффициентов спектра будут совпадать со спектром «обычного» действительного ДПФ, представленным в «комплексном» виде, а остальные коэффициенты будут их симметричным отражением относительно



Понравилась статья? Поделитесь с друзьями!