Дискретное преобразование фурье литература. Дискретное преобразование Фурье на VB.NET
Линейная фильтрация изображений может осуществляться как в пространственной, так и в частотной области. При этом считается, что "низким" пространственным частотам соответствует основное содержание изображения - фон и крупноразмерные объекты, а "высоким" пространственным частотам - мелкоразмерные объекты, мелкие детали крупных форм и шумовая компонента.
Традиционно для перехода в область пространственных частот используются методы, основанные на $\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)$. Этот алгоритм эффективно фильтрует
засоренные однородные поверхности (области) фона. Однако при попадании в
скользящее окно обработки неоднородных участков фона импульсная
характеристика фильтра сужается ввиду резкого изменения локальных статистик,
и эти неоднородности (контуры, пятна) передаются практически без
расфокусировки, свойственной неадаптивным методам линейной фильтрации. К достоинствам методов линейной фильтрации следует отнести их ясный
физический смысл и простоту анализа результатов. Однако при резком ухудшении
соотношения сигнал/шум, при возможных вариантах площадного зашумления и
наличии высокоамплитудного импульсного шума линейные методы предварительной
обработки могут оказаться недостаточными. В этой ситуации значительно более
мощными оказываются нелинейные методы. Даётся программный код для прямого и обратного преобразования Фурье. Рассматривается быстрое преобразование Фурье. Дискретное преобразование Фурье (ДПФ) - это мощный инструмент анализа, который широко используется в области цифровой обработки сигналов (ЦОС). Существуют прямое и обратное преобразования Фурье. Прямое дискретное преобразование Фурье переводит сигнал из временной области в частотную и служит для анализа частотного спектра сигнала. Обратное преобразование делает ровно противоположное: по частотному спектру сигнала восстанавливает сигнал во временной области. Для расчёта преобразования Фурье обычно используется ускоренная процедура расчёта - т.н. быстрое преобразование Фурье (БПФ). Это позволяет в значительной мере сократить процессорное время на достаточно сложные и ресурсоёмкие математические расчёты. Для начала нам потребуется вспомогательный класс, который будет описывать комплексные числа. Комплексные числа - это особый вид чисел в математике. Каждое комплексное число состоит из двух частей - действительной и мнимой. Сейчас нам достаточно знать о комплексных числах применительно к ДПФ то, что действительная часть комплексного числа хранит информацию об амплитуде сигнала, а мнимая - о фазе. На вход функции передаётся массив комплексных чисел. Действительная часть которого представляет произвольный дискретный сигнал, с отсчётами через равные промежутки времени. Мнимая часть содержит нули. Число отсчётов в сигнале должно равняться степени двойки. Если ваш сигнал короче, то дополните его нулями до числа, кратного степени 2: 256, 512, 1024 и т.д. Чем длиннее сигнал, тем у рассчитанного спектра будет выше разрешение по частоте. Обратное дискретное преобразование Фурье (ОДПФ) одним из этапов расчёта включает в себя прямое ДПФ на массиве комплексных чисел, где мнимая часть - это инверсия относительно оси X мнимой части спектра. Ну и конечно же, опишем использовавшийся метод, который проверяет число элементов переданного массива:
""" Теперь давайте проверим, что наши функции работают. Для этого пропустим произвольный сигнал через механизм прямого преобразования Фурье, а затем «соберём» его обратно с помощью обратного преобразования Фурье. Восстановленный сигнал должен практически совпадать с исходным. Ошибки округления, возникающие при работе с числами в компьютере, имеют место быть, поэтому сигналы не будут идентичны полностью, но их отклонение друг от друга должно быть пренебрежимо малым. Для примера в качестве исходного сигнала возьмём функцию синуса и сформируем данные длиной 128 отсчётов вот таким образом: Dim cn(127) As ComplexNumber
For i As Integer = 0 To cn.Length - 1
cn(i) = New ComplexNumber(Math.Sin(i * 3 * Math.PI / 180))
Next
Получим вот такой сигнал: Здесь по оси X - номера отсчётов во временной области, по оси Y - амплитуда. Обратим внимание, что сигнал состоит только из действительных частей, а мнимая часть на всём отрезке равна "0". Теперь передадим этот сигнал на вход функции FFT(). По полученным в ходе прямого преобразования Фурье массивам комплексных чисел построим два графика - действительной (Re) и мнимой (Im) частей спектра: Здесь по оси X - отсчёты в частотной области, по оси Y - амплитуда. Чтобы получить реальные значения частоты, необходимо рассчитать их, учитывая, что "0" оси Y соответствует нулевой частоте, максимум оси Y соответствует частоте дискретизации. Полученный спектр сигнала передадим функции обратного преобразования Фурье IFFT(). Получим массив комплексных чисел, где действительная часть будет содержать восстановленный сигнал: Как видно, восстановленный сигнал полностью повторяет исходный. Это одно из преобразований Фурье, широко применяемых в алгоритмах цифровой обработки сигналов (его модификации применяются в сжатии звука в MP3, сжатии изображений в JPEG и др.), а также в других областях, связанных с анализом частот в дискретном (к примеру, оцифрованном аналоговом) сигнале. Дискретное преобразование Фурье требует в качестве входа дискретную функцию. Такие функции часто создаются путём дискретизации (выборки значений из непрерывных функций). Дискретные преобразования Фурье помогают решать частные дифференциальные уравнения и выполнять такие операции, как свёртки. Дискретные преобразования Фурье также активно используются в статистике, при анализе временных рядов. Преобразования бывают одномерные, двумерные и даже трёхмерные. Прямое преобразование: Обратное преобразование: Обозначения: § N
- количество значений сигнала, измеренных за период, а также количество компонент разложения; § - измеренные значения сигнала (в дискретных временных точках с номерами , которые являются входными данными для прямого преобразования и выходными для обратного; § - N
комплексных амплитуд синусоидальных сигналов, слагающих исходный сигнал; являются выходными данными для прямого преобразования и входными для обратного; поскольку амплитуды комплексные, то по ним можно вычислить одновременно и амплитуду, и фазу; § - обычная (вещественная) амплитуда k-го синусоидального сигнала; § arg(X k
) - фаза k-го синусоидального сигнала (аргумент комплексного числа); § k
- частота k-го сигнала, равная , где T
- период времени, в течение которого брались входные данные. Из последнего видно, что преобразование раскладывает сигнал на синусоидальные составляющие (которые называются гармониками) с частотами от N колебаний за период до одного колебания за период. Поскольку частота дискретизации сама по себе равна N отсчётов за период, то высокочастотные составляющие не могут быть корректно отображены - возникает муаров эффект. Это приводит к тому, что вторая половина из N комплексных амплитуд, фактически, является зеркальным отображением первой и не несёт дополнительной информации. Рассмотрим некоторый периодический сигнал x
(t
) c периодом равным T. Разложим его в ряд Фурье:
Проведем дискретизацию сигнала так, чтобы на периоде было N отсчетов. Дискретный сигнал представим в виде отсчетов: x n
= x
(t n
), где , тогда эти отсчеты через ряд Фурье запишутся следующим образом:
Используя соотношение: , получаем:
Таким образом, мы получили обратное дискретное преобразование Фурье.
Умножим теперь скалярно выражение для x n
на и получим: Здесь использованы: а) выражение для суммы конечного числа членов (экспонент) геометрической прогрессии, и б) выражение символа Кронекера как предела отношения функций Эйлера для комплексных чисел. Отсюда следует, что: Эта формула описывает прямое дискретное преобразование Фурье
. В литературе принято писать множитель в обратном преобразовании, и поэтому обычно пишут формулы преобразования в следующем виде: Дискретное преобразование Фурье является линейным преобразованием, которое переводит вектор временных отсчётов в вектор спектральных отсчётов той же длины. Таким образом, преобразование может быть реализовано как умножение квадратной матрицы на вектор:
Преобразование Фурье
При использовании преобразований Фурье изображение представляется в виде суммы сложных показательных функций переменных амплитуды, частоты и фазы. Преобразование Фурье играет очень важную роль во многих областях обработки изображений, включая улучшение, анализ, восстановление и сжатие. Основные определения преобразования Фурье
Если ƒ(m,n)
представляет собой функцию двух дискретных пространственных переменных m и n, тогда двумерное преобразование Фурье функции ƒ(m,n)
может быть представлено следующим выражением Переменные представляют собой угловые частоты. Таким образом, представляет собой функцию ƒ(m,n)
в частотной области. является комплекснозначной функцией с соответствующими частотами . Частоты находятся в пределах диапазона , . Отметим, что F
(0,0) представляется в виде суммы всех переменных ƒ(m,n)
. По этой причине F
(0,0) часто называют постоянной составляющей преобразования Фурье. Обратное двумерное преобразование Фурье представляется выражением Т.е. это выражение представляет ƒ(m,n)
в виде суммы бесконечного числа сложных экспоненциальных функций (синусоид) с различными частотами. Амплитуда и фаза определяют вклад частот в представление . Визуализация Фурье-преобразования
При иллюстрации Фурье-преобразования допустим, что функция ƒ(m,n)
равна 1 и представлена в виде прямоугольника. Для упрощения диаграммы, функция ƒ(m,n)
будет представляться непрерывной функцией двух дискретных переменных m
и n
. На рисунке внизу с использованием функции mesh визуализировано значения амплитуд, которые получены при Фурье-преобразовании прямоугольной функции, представленной на предыдущем рисунке. Визуализацию амплитуды еще называют визуализацией преобразований Фурье. Пик функции находится в центре и отображает значение F
(0,0), которое является суммой всех значений ƒ(m,n)
. Все остальные составляющие представляют собой распределение энергии по вертикальным и горизонтальным частотам. Другой путь визуализации Фурье-преобразования заключается в отображении значений в виде изображения. Рассмотрим примеры преобразования Фурье функций различных простых форм. Дискретное косинусное преобразование
Дискретные косинусные преобразования представляют изображение в виде суммы синусоид с различной амплитудой и частотой. Функция dct2 в приложении Image Processing Toolbox реализует двумерные дискретные косинусные преобразования изображений. Одна из особенностей дискретного преобразования Фурье состоит в том, что некоторые локальные участки изображения можно охарактеризовать небольшим количеством коэффициентов дискретного преобразования Фурье. Это свойство очень часто используется при разработке методов сжатия изображений. Например, дискретное косинусное преобразование является основой международного стандарта, который используется в алгоритме сжатия изображений с потерями JPEG. Название формата “JPEG” состоит из первых букв названия рабочей группы, которая принимала участие в разработке этого стандарта (Joint Photographic Experts Group). Двумерное дискретное косинусное преобразование матрицы A
с размерами реализуется согласно следующему выражению Значения B pq
называют коэффициентами дискретного косинусного преобразования матрицы A
. (Следует отметить, что индексы матрицы в MATLAB всегда начинаются с 1, а не с 0. Поэтому элементы матрицы, которые представлены в MATLAB как A(1,1) и B(1,1), будут соответствовать элементам A 00
и B 00
из приведенной выше формулы.) Обратное дискретное косинусное преобразование реализуется согласно выражениям Выражение обратного дискретного косинусного преобразования может интерпретироваться как представление матрицы A
с размерами в виде суммы следующих функций Эти функции называются основными (базовыми) функциями дискретного косинусного преобразования. Коэффициенты дискретного косинусного преобразования B pq
можно рассматривать как весовые при каждой базовой функции. Например, для матрицы с размером элементов существует 64 базовые функции, что продемонстрировано на изображении. Горизонтальные частоты увеличиваются слева направо, а вертикальные – сверху вниз. Матрица дискретных косинусных преобразований
Приложение Image Processing Toolbox предлагает два различных пути реализации дискретных косинусных преобразований. Первый метод реализован в функции dct2. Функция dct2 использует быстрое преобразования Фурье для ускорения вычислений. Второй метод использует матрицу дискретных косинусных преобразований, которая возвращается функцией dctmtx. Матрица преобразований T формируется согласно следующего выражения Для матрицы A
с размерами представляет собой матрицу с размерами , где каждый столбец содержит одномерное дискретное косинусное преобразование A
. Двумерное дискретное косинусное преобразование A
вычисляется как B=T*A*T’
. Обратное двумерное дискретное косинусное преобразование B
вычисляется как T’*B*T.
Дискретные косинусные преобразования и сжатие изображений
В алгоритме сжатия изображений JPEG исходное изображение разделяется на блоки с размерами или элементов. Далее для каждого блока вычисляется двумерное дискретное косинусное преобразование. Коэффициенты дискретных косинусных преобразований квантируются, кодируются и передаются. Получатель JPEG–данных декодирует коэффициенты дискретного косинусного преобразования, вычисляет обратное двумерное дискретное косинусное преобразование в каждом блоке и далее совмещает их вместе в одно изображение. Рассмотрим пример вычисления двумерных дискретных косинусных преобразований в блоках с размерами элементов исходного изображения. Далее при реконструкции изображения будем учитывать только 10 коэффициентов из каждого блока, остальные приравняем к нулю. При проведении описанных вычислений будет использоваться также матрица преобразований. I = imread("cameraman.tif");
I = im2double(I);
T = dctmtx(8);
B = blkproc(I,,"P1*x*P2",T,T");
mask = ;
B2 = blkproc(B,,"P1.*x",mask);
I2 = blkproc(B2,,"P1*x*P2",T",T);
imshow(I);
figure, imshow(I2)
На рисунке представлено два изображения – исходное и реконструированное. При реконструкции изображения использовалось только 15 % коэффициентов дискретных косинусных преобразований. Однако, следует отметить, что качество реконструированного изображения является довольно приемлемым. Для просмотра других свойств дискретного косинусного преобразования см. функцию dctdemo. Преобразования Радона
Функция radon в приложении Image Processing Toolbox вычисляет матрицу проекций изображения вдоль заданных направлений. Проекция двумерной функции f(x,y) равна интегралу вдоль указанной линии. Функция Радона представляет собой вычисление проекций изображения на оси, которые задаются углами в градусах относительно горизонтали против часовой стрелки. На рисунке показана проекция некоторой фигуры под указанным углом На рисунке внизу показаны горизонтальные и вертикальные проекции для простой двумерной функции. Проекции могут вычисляться вдоль произвольного угла theta. Встроенная в приложение Image Processing Toolbox функция radon вычисляет проекции изображения вдоль определенных направлений. Проекция двумерной функции f(x,y) на ось x’ представляет собой линейный интеграл Таким образом, оси x’ y’ задаются поворотом на угол против часовой стрелки. На изображении внизу проиллюстрировано геометрию преобразования Радона. Визуализация преобразований Радона
При проведении преобразований Радона необходимо указать исходное изображение и вектор углов theta. Radon(I,theta);
R
представляет собой матрицу, в которой каждый столбец является преобразованием Радона для одного из углов, который содержится в векторе theta. Вектор xp содержит соответствующие координаты вдоль оси x. Центральный пиксель I определяется согласно выражению floor((size(I)+1)/2). Рассмотрим, как в преобразованиях Радона вычисляются проекции. Рассмотрим проекции под углом 0° и 45°. I = zeros(100,100);
I(25:75, 25:75) = 1;
imshow(I)
Radon(I,);
figure; plot(xp,R(:,1)); title("R_{0^o} (x\prime)")
Преобразования Радона при 0° Figure; plot(xp,R(:,2)); title("R_{45^o} (x\prime)")
Преобразования Радона при большом числе углов часто отображается в виде изображения. В данном примере рассмотрено преобразования Радона для изображения в виде квадрата при диапазоне углов от 0° до 180° с дискретностью 1°. Theta = 0:180;
= radon(I,theta);
imagesc(theta,xp,R);
title("R_{\theta} (X\prime)");
xlabel("\theta (degrees)");
ylabel("X\prime");
set(gca,"XTick",0:20:180);
colormap(hot);
colorbar
Использование преобразований Радона при детектировании линий
Преобразования Радона аналогичны другим известным операциям, которые известны как преобразования Хоха. Функцию radon можно применять для детектирования прямых линий. Рассмотрим основные этапы этого процесса. Наибольший пик в матрице R
соответствует =1° и x´= -80. Из центра исходного изображения проводится линия под углом на расстояние x’. Перпендикулярно к этой линии проводится прямая, которая соответствует прямой на исходном изображении. Кроме того, на изображении присутствуют и другие линии, которые представлены в матрице R
соответствующими пиками. Современную технику связи невозможно представить без спектрального анализа. Представление сигналов в частотной области необходимо как для анализа
их характеристик, так и для анализа блоков и узлов приемопередатчиков систем радиосвязи. Для преобразования сигналов в частотную область применяется
прямое преобразование Фурье. Обобщенная формула прямого преобразования Фурье записывается следующим образом: Как видно из этой формулы для частотного анализа производится вычисление корреляционной зависимости между сигналом, представленным во временной
области и комплексной экспонентой с заданной частотой. При этом по формуле Эйлера комплексная экспонента разлагается на реальную и мнимую часть: Сигнал, представленный в частотной области можно снова перевести во временное представление при помощи обратного преобразования Фурье. Обобщенная
формула обратного преобразования Фурье записывается следующим образом: В формуле прямого преобразования Фурье используется интегрирование по времени от минус бесконечности до бесконечности. Естественно это является
математической абстракцией. В реальных условиях мы можем провести интегрирование от данного момента времени, который мы можем обозначить за 0, до
момента времени T. Формула прямого преобразования Фурье при этом будет преобразована к следующему виду: В результате существенно меняются свойства преобразования Фурье
. Спектр сигнала вместо непрерывной функции становится дискретным рядом
значений
. Теперь минимальной частотой и одновременно шагом частотных значений спектра сигнала становится: Только функции sin и cos c частотами k/T
будут взаимно ортогональны, а это является непременным условием преобразования Фурье. Набор первых
функций разложения в ряд Фурье приведен на рисунке 1. При этом длительность функций совпадает с длительностью анализа T
. Теперь спектр сигнала будет выглядеть так, как это показано на рисунке 2. В данном случае формула вычисления прямого преобразования Фурье (4) преобразуется к следующему виду: Формула обратного преобразования Фурье для случая определения спектра на ограниченном отрезке времени будет выглядеть следующим образом: Подобным образом можно определить формулу прямого преобразования Фурье для цифровых отсчетов сигнала. Учитывая, что вместо непрерывного сигнала
используются его цифровые отсчеты, в выражении (6) интеграл заменяется на сумму. В данном случае длительность анализируемого сигнала определяется
количеством цифровых отсчетов N
. Преобразование Фурье для цифровых отсчетов сигнала называется дискретным преобразованием Фурье
и
записывается следующим образом: Теперь рассмотрим как изменились свойства дискретного преобразования Фурье (ДПФ) по сравнению с прямым преобразованием Фурье на ограниченном
интервале времени. Когда мы рассматривали дискретизацию аналогового сигнала, мы выяснили, что спектр входного сигнала должен быть ограничен по
частоте. Это требование ограничивает количество дискретных составляющих спектра сигнала. Первоначально может показаться, что мы можем ограничить
спектр сигнала частотой f
д /2, что соответствует количеству частотных составляющих K = N
/2 . Однако это не так.
Несмотря на то, что спектр сигнала для действительных отсчетов сигнала для положительных частот и отрицательных частот симметричен относительно 0,
отрицательные частоты могут потребоваться для некоторых алгоритмов работы со спектрами, например, для . Еще больше отличие получается при выполнении дискретного преобразования Фурье над комплексными отсчетами
входного сигнала. В результате для полного описания спектра цифрового сигнала требуется N
частотных отсчетов (k = 0, ..., N/2
).Комплексное представление преобразования Фурье.
Быстрое преобразование Фурье.
Двумерное преобразование Фурье.
Свертка с использованием преобразования Фурье.
Фильтрация изображений в частотной области.
1
Комплексные
числа
2
Прямое дискретное быстрое
преобразование Фурье
3
Обратное дискретное быстрое
преобразование Фурье
4
Проверка прямого и обратного
преобразования Фурье
где 
Прямоугольная функция
Амплитуда изображения прямоугольной функции
Логарифмическое представление Фурье-преобразования прямоугольной функции
Примеры Фурье преобразования функций различных простых форм
64 базовые функции, которые получены для матрицы с размерами элементов
Параллельно-лучевая проекция с углом поворота theta
Горизонтальная и вертикальная проекции некоторой простой функции
Геометрия преобразования Радона
Преобразования Радона при 45°
Преобразования радона с использованием 180 проекций
Геометрия преобразования Радона при детектировании прямых линий
(3)
(4)
, (5)
Рисунок 1. Функции разложения в ряд Фурье
Рисунок 2. Спектр функции x
(t
) при анализе на ограниченном интервале времени
(6)
(7)
Смерть королевы Елизаветы II станет катастрофой для Великобритании
Образование нашей планеты произошло
Как найти свое призвание Интеллектуальная игра как найти свое призвание