Описание Image Processing Toolbox. Преобразование Фурье
Линейная фильтрация изображений может осуществляться как в пространственной, так и в частотной области. При этом считается, что "низким" пространственным частотам соответствует основное содержание изображения - фон и крупноразмерные объекты, а "высоким" пространственным частотам - мелкоразмерные объекты, мелкие детали крупных форм и шумовая компонента.
Традиционно для перехода в область пространственных частот используются методы, основанные на $\textit{преобразовании Фурье}$. В последние годы все большее применение находят также методы, основанные на $\textit{вейвлет-преобразовании (wavelet-transform)}$.
Преобразование Фурье.
Преобразование Фурье позволяет представить практически любую функцию или набор данных в виде комбинации таких тригонометрических функций, как синус и косинус, что позволяет выявить периодические компоненты в данных и оценить их вклад в структуру исходных данных или форму функции. Традиционно различаются три основные формы преобразования Фурье: интегральное преобразование Фурье, ряды Фурье и дискретное преобразование Фурье.
Интегральное преобразование Фурье переводит вещественную функцию в пару вещественных функций или одну комплексную функцию в другую.
Вещественную функцию $f(x)$ можно разложить по ортогональной системе тригонометрических функций, то есть представить в виде
$$ f\left(x \right)=\int\limits_0^\infty {A\left(\omega \right)} \cos \left({2\pi \omega x} \right)d\omega -\int\limits_0^\infty {B\left(\omega \right)} \sin \left({2\pi \omega x} \right)d\omega , $$
где $A(\omega)$ и $B(\omega)$ называются интегральными косинус- и синус-преобразованиями:
$$ A\left(\omega \right)=2\int\limits_{-\infty }^{+\infty } {f\left(x \right)} \cos \left({2\pi \omega x} \right)dx; \quad B\left(\omega \right)=2\int\limits_{-\infty }^{+\infty } {f\left(x \right)} \sin \left({2\pi \omega x} \right)dx. $$
Ряд Фурье представляет периодическую функцию $f(x)$, заданную на интервале $$, в виде бесконечного ряда по синусам и косинусам. То есть периодической функции $f(x)$ ставится в соответствие бесконечная последовательность коэффициентов Фурье
$$ f\left(x \right)=\frac{A_0 }{2}+\sum\limits_{n=1}^\infty {A_n } \cos \left({\frac{2\pi xn}{b-a}} \right)+\sum\limits_{n=1}^\infty {B_n \sin \left({\frac{2\pi xn}{b-a}} \right)} , $$
$$ A_n =\frac{2}{b-a}\int\limits_a^b {f\left(x \right)} \cos \left({\frac{2\pi nx}{b-a}} \right)dx; \quad B_n =\frac{2}{b-a}\int\limits_a^b {f\left(x \right)} \sin \left({\frac{2\pi nx}{b-a}} \right)dx. $$
Дискретное преобразование Фурье переводит конечную последовательность вещественных чисел в конечную последовательность коэффициентов Фурье.
Пусть $\left\{ {x_i } \right\}, i= 0,\ldots, N-1 $ - последовательность вещественных чисел - например, отсчеты яркости пикселов по строке изображения. Эту последовательность можно представить в виде комбинации конечных сумм вида
$$ x_i =a_0 +\sum\limits_{n=1}^{N/2} {a_n } \cos \left({\frac{2\pi ni}{N}} \right)+\sum\limits_{n=1}^{N/2} {b_n \sin \left({\frac{2\pi ni}{N}} \right)} , $$
$$ a_0 =\frac{1}{N}\sum\limits_{i=0}^{N-1} {x_i } , \quad a_{N/2} =\frac{1}{N}\sum\limits_{i=0}^{N-1} {x_i } \left({-1} \right)^i, \quad a_k =\frac{2}{N}\sum\limits_{i=0}^{N-1} {x_i \cos \left({\frac{2\pi ik}{N}} \right)}, $$
$$
b_k =\frac{2}{N}\sum\limits_{i=0}^{N-1} {x_i \sin \left({\frac{2\pi ik}{N}} \right)},
\quad
i\le k Основное отличие между тремя формами преобразования Фурье заключается в том,
что если интегральное преобразование Фурье определено по всей области
определения функции $f(x)$, то ряд и дискретное преобразование Фурье
определены только на дискретном множестве точек, бесконечном для ряда Фурье
и конечном для дискретного преобразования. Как видно из определений преобразования Фурье, наибольший интерес для систем
цифровой обработки сигналов представляет дискретное преобразование Фурье.
Данные, получаемые с цифровых носителей или источников информации,
представляют собой упорядоченные наборы чисел, записанные в виде векторов
или матриц. Обычно принимается, что входные данные для дискретного преобразования
представляют собой равномерную выборку с шагом $\Delta $, при этом величина
$T=N\Delta $ называется длиной записи, или основным периодом. Основная
частота равна $1/T$. Таким образом, в дискретном преобразовании Фурье
производится разложение входных данных по частотам, которые являются целым
кратным основной частоты. Максимальная частота, определяемая размерностью
входных данных, равна $1/2 \Delta $ и называется $\it{частотой Найквиста}$. Учет частоты Найквиста
имеет важное значение при использовании дискретного преобразования. Если
входные данные имеют периодические составляющие с частотами, превышающими
частоту Найквиста, то при вычислении дискретного преобразования Фурье
произойдет подмена высокочастотных данных более низкой частотой, что может
привести к ошибкам при интерпретации результатов дискретного преобразования. Важным инструментом анализа данных является также $\it{энергетический спектр}$. Мощность сигнала на
частоте $\omega $ определяется следующим образом: $$
P \left(\omega \right)=\frac{1}{2}\left({A \left(\omega \right)^2+B \left(\omega \right)^2} \right) .
$$ Эту величину часто называют $\it{энергией сигнала}$ на частоте $\omega $. Согласно теореме
Парсеваля общая энергия входного сигнала равна сумме энергий по всем
частотам. $$
E=\sum\limits_{i=0}^{N-1} {x_i^2 } =\sum\limits_{i=0}^{N/2} {P \left({\omega _i } \right)} .
$$ График зависимости мощности от частоты называется энергетическим спектром
или спектром мощности. Энергетический спектр позволяет выявлять скрытые
периодичности входных данных и оценивать вклад определенных частотных
компонент в структуру исходных данных. Кроме тригонометрической формы записи дискретного преобразования Фурье
широко используется $\it{комплексное представление}$. Комплексная форма записи преобразования Фурье широко
используется в многомерном анализе и в частности при обработке изображений. Переход из тригонометрической в комплексную форму осуществляется на
основании формулы Эйлера $$
e^{j\omega t}=\cos \omega t+j\sin \omega t,
\quad
j=\sqrt {-1} .
$$ Если входная последовательность представляет собой $N$ комплексных чисел, то
ее дискретное преобразование Фурье будет иметь вид $$
G_m =\frac{1}{N}\sum\limits_{n=1}^{N-1} {x_n } e^{\frac{-2\pi jmn}{N}},
$$ а обратное преобразование $$
x_m =\sum\limits_{n=1}^{N-1} {G_n } e^{\frac{2\pi jmn}{N}}.
$$ Если входная последовательность представляет собой массив вещественных
чисел, то для нее существует как комплексное, так и синусно-косинусное
дискретное преобразование. Взаимосвязь этих представлений выражается
следующим образом: $$
a_0 =G_0 ,
\quad
G_k =\left({a_k -jb_k } \right)/2,
\quad
1\le k\le N/2;
$$ остальные $N/2$ значений преобразования являются комплексно сопряженными и не
несут дополнительной информации. Поэтому график спектра мощности дискретного
преобразования Фурье симметричен относительно $N/2$. Простейший способ вычисления дискретного преобразования Фурье (ДПФ) -
прямое суммирование, оно приводит к $N$ операциям на каждый коэффициент. Всего
коэффициентов $N$, так что общая сложность $O\left({N^2} \right)$. Такой
подход не представляет практического интереса, так как существуют гораздо
более эффективные способы вычисления ДПФ, называемые быстрым преобразованием
Фурье (БПФ), имеющее сложность $O (N\log N)$.
БПФ применяется только к последовательностям, имеющим длину (число
элементов), кратную степени 2. Наиболее общий принцип, заложенный в алгоритм
БПФ, заключается в разбиении входной последовательности на две
последовательности половинной длины. Первая последовательность заполняется
данными с четными номерами, а вторая - с нечетными. Это дает возможность
вычисления коэффициентов ДПФ через два преобразования размерностью $N/2$. Обозначим $\omega _m =e^{\frac{2\pi j}{m}}$, тогда $G_m =\sum\limits_{n=1}^{(N/2)-1} {x_{2n} } \omega
_{N/2}^{mn} +\sum\limits_{n=1}^{(N/2)-1} {x_{2n+1} }
\omega _{N/2}^{mn} \omega _N^m $. Для $m < N/2$ тогда можно записать $G_m =G_{\textrm{even}} \left(m \right)+G_{\textrm{odd}}
\left(m \right)\omega _N^m $. Учитывая, что элементы ДПФ с индексом
б ольшим, чем $N/2$, являются комплексно сопряженными к элементам с индексами
меньшими $N/2$, можно записать $G_{m+(N/2)} =G_{\textrm{even}} \left(m \right)-G_{\textrm{odd}}
\left(m \right)\omega _N^m $. Таким образом, можно вычислить БПФ длиной $N$,
используя два ДПФ длиной $N/2$. Полный алгоритм БПФ заключается в рекурсивном
выполнении вышеописанной процедуры, начиная с объединения одиночных элементов в пары,
затем в четверки и так до полного охвата исходного массива данных. Дискретное преобразование Фурье для двумерного массива чисел размера
$M\times N$ определяется следующим образом: $$
G_{uw} =\frac{1}{NM}\sum\limits_{n=1}^{N-1} {\sum\limits_{m=1}^{M-1} {x_{mn}
} } e^{{-2\pi j\left[ {\frac{mu}{M}+\frac{nw}{N}} \right]} },
$$ а обратное преобразование $$
x_{mn} =\sum\limits_{u=1}^{N-1} {\sum\limits_{w=1}^{M-1} {G_{uw} } }
e^{ {2\pi j\left[ {\frac{mu}{M}+\frac{nw}{N}}
\right]} }.
$$ В случае обработки изображений компоненты двумерного преобразования Фурье
называют $\textit{пространственными частотами}$. Важным свойством двумерного преобразования Фурье является возможность его
вычисления с использованием процедуры одномерного БПФ: $$
G_{uw} =\frac{1}{N}\sum\limits_{n=1}^{N-1} { \left[
{\frac{1}{M}\sum\limits_{m=0}^{M-1} {x_{mn} e^{\frac{-2\pi jmw}{M}}} } \right] } e^{\frac{-2\pi jnu}{N}},
$$ Здесь выражение в квадратных скобках есть одномерное преобразование строки
матрицы данных, которое может быть выполнено с одномерным БПФ. Таким
образом, для получения двумерного преобразования Фурье нужно сначала
вычислить одномерные преобразования строк, записать результаты в исходную
матрицу и вычислить
одномерные преобразования для столбцов полученной
матрицы. При вычислении двумерного преобразования Фурье низкие частоты будут
сосредоточены в углах матрицы, что не очень удобно для дальнейшей обработки
полученной информации. Для перевода получения представления двумерного
преобразования Фурье, в котором низкие частоты сосредоточены в центре
матрицы, можно выполнить простую процедуру, заключающуюся в умножении
исходных данных на $-1^{m+n}$. На рис. 16 показаны исходное изображение и его Фурье-образ. Полутоновое изображение и его Фурье-образ (изображения получены в системе LabVIEW) Свертка функций $s(t)$ и $r(t)$ определяется как $$
s\ast r\cong r\ast s\cong \int\limits_{-\infty }^{+\infty } {s(\tau)} r(t-\tau)d\tau .
$$ На практике приходится иметь дело с дискретной сверткой, в которой
непрерывные функции заменяются наборами значений в узлах равномерной сетки
(обычно берется целочисленная сетка): $$
(r\ast s)_j \cong \sum\limits_{k=-N}^P {s_{j-k} r_k }.
$$ Здесь $-N$ и $P$ определяют диапазон, за пределами которого $r(t) = 0$. При вычислении свертки с помощью преобразования Фурье используется свойство
преобразования Фурье, согласно которому произведение образов функций в
частотной области эквивалентно свертке этих функций во временн ой области. Для вычисления сверки необходимо преобразовать исходные данные в частотную
область, то есть вычислить их преобразование Фурье, перемножить результаты
преобразования и выполнить обратное преобразование Фурье, восстановив
исходное представление. Единственная тонкость в работе алгоритма связана с тем, что в случае
дискретного преобразования Фурье (в отличие от непрерывного) происходит
свертка двух периодических функций, то есть наши наборы значений задают
именно периоды этих функций, а не просто значения на каком-то отдельном
участке оси. То есть алгоритм считает, что за точкой $x_{N }$ идет не ноль, а
точка $x_{0}$, и так далее по кругу. Поэтому, чтобы свертка корректно
считалась, необходимо приписать к сигналу достаточно длинную
последовательность нулей. Линейные методы фильтрации относятся к числу хорошо структурированных
методов, для которых разработаны эффективные вычислительные схемы,
основанные на быстрых алгоритмах свертки и спектральном анализе. В общем
виде линейные алгоритмы фильтрации выполняют преобразование вида $$
f"(x,y) = \int\int f(\zeta -x, \eta -y)K (\zeta , \eta) d \zeta d \eta ,
$$ где $K(\zeta ,\eta)$ - ядро линейного преобразования. При дискретном представлении сигнала интеграл в данной формуле вырождается
во взвешенную сумму отсчетов исходного изображения в пределах некоторой
апертуры. При этом выбор ядра $K(\zeta ,\eta)$ в соответствии с тем или
иным критерием оптимальности может привести к ряду полезных свойств
(гауссовское сглаживание при регуляризации задачи численного
дифференцирования изображения и др.). Наиболее эффективно линейные методы обработки реализуются в частотной
области. Использование Фурье-образа изображения для выполнения операций фильтрации
обусловлено прежде всего более высокой производительностью таких операций.
Как правило, выполнение прямого и обратного двумерного преобразования Фурье
и умножение на коэффициенты Фурье-образа фильтра занимает меньше времени,
чем выполнение двумерной свертки исходного изображения. Алгоритмы фильтрации в частотной области основываются на теореме о свертке.
В двумерном случае преобразование свертки выглядит следующим образом: $$
G\left({u,v} \right)=H\left({u,v} \right)F\left({u,v} \right),
$$ где $G$ - Фурье-образ результата свертки, $H$ - Фурье-образ фильтра, а
$F$ - Фурье-образ исходного изображения. То есть в частотной области двумерная
свертка заменяется поэлементным перемножением образов исходного изображения
и соответствующего фильтра. Для выполнения свертки необходимо выполнить следующие действия. Связь между функцией фильтра в частотной и пространственной области можно
определить, используя теорему о свертке $$
\Phi \left[ {f\left({x,y} \right)\ast h(x,y)} \right]=F\left({u,v}
\right)H\left({u,v} \right),
$$ $$
\Phi \left[ {f\left({x,y} \right)h(x,y)} \right]=F\left({u,v} \right)\ast
H\left({u,v} \right).
$$ Свертка функции с импульсной функцией может быть представлена следующим
образом: $$
\sum\limits_{x=0}^M {\sum\limits_{y=0}^N {s\left({x,y} \right)}
} \delta \left({x-x_0 ,y-y_0 } \right)=s(x_0 ,y_0).
$$ Фурье-преобразование импульсной функции $$
F\left({u,v} \right)=\frac{1}{MN}\sum\limits_{x=0}^M
{\sum\limits_{y=0}^N {\delta \left({x,y} \right) } } e^{ {-2\pi
j\left({\frac{ux}{M}+\frac{vy}{N}}
\right)} }
=\frac{1}{MN}.
$$ Пусть $f(x,y) = \delta (x,y)$, тогда свертка $$
f\left({x,y} \right)\ast h(x,y)=\frac{1}{MN}h\left({x,y} \right),
$$ $$
\Phi \left[ {\delta \left({x,y} \right)\ast h(x,y)} \right]=\Phi
\left[ {\delta \left({x,y} \right)} \right]H\left({u,v}
\right)=\frac{1}{MN}H\left({u,v} \right).
$$ Из этих выражений видно, что функции фильтра в частотной и пространственной
областях взаимосвязаны через преобразование Фурье. Для данной функции
фильтра в частотной области всегда можно найти соответствующий фильтр в
пространственной области, применив обратное преобразование Фурье. То же верно
и для обратного случая. Используя данную взаимосвязь, можно определить
процедуру синтеза пространственных линейных фильтров.
{$\textit{Идеальный фильтр низких частот}$} $H(u,v)$ имеет вид
$$H(u,v) = 1, \quad \mbox{если }D(u,v) < D_0 ,$$
$$H(u,v) = 0, \quad \mbox{если }D(u,v) \ge D_0 ,$$
где $D\left({u,v} \right)=\sqrt {\left({u-\frac{M}{2}} \right)^2+\left({v-\frac{N}{2}} \right)^2}$ - расстояние от центра частотной плоскости. {$\textit{Идеальный высокочастотный фильтр}$} получается путем
инверсии идеального низкочастотного фильтра: $$
H"(u,v) = 1-H(u,v).
$$ Здесь происходит полное подавление низкочастотных компонент при сохранении
высокочастотных. Однако как и в случае идеального низкочастотного фильтра,
его применение чревато появлением существенных искажений. Для синтеза фильтров с минимальными искажениями используются различные
подходы. Одним из них является синтез фильтров на основе экспоненты. Такие
фильтры привносят минимальные искажения в результирующее изображение и
удобны для синтеза в частотной области. Широко используемым при обработке изображений является семейство фильтров на
основании вещественной функции Гаусса. $\textit{Низкочастотный гауссовский фильтр}$ имеет вид $$
h\left(x \right)=\sqrt {2\pi } \sigma Ae^{-2\left({\pi \sigma
x} \right)^2} \mbox{ и } H\left(u \right)=Ae^{-\frac{u^2}{2\sigma ^2}}
$$ Чем уже профиль фильтра в частотной области (чем больше $\sigma $), тем он
шире в пространственной. {$\textit{Высокочастотный гауссовский фильтр}$} имеет вид $$
h\left(x \right)=\sqrt {2\pi } \sigma _A Ae^{-2\left({\pi \sigma _A x}
\right)^2}-\sqrt {2\pi } \sigma _B Be^{-2\left({\pi \sigma _B x} \right)^2 },
$$ $$
H\left(u \right)=Ae^{-\frac{u^2}{2\sigma _A^2 }}-Be^{-\frac{u^2}{2\sigma
_B^2 }}.
$$ В двумерном случае {$\it{низкочастотный}$} фильтр гаусса
выглядит следующим образом: $$
H\left({u,v} \right)=e^{-\frac{D^2\left({u,v} \right)}{2D_0^2 }}.
$$ {$\it{Высокочастотный}$} гауссовский фильтр имеет вид $$
H\left({u,v} \right)=1-e^{-\frac{D^2\left({u,v} \right)}{2D_0^2 }}.
$$ Рассмотрим пример фильтрации изображения (рис. 1) в частотной области
(рис. 17 - 22). Заметим, что частотная фильтрация изображения может
иметь смысл как сглаживания ($\textit{низкочастотная фильтрация}$), так и выделения контуров и мелкоразмерных
объектов ($\textit{высокочастотная фильтрация}$). Как видно из рис. 17, 19, по мере нарастания "мощности" фильтрации в низкочастотной
составляющей изображения все сильнее проявляется эффект "кажущейся расфокусировки" или
$\it{размытия}$ изображения. В то же время в высокочастотную составляющую, где в начале
наблюдаются лишь контура объектов, постепенно переходит большая часть
информационного содержания изображения (рис. 18, 20 - 22). Рассмотрим теперь поведение высокочастотных и низкочастотных фильтров
(рис. 23 - 28)
в присутствии аддитивного гауссовского шума на изображении (рис. 7). Как видно из рис. 23, 25, свойства низкочастотных фильтров по
подавлению аддитивной случайной помехи аналогичны свойствам ранее
рассмотренных линейных фильтров - при достаточной мощности фильтра помехи
подавляются, однако платой за это является сильное размытие контуров и
"расфокусировка" всего изображения. Высокочастотная составляющая
зашумленного изображения перестает быть информативной, так как помимо
контурной и объектовой информации там теперь также полностью присутствует и
шумовая компонента (рис. 27, 28). Применение частотных методов наиболее целесообразно в случае, когда известны
статистическая модель шумового процесса или/и оптическая передаточная
функция канала передачи изображения. Учесть такие априорные данные удобно,
выбрав в качестве восстанавливающего фильтра обобщенный управляемый
(параметрами $\sigma$ и $\mu$) фильтр следующего вида: $$
F(w_1,w_2)=
\left[
{ \frac
{1}
{P(w_1,w_2)}
}\right]
\cdot
\left[
{\frac
{{\vert P(w_1,w_2) \vert }^2}
{\vert P(w_1,w_2) \vert ^2 + \alpha \vert Q(w_1,w_2) \vert ^2}
}\right].
$$ где $0 < \sigma < 1$, $0 < \mu < 1$ - назначаемые параметры фильтра,
$P(w_{1}$, $w_{2})$ - передаточная функция системы, $Q(w_{1}$, $w_{2})$ -
стабилизатор фильтра, согласованный с энергетическим спектром фона. Выбор
параметров $\sigma = 1$, $\mu = 0$ приводит к чисто инверсной фильтрации,
$\sigma =\mu = 1$ к \it{винеровской фильтрации}, что позволяет получить изображение, близкое к
истинному в смысле минимума СКО при условии, что спектры
плотности мощности
изображения и его шумовой компоненты априорно известны. Для дальнейшего
улучшения эффекта сглаживания в алгоритм линейной (винеровской) фильтрации
вводят адаптацию, основанную на оценке локальных статистик: математического
ожидания $M(P)$ и дисперсии $\sigma (P)$. Этот алгоритм эффективно фильтрует
засоренные однородные поверхности (области) фона. Однако при попадании в
скользящее окно обработки неоднородных участков фона импульсная
характеристика фильтра сужается ввиду резкого изменения локальных статистик,
и эти неоднородности (контуры, пятна) передаются практически без
расфокусировки, свойственной неадаптивным методам линейной фильтрации. К достоинствам методов линейной фильтрации следует отнести их ясный
физический смысл и простоту анализа результатов. Однако при резком ухудшении
соотношения сигнал/шум, при возможных вариантах площадного зашумления и
наличии высокоамплитудного импульсного шума линейные методы предварительной
обработки могут оказаться недостаточными. В этой ситуации значительно более
мощными оказываются нелинейные методы. Я полагаю что все в общих чертах знают о существовании такого замечательного математического инструмента как преобразование Фурье. Однако в ВУЗах его почему-то преподают настолько плохо, что понимают как это преобразование работает и как им правильно следует пользоваться сравнительно немного людей. Между тем математика данного преобразования на удивление красива, проста и изящна. Я предлагаю всем желающим узнать немного больше о преобразовании Фурье и близкой ему теме того как аналоговые сигналы удается эффективно превращать для вычислительной обработки в цифровые. Без использования сложных формул и матлаба я постараюсь ответить на следующие вопросы: Я буду исходить из предположения что читатель понимает что такое интеграл , комплексное число (а так же его модуль и аргумент), свертка функций , плюс хотя бы “на пальцах” представляет себе что такое дельта-функция Дирака . Не знаете - не беда, прочитайте вышеприведенные ссылки. Под “произведением функций” в данном тексте я везде буду понимать “поточечное умножение” Начать надо, наверное, с того что обычное преобразование Фурье - это некая такая штука которая, как можно догадаться из названия, преобразует одни функции в другие, то есть ставит в соответствие каждой функции действительного переменного 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 полностью определяется её значениями на отрезке
=
∑
x[
n]
e−
jnk
(2
π
N
)
n=
0
N −
1
∑
X
[
k
]
e
jnk(2
π
N)
N k
=
0
Если по этим формулам разложить в спектр действительный сигнал, то первые N
/2+1 комплексных коэффициентов спектра будут совпадать со спектром «обычного» действительного ДПФ, представленным в «комплексном» виде, а остальные коэффициенты будут их симметричным отражением относительно Современную технику связи невозможно представить без спектрального анализа. Представление сигналов в частотной области необходимо как для анализа
их характеристик, так и для анализа блоков и узлов приемопередатчиков систем радиосвязи. Для преобразования сигналов в частотную область применяется
прямое преобразование Фурье. Обобщенная формула прямого преобразования Фурье записывается следующим образом: Как видно из этой формулы для частотного анализа производится вычисление корреляционной зависимости между сигналом, представленным во временной
области и комплексной экспонентой с заданной частотой. При этом по формуле Эйлера комплексная экспонента разлагается на реальную и мнимую часть: Сигнал, представленный в частотной области можно снова перевести во временное представление при помощи обратного преобразования Фурье. Обобщенная
формула обратного преобразования Фурье записывается следующим образом: В формуле прямого преобразования Фурье используется интегрирование по времени от минус бесконечности до бесконечности. Естественно это является
математической абстракцией. В реальных условиях мы можем провести интегрирование от данного момента времени, который мы можем обозначить за 0, до
момента времени T. Формула прямого преобразования Фурье при этом будет преобразована к следующему виду: В результате существенно меняются свойства преобразования Фурье
. Спектр сигнала вместо непрерывной функции становится дискретным рядом
значений
. Теперь минимальной частотой и одновременно шагом частотных значений спектра сигнала становится: Только функции sin и cos c частотами k/T
будут взаимно ортогональны, а это является непременным условием преобразования Фурье. Набор первых
функций разложения в ряд Фурье приведен на рисунке 1. При этом длительность функций совпадает с длительностью анализа T
. Теперь спектр сигнала будет выглядеть так, как это показано на рисунке 2. В данном случае формула вычисления прямого преобразования Фурье (4) преобразуется к следующему виду: Формула обратного преобразования Фурье для случая определения спектра на ограниченном отрезке времени будет выглядеть следующим образом: Подобным образом можно определить формулу прямого преобразования Фурье для цифровых отсчетов сигнала. Учитывая, что вместо непрерывного сигнала
используются его цифровые отсчеты, в выражении (6) интеграл заменяется на сумму. В данном случае длительность анализируемого сигнала определяется
количеством цифровых отсчетов N
. Преобразование Фурье для цифровых отсчетов сигнала называется дискретным преобразованием Фурье
и
записывается следующим образом: Теперь рассмотрим как изменились свойства дискретного преобразования Фурье (ДПФ) по сравнению с прямым преобразованием Фурье на ограниченном
интервале времени. Когда мы рассматривали дискретизацию аналогового сигнала, мы выяснили, что спектр входного сигнала должен быть ограничен по
частоте. Это требование ограничивает количество дискретных составляющих спектра сигнала. Первоначально может показаться, что мы можем ограничить
спектр сигнала частотой f
д /2, что соответствует количеству частотных составляющих K = N
/2 . Однако это не так.
Несмотря на то, что спектр сигнала для действительных отсчетов сигнала для положительных частот и отрицательных частот симметричен относительно 0,
отрицательные частоты могут потребоваться для некоторых алгоритмов работы со спектрами, например, для . Еще больше отличие получается при выполнении дискретного преобразования Фурье над комплексными отсчетами
входного сигнала. В результате для полного описания спектра цифрового сигнала требуется N
частотных отсчетов (k = 0, ..., N/2
). Многие сигналы удобно анализировать, раскладывая их на синусоиды (гармоники). Тому есть несколько причин. Например, подобным образом работает человеческое ухо. Оно раскладывает звук на отдельные колебания различных частот. Кроме того, можно показать, что синусоиды являются "собственными функциями" линейных систем (т.к. они проходят через линейные системы, не изменяя формы, а могут изменять лишь фазу и амплитуду). Еще одна причина в том, что теорема Котельникова формулируется в терминах спектра сигнала. Преобразование Фурье (Fourier transform) - это разложение функций на синусоиды (далее косинусные функции мы тоже называем синусоидами, т.к. они отличаются от "настоящих" синусоид только фазой). Существует несколько видов преобразования Фурье. 1. Непериодический непрерывный сигнал можно разложить в интеграл Фурье. 2. Периодический непрерывный сигнал можно разложить в бесконечный ряд Фурье. 3. Непериодический дискретный сигнал можно разложить в интеграл Фурье. 4. Периодический дискретный сигнал можно разложить в конечный ряд Фурье. Компьютер способен работать только с ограниченным объемом данных, следовательно, реально он способен вычислять только последний вид преобразования Фурье. Рассмотрим его подробнее. До сих пор мы рассматривали ДПФ от действительных сигналов. Обобщим теперь ДПФ на случай комплексных сигналов. Пусть x[n], 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 - его спектр): Непосредственное вычисление двумерного ДПФ по приведенным формулам требует огромных вычислительных затрат. Однако можно доказать, что двумерное ДПФ обладает свойством сепарабельности, т.е. его можно вычислить последовательно по двум измерениям. Для вычисления двумерного ДПФ достаточно вычислить одномерные комплексные ДПФ всех строк изображения, а затем вычислить в результирующем "изображении" одномерные комплексные ДПФ всех столбцов. При этом результаты всех одномерных комплексных ДПФ нужно записывать на место исходных данных для этих ДПФ. Например, при вычислении одномерного ДПФ первой строки изображения нужно результат ДПФ записать в первую строку этого изображения (он имеет тот же размер). Для этого нужно каждый "пиксель" хранить в виде комплексного числа. Таким образом, эффективный алгоритм вычисления ДПФ изображения заключается в вычислении одномерных БПФ сначала от всех строк, а потом - от всех столбцов изображения.Комплексное представление преобразования Фурье.
Быстрое преобразование Фурье.
Двумерное преобразование Фурье.
Свертка с использованием преобразования Фурье.
Фильтрация изображений в частотной области.
(3)
(4)
, (5)
Рисунок 1. Функции разложения в ряд Фурье
Рисунок 2. Спектр функции x
(t
) при анализе на ограниченном интервале времени
(6)
(7)
Комплексное ДПФ
Двумерное ДПФ