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

         }

    }

}

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

Внешний цикл - это основные итерации. На каждой из них 2Nmax/N ДПФ (длиной по N/2 элементов каждое) преобразуются в Nmax/N ДПФ  (длиной по N элементов каждое).

Следующий цикл по k представляет собой цикл синхронного  вычисления элементов с индексами k и k + N/2 во всех результирующих ДПФ.

Самый внутренний цикл перебирает Nmax/N штук ДПФ одно за другим.

Именно так, а  не иначе: сначала вычисляются элементы с номерами 0 и N/2, во всех ДПФ в  количестве Nmax/N штук, потом с номерами 1 и N/2 + 1 и так далее. На рис.4 показана последовательность вычислений для N = 8. Такая последовательность обеспечивает однократное вычисление .

Обратите внимание, что x[mpNd2] вычисляется раньше, чем x[k], чтобы значение x[k] не было модифицировано прежде, чем будет использовано.

Алгоритм обратного дискретного преобразования Фурье отличается только тем, что вместо -j надо использовать +j и после всех вычислений поделить каждое x[k] на Nmax.

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

(x,y)(x',y')=(xx'-yy')(xy'+x'y)  
(x,y)+(x',y')=(x+x')(y+y')  
e^-jv=(cos v, - sinv)

где v - действительное число.

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

.    (17)

Это позволило  бы избавиться от возведений в степень и обойтись одними умножениями, если заранее вычислить W_N для N = 2, 4, 8,…, Nmax. Но то, что можно делать в математике, далеко не всегда можно делать в программах. По мере увеличения k поворачивающий множитель будет изменяться, но вместе с тем будет расти погрешность. Ведь вычисления с плавающей точкой на реальном компьютере совсем без погрешностей невозможны. И после N/2 подряд умножений в поворачивающем множителе может накопиться огромное отклонение от точного значения. Вспомним правило: при умножении величин их относительные погрешности складываются. Так что, если погрешность одного умношения равна s%, то после 1000 умножений она может достигнуть в худшем случае 1000s%.

Этого можно  было бы избежать, будь число W_N целым, но оно не целое при N > 2, так как вычисляется через синус и косинус:

Как же избежать множества обращений к весьма медленным функциям синуса и косинуса? Здесь нам приходит на помощь давно  известный алгоритм вычисления степени  через многократные возведения в квадрат и умножения. Например:

В данном случае нам понадобилось всего 5 умножений (учитывая, что   не нужно вычислять дважды) вместо 13. В худшем случае для возведения в степень от 1 до N/2-1 нужно log₂N умножений вместо N/2, что дает вполне приемлемую погрешность для большинства практических задач.

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

,

для хранения которых  нужно дополнительно log₂(Nmax) комплексных ячеек памяти:

Если очередное  не было вычислено предварительно, то берется двоичное представление k и анализируется. Каждому единичному биту соответствует ровно один множитель. В общем случае единице в бите с номером b (младший бит имеет номер 0) соответствует множитель , который хранится в b-й ячейке упомянутого выше массива.

Есть способ уменьшить количество умножений  для вычисления  до одного на два цикла. Но для этого нужно отвести N/2 комплексных ячеек для хранения всех предыдущих . Алгоритм достаточно прост. Нечетные элементы вычисляются по формуле (17). Четные вычисляются по формуле:

То есть, ничего не вычисляется, а берется одно из значений, вычисленных на предыдущем шаге, когда N было вдвое меньше. Чтобы  не нужно было копировать величины на новые позиции достаточно их сразу располагать в той позиции, которую они займут при N = Nmax и вводить простую поправку Skew (см. листинг программы).

И последние  пояснения относительно листинга.

Во-первых, здесь  присутствует реализация простейших операций над комплексными числами (классы Complex и ShortComplex), оптимизированная под данную задачу. Обычно та или иная реализация уже есть для многих компиляторов, так что можно использовать и ее.

Во-вторых, массив W2n содержит заранее вычисленные коэффициенты W₂, W₄, W₈,...,W_Nmax.

В-третьих, для  вычислений используются наиболее точное представление чисел с плавающей  точкой в C++: long double размером в 10 байт на платформе Intel. Для хранения результатов  в массивах используется тип double длиной 8 байт. Причина - не в желании сэкономить 20% памяти, а в том, что 64-битные и 32-битные процессоры лучше работают с выровненными по границе 8 байт данными, чем с выровненными по границе 10 байт.

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

Пусть у нас  есть функция синуса x = sin(t).

Максимальная  амплитуда этого колебания равна 1. Если умножить его на некоторый  коэффициент A, то получим тот же график, растянутый по вертикали в A раз: x = Asin(t).

Период колебания  равен 2π. Если мы хотим увеличить  период до T, то надо умножить переменную t на коэффициент. Это вызовет растяжение графика по горизонтали: x = A sin(2πt/T).

Частота колебания  обратна периоду: ν = 1/T. Также говорят о круговой частоте, которая вычисляется по формуле: ω= 2πν = 2πT. Откуда: x = A sin(ωt).

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

Очень похоже выглядит и выражение гармоники через  косинус:

Большой разницы  нет. Достаточно изменить фазу на π/2, чтобы  перейти от синуса к косинусу и  обратно. Далее будем подразумевать  под гармоникой функцию косинуса:

x = A cos(2πt/T + φ) = A cos(2πνt + φ) = A cos(ωt + φ)     (18)

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

Преобразуем (18) по формуле косинуса суммы:

x = A cos φ cos(2πt/T) - A sin φ sin(2πt/T)     (19)

Выделим в (19) элементы, независимые от t, и обозначим  их как Re и Im:

x = Re cos(2πt/T) - Im sin(2πt / T)    (20)

Re = A cos φ, Im = A sin φ

По величинам Re и Im можно однозначно восстановить амплитуду и фазу исходной гармоники:

   и        (21)

Рассмотрим очень  распространенную практическую ситуацию. Пусть у нас есть звуковое или  какое-то иное колебание в виде функции x = f(t). Пусть это колебание было записано в виде графика для отрезка времени [0, T]. Для обработки компьютером нужно выполнить дискретизацию. Отрезок делится на N-1 равных частей, границы частей обозначим t_n = nT/N. Сохраняются N значений функции на границах частей: x_n = f(t_n) = { x₀, x₁, x₂,..., x_N }.

     

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

    (22)

Теперь возьмем  обратное преобразование Фурье:

    (23)

Выполним над  этой формулой следующие действия: разложим каждое комплексное X_k на мнимую и действительную составляющие X_k = Re_k + j Im_k; разложим экспоненту по формуле Эйлера на синус и косинус действительного аргумента; перемножим; внесем 1/N под знак суммы и перегруппируем элементы в две суммы:

    (24)

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

    (25)

Поскольку при дискретизации мы брали t_n = nT/N и , то можем выполнить замену: n = t_nN/T. Следовательно, в синусе и косинусе вместо 2πkn/N можно написать 2πkt_n/T. В результате получим:

    (26)

Сопоставим эту  формулу с формулой (20) для гармоники:

x = Re cos(2πt/T) - Im sin(2πt / T)    (20)

Слагаемые суммы (26) аналогичны формуле (20), а формула (20) описывает гармоническое колебание. Значит сумма (26) представляет собой  сумму из N гармонических колебаний  разной частоты, фазы и амплитуды.

Выше объяснялось, каким образом формула вида (20) может быть преобразована в формулу вида (18):

x = A cos(2πt/T + φ)    (18)

Выполним такое  же преобразование для слагаемых  суммы (26), преобразуем их из вида (20) в вид (18). Получим:

    (27)

Далее будем  функцию 

G_k(t) = A_k cos(2πtk/T + φ_k)     (28)

называть  k-й гармоникой.

Для вычисления A_k и φ_k надо использовать формулу (21). Теперь выпишем в одном месте все формулы, которые связывают амплитуду, фазу, частоту и период каждой из гармоник с коэффициентами X_k:

    (29)

Итак.

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

Что такое зеркальный эффект.

Предположим, что  исходный сигнал состоял из суммы  гармоник. f_s(t) = A_s cos(2πtm_s / T + φ_s). Пусть мы этот сигнал подвергли дискретизации, выполнили над ним прямое и обратное преобразование Фурье. Представили в виде суммы гармоник G_k(t) = A_k cos(2πtk / T + φ_k), как это описано в предыдущей главе. Спрашивается, эти гармоники G_k - те же самые, что и исходные гармоники f_s или нет? Оказывается, нет, не те. Но кое-какую информацию об исходных гармониках все же можно попытаться восстановить.

Эта задача имеет  практический интерес. Пусть нам  дан некий сигнал, который физически  получился как сумма гармонических  колебаний (или близких к ним). Простейший пример: кто-то сыграл аккорд, аккорд записан как звуковое колебание  в виде mp3 или wav-файла; и надо восстановить, из каких нот аккорд состоял. Или другой случай. Во время испытаний самолета возик флаттер (разрушительные нарастающие колебания), самолет разбился, но самописцы в черном ящике записали изменение перегрузки (суммарное механическое колебание). Надо посмотреть, из каких гармоник состояло это колебание. Каждой гармонике соответствует некоторая часть конструкции. В результате можно понять, какие части самолета колебались сильнее всего и вызвали флаттер.