Болезни Военный билет Призыв

Обратное дискретное преобразование фурье. Преобразование Фурье. Линейная фильтрация в частотной области

Линейная фильтрация изображений может осуществляться как в пространственной, так и в частотной области. При этом считается, что "низким" пространственным частотам соответствует основное содержание изображения - фон и крупноразмерные объекты, а "высоким" пространственным частотам - мелкоразмерные объекты, мелкие детали крупных форм и шумовая компонента.

Традиционно для перехода в область пространственных частот используются методы, основанные на $\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$ - Фурье-образ исходного изображения. То есть в частотной области двумерная свертка заменяется поэлементным перемножением образов исходного изображения и соответствующего фильтра.

Для выполнения свертки необходимо выполнить следующие действия.

  1. Умножить элементы исходного изображения на $-1^{m+n}$, для центрирования Фурье-образа.
  2. Вычислить Фурье образ $F(u,v)$, используя БПФ.
  3. Умножить Фурье образ $F(u,v)$ на частотную функцию фильтра $H(u,v)$.
  4. Вычислить обратное преобразование Фурье.
  5. Умножить вещественную часть обратного преобразования на $-1^{m+n}$.

Связь между функцией фильтра в частотной и пространственной области можно определить, используя теорему о свертке

$$ \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). $$

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

  1. Определяем требуемые характеристики (форму) фильтра в частотной области.
  2. Выполняем обратное преобразование Фурье.
  3. Полученный фильтр можно использовать как маску для пространственной свертки, при этом размеры маски можно уменьшить по сравнению с размерами исходного фильтра.

{$\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)$. Этот алгоритм эффективно фильтрует засоренные однородные поверхности (области) фона. Однако при попадании в скользящее окно обработки неоднородных участков фона импульсная характеристика фильтра сужается ввиду резкого изменения локальных статистик, и эти неоднородности (контуры, пятна) передаются практически без расфокусировки, свойственной неадаптивным методам линейной фильтрации.

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

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

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

(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 ).

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

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

Преобразование Фурье (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 комплексных коэффициентов спектра будут совпадать со спектром «обычного» действительного ДПФ, представленным в «комплексном» виде, а остальные коэффициенты будут их симметричным отражением относительно

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

Пусть имеется дискретный сигнал х Д (t) , состоящий из N отсчетов х к , и дискретный сигнал у д (Г), состоящий из N отсчетов у к, тогда дискретной сверткой этих двух сигналов называется сигнал z A (t) , для которого

Дискретные сигналы получили широкое распространение при создании систем с импульсной модуляцией.

Устройство дискретизации в простейшем случае представляет собой стробируемый каскад (ключ), открывающийся на время т и с периодом А (рис. 4.7).


Рис. 4.

Интервал дискретизации А может быть постоянным (равномерная дискретизация) или переменным (адаптивная дискретизация). Наиболее распространенной формой дискретизации является равномерная, в основе которой лежит теорема Котельникова.

Импульсный модулятор - это устройство с двумя входами, на один из которых подается аналоговый сигнал, а на второй поступают короткие синхронизирующие импульсы с периодом повторения А. При этом в момент поступления синхроимпульса происходит измерение мгновенного значения сигнала лс(г). На выходе модулятора возникает последовательность импульсов, каждый из которых имеет площадь, пропорциональную соответствующему отсчетному значению аналогового сигнала (рис. 4.7).

Сигнал Хмпн (t ) на выходе импульсного модулятора называют модулированной импульсной последовательностью (МИП). Математически МИП записывается так

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

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

Рассмотрим особенности спектрального представления дискретного сигнала, заданного на интервале своими отсчетами x 0 ,x x ,...,x N _ x . Полное число отсчетов N - Т / А.

Методика изучения таких дискретных сигналов состоит в том, что полученная выборка отсчетных значений мысленно повторяется бесконечное число раз. В результате сигнал становится периодическим (рис. 4.8).

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


Рис. 4.8.

Запишем модель ограниченного периодического сигнала в виде последовательности дельта-импульсов:

Разложим сигнал Хмип (0 в ряд Фурье:

здесь замена переменных? = f / А. Окончательно получаем

Эта формула определяет последовательность коэффициентов, образующих дискретное преобразование Фурье (ДПФ) рассматриваемого сигнала.

ДПФ обладает следующими свойствами:

1. ДПФ есть линейное преобразование, т. е. если z k = а х к + /? у к, то

С"Z П ~ ^ С Х п Р Су п .

2. Число различных коэффициентов Cq,Ci,...,C n _i равно числу N отсчетов за период, при n = N коэффициент C N = С 0 .

3. С 0 является средним значением всех отсчетов С 0 = - к.

N к

  • 4. Если N- четное число, то С N = -^(-1) к х к.
  • 7 ^ ?=о
  • 5. Если отсчеты х к - вещественные числа и N - четное число, то C N = C* N , / = 0; Л/7 2 -1.
  • -+i - -i
  • 6. Если y k =x k+m , m = l;JV-l,TO C, t =C, * e ~ j2rrkm,N .
  • 2 tf-l
  • 7. Если z k = - > T0 C z к =C X k C y k

iy/ i =0

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

Если на основе отсчетов x 0 ,x l ,...,x N _ l некоторого сигнала найдены коэффициенты ДПФ C 0 ,Ci,... 9 C n/2 , то по ним можно восстановить аналоговый сигнал с ограниченным спектром x(t). Ряд Фурье такого сигнала имеет вид (при четном N)

где |Q| - модуль коэффициентов ДПФ; =arg - фазовый угол (аргумент)

коэффициентов ДПФ. Частота первой гармоники: f= - / в = - = -/i- нечетном N последнее слагаемое в формуле (4.17) равно:

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

Эта формула носит название обратного дискретного преобразования Фурье {ОДПФ).

Обозначим через

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

. (3.21)

Если сигнал существует только внутри прямоугольника со сторонами элементов (рис. 3.4.а), то сигнал определен на всей плоскости и является на ней прямоугольно-периодическим (рис. 3.4.б).

Рис. 3.4. Реальное (а) и периодически продолженное (б) изображения

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

Базисные функции этого двумерного представления - двумерные комплексные экспоненты (иногда называемые комплексными синусоидами)

(3.23)

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

Коэффициенты Фурье ряда (3.22) образуют двумерный частотный спектр сигнала и определяются формулой прямого преобразования Фурье:

(3.24)

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

Заметим, что для точного представления дискретного сигнала с двумерным периодом элементов согласно формулам БПФ достаточно конечного числа базисных функций (3.23) - ряд (3.22) является конечным. Это и понятно, поскольку сам представляемый сигнал содержит в одном периоде конечное число точек, т.е. имеет конечное число степеней свободы. Ясно, что число степеней свободы в спектре не может отличаться от числа степеней свободы в самом сигнале.

Остановимся на наиболее существенных свойствах двумерного дискретного спектра Фурье. Вычислим спектральные коэффициенты (3.24) в частотных точках :

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

,

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

Преобразуем соотношение (3.25) следующим образом. Во-первых, вместо частот подставим частоты . Во-вторых, выполним комплексное сопряжение обеих частей, что не нарушит равенства. В результате нетрудно получить выражение:

,

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

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

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