Разработка анализатора спектра ультразвукового сигнала

n, к – диапазон последовательности;

N – количество отсчетов спектра;

К = 0, 1, …, N -1.

Обратное преобразование в таком случае будет иметь вид, как описано в формуле 20.

,                  (20)
 

При внимательном рассмотрении можно заметить, что индекс при Hn  принимает N+1 значение, в то время как при hk  - только N значений. Таким образом, как будто бы получается, что функция H (f) содержит в себе больше информации, чем h (t). На самом деле это не так, поскольку значения H-N/2  и HN/2  совпадают.

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

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

Наиболее популярным из алгоритмов ускоренного вычисления ДПФ является  метод Cooley-Tukey, позволяющий вычислить ДПФ для числа отсчетов N = 2 k за время порядка Nlog2 N (отсюда и название - быстрое преобразование Фурье, БПФ). Этот способ чем-то  напоминает быструю сортировку. В ходе работы алгоритма также проводится рекурсивное разбиение массива чисел на два подмассива и сведение вычисления ДПФ от целого массива к вычислению ДПФ от подмассивов в отдельности.

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

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

Широко распространено ошибочное мнение о том, что метод Cooley-Tukey - единственный существующий метод выполнения БПФ, а само БПФ существует только для случая N равно 2 k. На самом деле это не так - существуют алгоритмы БПФ для любого числа отсчетов. Метод Винограда позволяет решить задачу для простого числа отсчетов N. Этот же алгоритм может быть легко обобщен на случай, когда N является степенью произвольного простого числа, а также на случай, когда число N является произведением степеней простых чисел, т.е. N является произвольным числом, чье разложение на простые множители нам известно.

Существуют алгоритмы БПФ для произвольного числа отсчетов, но наиболее широкое распространение получил только алгоритм для случая N равно 2 k, что является существенным ограничением.

Алгоритм, построенный по методу Cooley-Tukey, обладает рядом очень хороших технологических свойств. Структура алгоритма и его базовые операции не зависят от числа отсчетов. Алгоритм легко распараллеливается с использованием базовой операции и конвейеризуется, а также легко каскадируется (коэффициенты БПФ для 2N отсчетов могут быть легко получены преобразованием коэффициентов двух БПФ по N отсчетов, полученных "прореживанием" через один исходных 2N отсчетов). Алгоритм прост и компактен, не требует дополнительной оперативной памяти и допускает обработку данных "на месте". Существует целый ряд оптимизированных, именно для этого алгоритма, DSP-процессоров.

4.3.1 БПФ комплексной функции (прямое и обратное)

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

4.3.1.1Основные принципы

Прямое и обратное дискретные преобразования Фурье задаются формулами 21, 22.

,                            (21)

,                  (22)

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

4.3.2 Улучшения алгоритма быстрого преобразования Фурье

Улучшения носят технический характер. Первое из них позволяет избавиться от рекурсии и повысить скорость, за счет перестановки элементов массива с тем, чтобы было можно вычислять преобразование в цикле. Этот процесс называется bit - reversing, поскольку после перестановки новый номер элемента является записанным, наоборот, в двоичной системе исчисления старым номером. Причины этого можно понять, если попробовать в ходе рекурсивного разбиения переупорядочивать элементы массива так, чтобы сначала шли элементы, соответствующие слагаемым "левой" суммы, а потом - слагаемые "правой". При такой перестановке становится возможным проводить преобразование "на месте", без выделения дополнительного массива.

Улучшение касается использования тригонометрических функций. Алгоритм вызывает относительно медленную функцию sin лишь log2 N раз для инициализации рекуррентной последовательности, а все остальные требуемые значения получаются по специальной рекуррентной формуле, обладающей особой устойчивостью к вычислительным погрешностям.

Описание

На входе процедура получает параметры:

-nn - число значений функции. Должно быть степенью двойки.

-a - массив из 2*nn вещественных чисел, нумерация элементов от 0 до 2*nn-1. Если проводим прямое преобразование, этот массив задает nn значений функции.

Если преобразование обратное, то здесь хранятся соответствующие частоты.

-InverseFFT - направление преобразования. Если проводим прямое преобразование, то устанавливаем параметр в False. Устнавливаем его в True, если преобразование обратное.

На выходе:

-a - содержит результат преобразования.

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

Рассмотрим более подробно формат массива комплексных значений. Пусть у нас есть N значений функции, расположенных с интервалом Delta. Значения хранятся парами вещественных чисел, где первое число в паре соответствует действительной части, а второе - мнимой, как указано на рисунке 5. Каждому значению hk  соответствует определенный момент времени tk .

Рисунок 5 - формат массива комплексных значений

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

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

4.3.3 БПФ вещественной функции (прямое и обратное)

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

Описание

На входе процедура получает параметры:

-tnn - число значений функции. Должно быть степенью двойки;

-a - массив из tnn вещественных чисел, нумерация элементов от 0 до tnn-1.

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

-InverseFFT - направление преобразования.

Если проводим прямое преобразование, то устанавливаем параметр в False. Устанавливаем его в True, если преобразование обратное.

На выходе:

-a - содержит результат преобразования.

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

Рассмотрим более подробно формат массива. Пусть есть N значений функции, расположенных с интервалом Delta, показанные на рисунке 6. Каждому значению hk  соответствует определенный момент времени tk .

Рисунок 6 - формат массива вещественных значений

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

Если сравнить рисунок 5 и рисунок 6, то заметим, что частоты с f-1  по f-N/2+1  исчезли, а от частот f0  и fN/2  остались только вещественные части, занявшие прежнее место комплексной частоты f0 . Причиной этого являются свойства симметрии преобразования Фурье:

для вещественной функции h(t) приведена формула 23.

              H(-f) = H (f),              (23)

Таким образом, частоты с f-1  по f-N/2+1  уже не несут в себе новой информации, т.к. получаются комплексным сопряжением своих симметричных близнецов, а у частот f0  и fN/2  мнимые части равны нолю.

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

4.3.4 Быстрая свертка

Свертка функций s(t) и r(t) определяется формулой 24.

                                      ,                            (24)

где s(t), r(t) – функции.

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

                                            ,                               (25)                                                                                                                                                                                     

где  -N, P – диапазон, за пределами которого r(t) = 0.

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

Обычно функция отклика принимает ненулевые значения только на некотором отрезке [-N, P], включающем в себя ноль, как показано на

рисунке 7

Рисунок 7 - функция отклика

В то время, как сигнал тоже принимает ненулевые значения только на некотором отрезке [A, B], но положение этого отрезка не имеет жесткой привязки к началу координат. В таких случаях интересно значение свертки не на всей действительной оси, а только в тех точках, в которых измерили сигнал.

Если длина отклика мала (меньше, чем двоичный логарифм длины сигнала), можно провести корелляцию по определению, но в тех случаях, когда длина отклика сравнима с длиной сигнала, этот процесс становится слишком трудоемким. В таких случаях на помощь приходит быстрое преобразование Фурье. Дело в том, что при свертке функций происходит покомпонентное перемножение их Фурье - образов. Это свойство позволяет провести БПФ (за время порядка N·log2 N), покомпонентное перемножить полученные наборы комплексных чисел, затем провести обратное преобразование - и искомая свертка получена. Единственная тонкость в работе алгоритма связана с тем, что в случае дискретного преобразования Фурье (в отличие от непрерывного) происходит свертка двух периодических функций, т.е. наши наборы значений задают именно периоды этих функций, а не просто значения на каком-то отдельном участке оси. Т.е. алгоритм считает, что за точкой xN  идет не ноль, а точка x0 , затем x1  и так далее по кругу. Поэтому, чтобы свертка корректно считалась, алгоритм приписывает к сигналу достаточно длинную последовательность нулей.

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

Параметры процедуры

На входе:

-Signal - сигнал, с которым проводим свертку. Массив вещественных чисел, нумерация элементов от 0 до SignalLen-1;

-SignalLen - длина сигнала;

-Responce - функция отклика.

Нумерация элементов от 0 до NegativeLen+PositiveLen хранит набор значений, начиная с r(-NegativeLen) и заканчивая r(PositiveLen) (в порядке возрастания аргументов).

-NegativeLen - "отрицательная длина" отклика;

-PositiveLen - "положительная длина" отклика; За пределами

[-NegativeLen, PositiveLen] отклик равен нолю.

На выходе:

-Signal - содержит значения свертки в тех же точках, что были переданы в функцию.