Материал: Дискретное преобразование Фурье
Введение
Дискретное преобразование Фурье (в англоязычной литературе DFT, Discrete Fourier Transform) -- это одно из преобразований Фурье, широко применяемых в алгоритмах цифровой обработки сигналов (его модификации применяются в сжатии звука в MP3, сжатии изображений в JPEG и др.), а также в других областях, связанных с анализом частот в дискретном (к примеру, оцифрованном аналоговом) сигнале. Дискретное преобразование Фурье требует в качестве входа дискретную функцию. Такие функции часто создаются путём дискретизации (выборки значений из непрерывных функций). Дискретные преобразования Фурье помогают решать частные дифференциальные уравнения и выполнять такие операции, как свёртки. Дискретные преобразования Фурье также активно используются в статистике, при анализе временных рядов. Существуют многомерные дискретные преобразования Фурье.
Преобразование Фурье (Fouriertransform)- это разложение функций на синусоиды (далее косинусные функции мы тоже называем синусоидами, т.к. они отличаются от «настоящих» синусоид только фазой). Существует несколько видов преобразования Фурье.
- 1. Непериодический непрерывный сигнал можно разложить в интеграл Фурье.
- 2. Периодический непрерывный сигнал можно разложить в бесконечный ряд Фурье.
- 3. Непериодический дискретный сигнал можно разложить в интеграл Фурье.
- 4. Периодический дискретный сигнал можно разложить в конечный ряд Фурье.
Компьютер способен работать только с ограниченным объемом данных, следовательно, реально он способен вычислять только последний вид преобразования Фурье. Рассмотрим его подробнее.
ДПФ вещественного сигнала
Пусть дискретный сигнал x[n] имеет период N точек. В этом случае его можно представить в виде конечного ряда (т.е. линейной комбинации) дискретных синусоид:
Эквивалентная запись (каждый косинус раскладываем на синус и косинус, но теперь без фазы):
Рис. 1 - Базисные функции ряда Фурье для 8-точеченого дискретного сигнала. Слева - косинусы, справа - синусы. Частоты увеличиваются сверху вниз
Базисные синусоиды имеют кратные частоты. Первый член ряда (k=0) - это константа, называемая постоянной составляющей (DC offset) сигнала. Самая первая синусоида (k=1) имеет такую частоту, что ее период совпадает с перио-дом самого исходного сигнала. Самая высокочастотная составляющая (k=N/2) имеет такую частоту, что ее период равен двум отсчетам. Коэффициенты и называются спектром сигнала (spectrum). Они показывают амплитуды синусоид, из которых состоит сигнал. Шаг по частоте между двумя соседними синусоидами из разложения Фурье называется частотным разрешением спектра.
На рис. 1 показаны синусоиды, по которым происходит разложение дискретного сигнала из 8 точек. Каждая из синусоид состоит из 8 точек, то есть является обычным дискретным сигналом. Непрерывные синусоиды показаны на рисунке для наглядности.
Как мы увидим далее, для каждого сигнала можно однозначно определить коэффициенты и .
Зная эти коэффициенты, можно однозначно восстановить исходный сигнал, вычислив сумму ряда Фурье в каждой точке. Разложение сигнала на синусоиды (т.е. получение коэффициентов) называется прямым преобразованием Фурье. Обратный процесс - синтез сигнала по синусоидам - называется обратным преобразованием Фурье (inverseFouriertransform).
Алгоритм обратного преобразования Фурье очевиден (он содержится в формуле ряда Фурье; для проведения синтеза нужно просто подставить в нее коэффициенты). Рассмотрим алгоритм прямого преобразования Фурье, т.е. нахождения коэффициентов и .
Система функций:
от аргумента n является ортогональным базисом в пространстве периодических дискретных сигналов с периодом N. Это значит, что для разложения по ней любого элемента пространства (сигнала) нужно посчитать скалярные произведения этого элемента со всеми функциями системы, и полученные коэффициенты нормировать. Тогда для исходного сигнала будет справедлива формула разложения по базису с коэффициентами и .
Итак, коэффициенты и вычисляются как скалярные произведения (в непрерывном случае - интегралы от произведения функций, в дискретном случае - суммы от произведения дискретных сигналов):
Возникает вопрос: почему в исходном сигнале N чисел, а описывается он с помощью N+2 коэффициентов? Вопрос разрешается следующим образом: коэффициенты и всегда равны нулю (т.к. соответствующие им «базисные» сигналы тождественно равны нулю в дискретных точках), и их можно отбросить при вычислении обратного и прямого преобразований Фурье.
Итак, мы выяснили, что спектральное представление сигнала полностью эквивалентно самому сигналу. Между ними можно перемещаться, используя прямое и обратное преобразования Фурье. Алгоритм вычисления этих преобразований содержится в приведенных формулах.
Комплексное ДПФ
До сих пор мы рассматривали ДПФ от действительных сигналов. Обобщим теперь ДПФ на случай комплексных сигналов. Пусть x[n], n=0,…,N-1 - исходный комплексный сигнал, состоящий из N комплексных чисел. Обозначим X[k], k=0,…N-1 - его комплексный спектр, также состоящий из N комплексных чисел. Тогда справедливы следующие формулы прямого и обратного преобразований Фурье (здесь 1?=j):
Если по этим формулам разложить в спектр действительный сигнал, то первые N/2+1 комплексных коэффициентов спектра будут совпадать со спектром «обычного» действительного ДПФ, представленным в «комплексном» виде, а остальные коэффициенты будут их симметричным отражением относительно половины частоты дискретизации. Для косинусных коэффициентов отражение четное, а для синусных - нечетное.
Двумерное ДПФ
Для изображений, представляющих собой двумерный сигнал, спектром является также двумерный сигнал. Базисные функции преобразования Фурье имеют вид
причем фазы также могут быть различны. На изображении каждая из этих базисных функций представляют собой волну определенной частоты, определенной ориентации и определенной фазы.
Здесь N1xN2 - размер исходного сигнала, он же - размер спектра. k1 и k2 - это номера базисных функций (номера коэффициентов двумерного ДПФ, при которых эти функции находятся). Поскольку размер спектра равен размеру исходного сигнала, то k1 = 0,…,N1-1; k2 = 0,…,N2-1.
n1 и n2 - переменные-аргументы базисных функций. Поскольку область определения базисных функций совпадает с областью определения сигнала, то n1 = 0,…,N1-1; n2 = 0,…,N2-1.
Двумерное ДПФ (в комплексной форме) определяется следующими формулами (здесь x[n1,n2] - исходный сигнал, а X[k1,k2] - его спектр):
Непосредственное вычисление двумерного ДПФ по приведенным формулам требует огромных вычислительных затрат. Однако можно доказать, что двумерное ДПФ обладает свойством сеперабельности, т.е. его можно вычислить последовательно по двум измерениям. Для вычисления двумерного ДПФ достаточно вычислить одномерные комплексные ДПФ всех строк изображения, а затем вычислить в результирующем «изображении» одномерные комплексные ДПФ всех столбцов. При этом результаты всех одномерных комплексных ДПФ нужно записывать на место исходных данных для этих ДПФ. Например, при вычислении одномерного ДПФ первой строки изображения нужно результат ДПФ записать в первую строку этого изображения (он имеет тот же размер). Для этого нужно каждый «пиксель» хранить в виде комплексного числа.
Ряды Фурье
Разложению в ряды Фурье подвергаются периодические сигналы. Периодическую функцию любой формы, заданную на интервале одного периода Т = b-a и удовлетворяющую на этом интервале условиям Дирехле (ограниченная, кусочно-непрерывная, с конечным числом разрывов 1-го рода), можно представить в виде ряда Фурье:
s(t) =Snexp(jnDwt), Sn = S(nDw), Dw = 2p/T, (1)
где весовые коэффициенты Sn ряда определяются по формуле:
Sn = (1/T)s(t) exp(-jnDwt) dt.(2)
Ряд Фурье представляет собой ансамбль комплексных экспонент exp(jnDwt) с частотами, образующими арифметическую прогрессию. Функцию весовых коэффициентов S(nDw) принято называть комплексным спектром периодического сигнала или фурье-образом функции s(t). Спектр периодического сигнала является дискретной функцией, т.к. он определен только для целых значений n с шагом по частоте, обратным периоду: Dw = 2p/Т (или Df = 1/T). Первую частотную составляющую спектра при n = 1, равную w1 = 1Dw = 2p/T (или f1 = 1/T), называют основной частотой сигнала (первой гармоникой), остальные частоты дискретного спектра nDw при n>1 называют гармониками сигнала. Значения S(nDw) по положительным и отрицательным значениям n являются комплексно сопряженными. Шаг по частоте Dw между двумя соседними синусоидами называется частотным разрешением спектра.
С чисто математических позиций множество функций exp(jnDwt), -< n < образует бесконечномерный базис линейного пространства L2[a,b] ортогональных синус-косинусных функций, а коэффициенты Sn по (2) представляют собой проекции сигнала s(t) на эти базисные функции. Соответственно, сигнал s(t) в форме ряда Фурье (1) - это бесконечномерный вектор в пространстве L2[a,b], точка с координатами Sn по базисным осям пространства exp(jnDwt).
Коэффициенты Sn в (2) отображают функцию s(t) в новое пространство единственным образом. Если функция s(t) непрерывна, то ряд (1) сходится равномерно к s(t), при этом ошибка аппроксимации ||s(t)-sN(t)|| функции s(t) с усечением ряда (1) до ±N членов меньше ошибки аппроксимации любым другим рядом с тем же количеством членов. Если s(t) не является непрерывной (имеет разрывы), но конечна по энергии (квадратично интегрируема), то метрика ||s(t)-sN(t)|| стремится к нулю при N > ?, при этом в точках разрыва сумма ряда стремится к (s(t+)+s(t-))/2.
Подынтегральную функцию экспоненты в выражении (2) с использованием тождества Эйлера
exp(±jwt) = cos(wt) ± jsin(wt)
можно разложить на косинусную и синусную составляющие и выразить комплексный спектр в виде действительной и мнимой части:
Sn = (1/T)s(t) [cos(nDwt) - j sin(nDwt)] dt = Аn - jBn. (3)
An ? A(nDw)= (1/T)s(t) cos(nDwt) dt, (4)
Bn ? B(nDw)= (1/T) s(t) sin(nDwt) dt.(5)
На рис. 1 приведен пример периодического сигнала (прямоугольный импульс на интервале (1-3.3), повторяющийся с периодом Т=40) и форма действительной и мнимой части его спектра. Обратим внимание, что действительная часть спектра является четной относительно нуля функцией A(nDw) = A(-nDw), так как при вычислении значений A(nDw) по формуле (4) используется четная косинусная функция cos(nDwt) = cos(-nDwt). Мнимая часть спектра является нечетной функцией B(nDw) = -B(-nDw), так как для ее вычисления по (5) используется нечетная синусная функция sin(nDwt) = - sin(-nDwt).
Рис. 2 - Сигнал и его комплексный спектр
Комплексные числа дискретной функции (3) могут быть представлены в виде модулей и аргументов комплексной экспоненты, что дает следующую форму записи комплексного спектра:
Sn = Rnexp(jjn), (3')
Rn2 ? R2(nDw) = A2(nDw)+B2(nDw),
jn ? j(nDw) = arctg(-B(nDw)/A(nDw)).
Модуль спектра R(nDw) называют двусторонним спектром амплитуд или АЧХ - амплитудно-частотной характеристикой сигнала, а аргумент спектра (последовательность фазовых углов j(nDw)) - двусторонним спектром фаз или ФЧХ - фазочастотной характеристикой. Спектр амплитуд всегда представляет собой четную функцию: R(nDw) = R(-nDw), а спектр фаз нечетную: j(nDw) = -j(-nDw). Пример спектра в амплитудном и фазовом представлении для сигнала, показанного на рис.1, приведен на рис2. При рассмотрении спектра фаз следует учитывать периодичность 2p угловой частоты (при уменьшении фазового значения до величины менее -p происходит сброс значения -2p).
Рис. 3 - Модуль и аргумент спектра
Рис. 4 - Ортогональность функций
Если функция s(t) является четной, то все значения B(nDw) по (5) равны нулю, т.к. четные функции ортогональны синусным гармоникам и подынтегральное произведение s(t)·sin(nDwt) дает нулевой интеграл. Следовательно, спектр функции будет представлен только вещественными коэффициентами. Напротив, при нечетности функции s(t) обнуляются все значения коэффициентов А(nDw) (нечетные функции ортогональным косинусным гармоникам) и спектр является чисто мнимым. Этот фактор не зависит от выбора границ задания периода функции на числовой оси. На рис.3(А) можно наглядно видеть ортогональность первой гармоники синуса и четной функции, а на рис.3(В) соответственно косинуса и нечетной функции в пределах одного периода. Учитывая кратность частот последующих гармоник первой гармонике спектра, ортогональность сохраняется для всех гармоник ряда Фурье.