Fast Fourier Transform Algorithm. Часть 1. Простая теория

Fast Fourier Transform (FFT) или быстрое преобразование Фурье — один из самых важных алгоритмов в таких дисциплинах, как signal processing и data analysis.

Так получилось, что в нашем курсе математического анализа было уж чересчур формальное определение этого понятия, и поэтому я не придавал ему должного внимания. Мало ли, что есть в мат. анализе 🙂 Но все же краем уха слышал, что эта штука широко применяется при обработке сигналов. Тем не менее меня вряд ли бы «заинтересовало» преобразование Фурье, пока я лично не столкнулся с ним при обработке изображений.

В этой статье я постараюсь понятным языком объяснить смысл этого преобразования и привести его простейшую реализацию на языке Python.

1. Introduction: Why Fourier?

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

Основная идея, лежащая в этих «ужасных» формулах, очень проста, даже увлекательна 🙂 : любую функцию  f(x) можно представить как сумму ряда синусов и косинусов с увеличивающейся частотой. Другими словами, любые данные, изменяющиеся в пространстве или времени, могут быть представлены в другой области, называемой пространством частот (frequency space). Эта идея первому пришла в голову французскому математику и физику Joseph Fourier в 19 веке. С тех пор она оказалось полезной в различных областях, в основном в обработке сигнала.

 1.1 Frequency Space

Давайте поговорим немного о частотном пространстве перед тем, как углубимся в детали. Термин частота происходит из физики, как некоторое изменение во времени, описание характеристики некоторого периодического движения или поведения. Например, термин частота в компьютерном зрении обычно связан с вариацией яркости или цвета вдоль изображения, то есть является функцией пространственных координат, а не времени. Некоторые книги даже называют ее пространственной частотой (spatial frequency).

Например, если изображение, представленное в частотном пространстве, имеет высокие частоты, то это означает, что на нем присутствуют острые углы или детали. Давайте посмотрим на рисунки ниже, слева изображено исходное изображение, а справа — частотный график этой картинки.  

result_img1 result_img2 result_img3

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

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

2d-sine_result

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

1.2 So, What’s the Point?

Частотная область вносит существенные преимущества в обработку изображений (image processing). Она позволяет делать крупные операции фильтрации гораздо быстрее. Так как информация об изображении хранится в частотах, то это позволяет иногда отделить сигнал от шума, а также осуществлять различные вычисления, которые были бы очень трудоемкими в пространственной области. Кроме того Fourier Transform позволяет легко осуществлять переход из пространственной области в частотную и наоборот.

Например, допустим, у нас есть изображение с некоторым периодичным шумом, который мы хотели бы устранить. (Представьте себе отсканированную страницу с сероватыми полосами.) Если перевести данные изображения в частотное пространство, то любой периодический шум в исходной картинке будет отображаться как яркие пятна на диаграмме в частотном пространстве. Если мы устраним эти точки и применим inverse Fourier Transform, чтобы получить исходное изображение, мы сможем удалить большую часть шума и тем самым улучшить видимость изображения.

Начальное изображение:

noise_image

Посмотрим его диаграмму в частотном пространстве:

magnitude_spectrum

Так как у нас простой периодический шум, то он характеризуется, скорее всего, двумя пиками на диаграмме. Мы как раз видим два красных пика на оси абсцисс.

Посмотрим на нашу диаграмму в сером цвете:

magnitude_spectrum_gray

Удалим из диаграммы наш шум

magnitude_spectrum_cut

Сделаем обратное преобразование Фурье и посмотрим результат:

denoised_image

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

Дополнительные преимущества методов Фурье и их применение будут обсуждаться позднее в моих статьях.

2 Basics

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

2.1 Representing Complex Numbers

Комплексное число может быть записано как

R+iI

где  R  и  I — действительные числа, и  i  эквивалентно  \sqrt{-1}.  R  обозначает вещественную часть, а  I  обозначает мнимую часть комплексного числа. Вещественные числа можно рассматривать как подмножество комплексных чисел, где  I = 0.

Пример:  1+2i, где  R = 1, I = 2.

Говоря геометрически, вещественные числа могут поместиться на бесконечно длинной прямой в 1 мерном пространстве. Если мы добавим еще одно измерение к нему, то мы сможем представить намного больше чисел, то есть числа, которые лежат выше и ниже линии действительных чисел. Таким образом, комплексные числа лежат на плоскости, а не на линии. Горизонтальная ось такого представления называется вещественной осью, а вертикальная — мнимой осью. Поэтому мнимому числу  R+iI  соответсвуют координаты  (R,I)  на этой плоскости.

Пример: число  1+2i  находится в точке  (1,2)  комплексной плоскости.

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

r(cos \theta + isin \theta),

где  r = \sqrt{R^2+I^2}  и  \theta = \arctan{\frac{I}{R}} . Факт, что  cos \theta = \frac{R}{r}  и  sin \theta = \frac{I}{r}, делает очевидным, что две формы представления обозначают одно и то же комплексное число.

Пример:  1+2i  может быть представлено как  \sqrt{3}(cos \theta + isin \theta) , где  tan \theta = \frac{2}{1} = 2 , что определят  \theta \approx 64.435  градусов. Магнитуда равна  \sqrt{3}  и угол — примерно  64.435  градусов или  1.107  радиан.

2.2 Euler’s Formula

Формула Эйлера:

e^{i\theta} = cos\theta + isin\theta,

где  e = 2.718281828...,  и  \theta  — угол, который может быть любым вещественным числом. Это дает нам еще одно представление комплексных чисел в виде:

re^{i\theta},

где  r  — модуль полярной формы комплексного числа, и  \theta  — угол.

Пример:  1+2i  может также быть представлено как  \sqrt{3}e^{i\theta} , где  \theta \approx 64.435  градусов или  1.107  радиан.

В некоторых учебниках комплексное число часто выражается в виде формулы Эйлера без явного указания на это. (Это как раз случай в формулах, связанных с преобразованием Фурье.) Если вы увидите что-то в форме  re^{i\theta} , будьте уверены, что это обычное комплексное число. Кроме того,  \theta  обычно означает угол в радианах, если не указано иное.

3 Fourier Transform

Сначала мы рассмотрим преобразование Фурье с чисто математической точки зрения, то есть мы будем говорить о «непрерывных» и «бесконечных» вещах. Я предполагаю, что вы знаете, что означают математические символы  \pi  и  e , и знакомы с комплексными числами. Предстоящий раздел содержит комплексную математику, так что если вы страдаете «интеграл-о-фобией», то быстро пробегите до следующего раздела и посмотрите на дискретное преобразование Фурье. Помните, что преобразование Фурье функции есть суммирование синусов и косинусов разной частоты. Суммирование может, в теории, состоять из бесконечного числа членов синусов и косинусов.

3.1 Equations

Пусть  f(x) — непрерывная функция вещественной переменной  x. Преобразование Фурье (Fourier Transform) функции  f(x)  определяется уравнением:

F(u)=\int^{-\infty}_{\infty}f(x)e^{-i2\pi ux}\,dx

где  i=\sqrt{-1}  и  u  часто называется частотной переменной. Суммирование синусов и косинусов может быть не таким очевидным на первый взгляд, но, применяя уравнение Эйлера, получаем:

F(u)=\int^{-\infty}_{\infty}f(x)(cos 2\pi ux -i sin 2\pi ux)\,dx

По данному F(u) мы можем вернуться назад и получить f(x), используя обратное преобразование Фурье (inverse Fourier Transform):

f(x)=\int^{-\infty}_{\infty}F(u)e^{i2\pi ux}\,dx

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

Некоторые могут спросить, что же такое  F(u).  F(u) это как раз данные в частотном пространстве, о котором мы говорили в первой секции. Даже если мы начинаем с вещественной функции  f(x)  в пространственной области, мы обычно заканчиваем комплексными значениями  F(u). Это происходит потому, что вещественное число, умноженное на комплексное, дает комплексное число.

Итак,  f(x)e^{i2\pi ux}  — комплексное, поэтому сумма этих членов должна также дать комплексное число, то есть  F(u). Поэтому,

F(u)=R(u)+iI(u)

где R(u) является вещественной компонентой (члены, не содержащие  i) и  I(u) есть комплексная компонента (члены, содержащие i).  ``(u)" написано, чтобы напомнить, что члены являются функциями от u. Как любые другие комплексные числа, мы также можем записать их в полярной форме, дающей:

F(u)=r(sin \theta + i cos \theta)=re^{i\theta}

В большинстве учебников эта форма Fourier Transform записывается в виде:

F(u)=\left|F(u)\right|e^{i\theta(u)}

но это по сути одно и то же.

Существует несколько часто используемых слов, когда речь идет о преобразовании Фурье.  Модуль  \left|F(u)\right|  называется спектром Фурье (Fourier spectrum) функции  f(x)  и  \theta(u)  это фазовый угол (phase angle). Квадрат спектра,  \left|F(u)\right|^2=R^2(u)+I^2(u)  часто обозначается как  P(u)  и называется спектром мощности (power spectrum) функции  f(x). Термин плотность спектра (spectral density) также широко используется для обозначения спектра мощности. Спектр Фурье часто строится от  u. Спектр Фурье полезен, потому что он может быть легко построен как функция от  u  на бумаге :).

fourier_plot

На картинке изображена простая функция и ее спектр Фурье.

Заметьте, что сами  F(u)  трудно построить от  u  в 2-D пространстве, так как они комплексные числа.

3.2 Discrete Fourier Transform

Теперь, когда вы знаете, несколько вещей о преобразовании Фурье, мы должны выяснить способ, как использовать его на практике. Возвращаясь к примеру, где мы преобразовывали изображение, брав значения яркости пикселей, яркость не была непрерывной функцией. (Помните, что преобразование Фурье из прошлой секции было для непрерывной функции  f(x).) Математики придумали для этого хорошее решение, а именно дискретное преобразование Фурье (discrete Fourier transform):

Учитываются  N  дискретных значений  f(x), отобранных с единым шагом,

F(u)=\frac{1}{N}\sum_{x=0}^{N-1}f(x)e^{-i2\pi ux/N}

для  u=0,1,2,...,N-1,  и

f(x)=\sum_{u=0}^{N-1}F(u)e^{i2\pi ux/N}

для  x=0,1,2,...,N-1.

 Обратите внимание, что интеграл заменяется на суммирование, которое является простым «for loop» при программировании. Для тех, кому любопытно, расчет внутри  \sum  это произведение  f(x)=R+iI  на  e^{i2\pi ux/N}=cos(p)-isin(p)  где  p=2\pi ux/N:

f(x)=e^{i2\pi ux/N}=(R+iI)*(cos(p)-isin(p))

=Rcos(p)-Risin(p)+Iicos(p)-Ii^2sin(p)

=Rcos(p)-Risin(p)+Iicos(p)+Isin(p)

=(Rcos(p)+Isin(p))+i(Icos(p)-Rsin(p)),

где  R,  I  — вещественные числа, и  I=0, когда f(x)  — вещественное число.

4 Simple Implementation

В конце приведу простейшую реализацию discrete Fourier Transform на языке Python, но по сути ее можно легко переписать на любой другой язык:

import numpy as np
import scipy as sp

def dft_slow(fx): 
    fx = np.asarray(fx, dtype=float)
    N = fx.shape[0]
    result = np.zeros(N,dtype=complex)
    for u in xrange(N):
        z = np.zeros(2)
        for k in xrange(N):
            p = 2 * sp.pi * k * u / N
            z[0] += fx[k][0] * sp.cos(p) + fx[k][1] * sp.sin(p)
            z[1] += fx[k][1] * sp.cos(p) - fx[k][0] * sp.sin(p)
        result[u] = z[0] + z[1]*1j
        
    return result

Собственно, обратное преобразование выглядит похожим образом:

def idft_slow(fu):
    fu = np.asarray(fu, dtype=complex)
    N = fu.shape[0]
    result = np.zeros((N,2), dtype=float)
    for k in xrange(N):
        z = np.zeros(2)
        for u in xrange(N):
            p = 2 * sp.pi * k * u / N
            z[0] += fu[u].real * sp.cos(p) - fu[u].imag * sp.sin(p) 
            z[1] += fu[u].real * sp.sin(p) + fu[u].imag * sp.cos(p)
        result[k] = z
        
    return result / N

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

Всем удачи и до скорых встреч на нашем блоге!

Fast Fourier Transform Algorithm. Часть 1. Простая теория: 4 комментария

Добавить комментарий

Заполните поля или щелкните по значку, чтобы оставить свой комментарий:

Логотип WordPress.com

Для комментария используется ваша учётная запись WordPress.com. Выход /  Изменить )

Google photo

Для комментария используется ваша учётная запись Google. Выход /  Изменить )

Фотография Twitter

Для комментария используется ваша учётная запись Twitter. Выход /  Изменить )

Фотография Facebook

Для комментария используется ваша учётная запись Facebook. Выход /  Изменить )

Connecting to %s