Преобразование Фурье. The Fast and the Furious

Зачастую при разработке алгоритмов мы упираемся в предел вычислительной сложности, который, казалось бы, преодолеть невозможно. Преобразование Фурье имеет сложность $O(n^2)$, а быстрый вариант, предложенный около 1805 года Гаусом1 (и переизобретенный в 1965 году Джеймсом Кули и Джоном Тьюки) $O(nlog(n))$. В данной статье хочу вам показать, что можно получить результаты преобразования за линейное время $O(n)$ или даже достичь константной сложности $O(1)$ при определенных условиях, которые встречаются в реальных задачах.


Когда передо мной встала задача написания программы анализа передаточных функций звуковых систем в реальном времени, я, как и все, сперва обратился к быстрому преобразованию. Всё было хорошо, но при больших размерах временного окна загрузка CPU становилась неприлично большой и с этим надо было что-то делать. Было решено сделать паузу и изучить преобразование еще раз, а заодно поискать способы решения проблемы. Вернемся к оригинальному преобразованию Жозефа Фурье2:

$f(x) = \sum\limits_{-\infty}^{+\infty} c_ke^{2\pi ikx/T} \\ c_k = \frac{1}{T}\int\limits_0^Tf(x)e^{-2\pi ikx / T}dx$



Посмотрим внимательно, что же тут происходит. Каждое выходное значение в частотной области $c_k$ является суммой всех входных значений сигнала $f(x)$ умноженных на $e^{-2\pi ikx / T}$. Чтобы произвести вычисления нам надо для каждого выходного значения пройтись по всем входным данным, т.е. выполнить те самые $n^2$ операций.

Избавляемся от n


Напомню, что изначально была задача анализа звуковых данных в реальном времени. Для этого выбранное временное окно (по сути буфер) размером N заполняется данными с частотой fd соответствующей частоте дискретизации. С периодом Т происходят преобразования входных данных из временного окна в частотное. Если посмотреть на реальные числа, то N варьируется от 214 (16 384) до 216 (65 536) сэмплов (значения остались в наследство от БПФ, где размер окна обязан быть степенью двойки). Время T = 80мс (12,5 раз в секунду), что позволяет видеть достаточно удобно изменения и не слишком нагружать CPU и GPU. Частота дискретизации fd стандартная и составляет 48кГц. Давайте посчитаем, сколько же данных во временном окне изменяется между измерениями. За время Т в буфер поступает $48000 * 0,08 = 3840$ сэмплов. Таким образом в окне обновляется только от 5% до 23% данных. В худшем случае 95% (а в лучшем 73%, что тоже немало!) обрабатываемых сэмплов будут снова и снова попадать в преобразование несмотря на то, что уже были обработаны в прошлых итерациях.

Внимательный читатель в этот момент поднимет руку и скажет: «постойте, а как же коэффициент $e^{-2\pi ikx / T}$? Ведь при каждом новом преобразовании те же данные будут располагаться на новых позициях ряда, а как следствие иметь разные коэффициенты?» Каждому пятерку за внимательность, давайте вспоминать немаловажную деталь в преобразовании о которой часто забывают. При исследовании значений функции $f(t)$ на интервале от 0 до t, функция считается периодической, что позволяет безболезненно производить сдвиг функции влево или вправо во времени. Как следствие, мы имеем полное право не вставлять в конец новое значение и удалять из начала старое, а циклически заменять данные в буфере.

Для наглядности, можно записать в табличном виде как будет изменяться буфер:
t=0 f(0) f(1) f(2) f(3) f(4) f(5) f(6) f(7) f(8) f(9)
t=1 f(10) f(1) f(2) f(3) f(4) f(5) f(6) f(7) f(8) f(9)
t=2 f(10) f(11) f(2) f(3) f(4) f(5) f(6) f(7) f(8) f(9)
t=3 f(10) f(11) f(12) f(3) f(4) f(5) f(6) f(7) f(8) f(9)
t=4 f(10) f(11) f(12) f(13) f(4) f(5) f(6) f(7) f(8) f(9)

Можно записать как изменяется преобразование во времени от t1 до t2:

$F_t = F_{t-1} + \Delta F \\ \Delta F : \Delta c_k = \frac{1}{T}\int\limits_{t_1}^{t_2}(f_t(x) - f_{t-1}(x))e^{-2\pi ikx / T}dx$


Значение $F_{t-1}(x)$ является результатом предыдущего преобразования, а сложность вычисления $\Delta f(x)$ не зависит от размера временного окна и, следовательно, является константным. В итоге сложность преобразования будет $O(n)$*, т.к. всё что нам остается, это один раз пройтись по частотному окну и применить изменения для изменившихся за время Т сэмплов. Отдельно хочу обратить ваше внимание, что коэффициенты $e^{-2\pi ikx / T}$ могут быть посчитаны заранее, что дает дополнительный выигрыш в производительности, а в цикле останется только две операции: вычитание вещественных чисел и умножение вещественного числа на комплексное, на практике обе этих операции просты и дёшевы.

Для полноты картины осталось лишь обозначить начальное состояние, но тут всё просто:

$F_0(x) = 0$


*- конечно, конечная сложность всего преобразования так и останется $O(n^2)$, но оно будет выполнено постепенно, за n итераций, пока обновляется буфер. $O(n)$ — это сложность обновления данных, но именно это нам и нужно (при использовании БПФ сложность каждого преобразования $O(nlog(n))$).

А что, если копнуть глубже. Или избавляемся от второго n


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

А теперь давайте проанализируем результирующий набор данных, учитывая условия задачи. Мы имеем набор комплексных чисел, каждое из которых описывает амплитуду и фазу колебаний на определенной частоте. Частоту можно определить по формуле: $f[j] = j \frac{fd}{N}$ для $j < \frac{N}{2}$. Давайте оценим шаг частотного окна на наших данных: $\Delta f = \frac{fd}{N}$ Для N = 214: 2,93Гц (а для 216: 0,73Гц). Таким образом в диапазоне от 1кГц до 2кГц мы получим 341 результат. Попробуйте самостоятельно оценить сколько данных будет в диапазоне от 8кГц до 16кГц для N = 65536. Много, правда? Очень много! Нужно ли нам столько данных? Конечно в задачах отображения частотных характеристик звуковых систем, ответ: нет. А с другой стороны, для анализа в области низких частот, маленький шаг очень полезен. Не стоит забывать и о том, что впереди еще графика, которая должна эти объемы ($\frac{N}{2}$) преобразовать к понятному человеку виду (усреднение, сплайн или сглаживание) и вывести их на экран. Да и на высоких частотах даже имея 4К экран и выводя график в полноэкранном режиме при логарифмической оси частот размер шага быстро окажется намного меньше 1 пикселя.

Опытным путем можно выяснить, что вполне достаточно иметь лишь 48 точек на октаву, а чтобы иметь данные чуть более сглаженными и усредненными, предлагаю остановиться на значении 96. В звуковом диапазоне частот от 20Гц до 20кГц нетрудно насчитать всего 10 октав: 20, 40, 80, 160, 320, 640, 1280, 2560, 5120, 10240, 20480, каждую из которых можно разделить на заданное количество поддиапазонов (не забывайте, что разбиение следует производить геометрически, а не арифметически), следовательно, более чем достаточно произвести преобразование только для 960 частот, чтобы получить результат, что в 16...65 раз меньше изначального варианта.

Таким образом, комбинируя оба подхода, получаем константную сложность выполнения алгоритма обновления данных $O(1)$.

Мёд в квадрате и ложка дёгтя


Вот теперь можно смело заявить, что от сложности $O(n^2)$ мы пришли к сложности $O(1)$ использовав два простых лайфхака:

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

Но, конечно, жизнь была бы действительно сказкой, если бы не одно но. Применение этих двух подходов позволило действительно разгрузить CPU так, что догадаться о том, что он рассчитывает преобразование Фурье и выводит результаты на экран даже при $N=2^{20}$ было сложно. Но кара не заставила себя ждать, когда ваши сигналы в реальности не периодичны (а это обязательно для получения верных результатов преобразования) и подобрать подходящий размер окна не представляется возможным, возникает необходимость использования различных оконных функций, что уже не позволяет использовать полноценно первый шаг. Практика показала, что применение оконных функций критично при исследовании сигналов с частотой меньше, чем $0,1f_d$. На больших частотах количество периодов, попавших во временное окно, значительно и ослабляет искажения, возникающие в следствии наличия разрыва первого рода (между f(0) и f(N-1)) в изначальной функции.

От второго шага я тоже в итоге отказался и вернулся обратно к БПФ, т.к. выигрыш в данной задаче уже был невелик.

В заключение


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

Увы, для меня в данной задаче, это осталось лишь небольшим математическим развлечением, но я надеюсь, что оно вдохновит вас заняться на праздниках исследованием других алгоритмов с точки зрения изменений входных данных во времени :)

Литература


  1. Cooley–Tukey FFT algorithm
  2. Е.А.Власова Ряды. Издательство МГТУ им Н.Э.Баумана. Москва. 2002

Изображение взято из манги Митио Сибуя «ЗАНИМАТЕЛЬНАЯ МАТЕМАТИКА. АНАЛИЗ ФУРЬЕ»
Поделиться публикацией
AdBlock похитил этот баннер, но баннеры не зубы — отрастут

Подробнее
Реклама

Комментарии 16

    0
    На моём ноутбуке FFTW на преобразование 65536 семплов тратит чуть более одной миллисекунды. Может, стоило оптимизировать сам алгоритм преобразования?
      0
      Этим я в итоге и занялся, но это уже совсем другая история, не столь интересная.
      0
      достаточно произвести преобразование только для 960 частот
      Если я правильно понял условия задачи — нет, недостаточно, так как нужно делать не преобразование 960 частот, а преобразование 960 частотных полос. А иначе получится не усреднение, а прореживание — и такой прореженный спектр будет сильно зашумлен.

      Делают так:
      1) выполняется прямое преобразование Фурье,
      2) обнуляется отрицательная часть частот, чтобы привести сигнал к аналитическому виду — в котором мнимая часть повёрнута на 90° относительно исходного,
      3) спектр разбивается на диапазоны в логарифмическом масштабе с использованием оконных функций,
      4) над каждым диапазоном выполняется обратное преобразование Фурье, причём минимального размера (если в диапазон попало 8 частот, то и размер FFT также будет 8),
      5) с получившего комплексного сигнала считается модуль (квадратный корень суммы квадратов действительной и мнимой части) — это и будет амплитуда.
        0
        Если я правильно понял условия задачи

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

        Делают так

        Спасибо!
          0
          Должен уточнить, что в предыдущем сообщении я пропустил слово «например». Существует несколько подходов к решению таких задач, в том числе и без использования преобразования Фурье вообще. Я описал вариант, предполагающий возможность обратного преобразования.
          0
          4) над каждым диапазоном выполняется обратное преобразование Фурье, причём минимального размера (если в диапазон попало 8 частот, то и размер FFT также будет 8),
          5) с получившего комплексного сигнала считается модуль (квадратный корень суммы квадратов действительной и мнимой части)

          Подскажите, пожалуйста, почему возникает необходимость делать iFFT, разве нет возможности вычислить это значение из частотных компонент?
            0
            Потому что когда мы говорим об амплитуде, мы говорим об огибающей сигнала. Если этот сигнал — синусоида, то его огибающая константна на протяжении всего времени, поэтому её и можно получить напрямую из спектра:



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



            А преобразования Фурье туда-сюда нужны, чтобы выделить интересуемый частотный диапазон и привести сигнал к комплексному виду:
          +1
          Можно ссылку на мангу? (И статью оформленную в таком стиле)
            0

            Честно пытался её прочесть справа налево

              0
              в первом ряду картинок даже получается, потом приходит осознание: «что-то тут не так...»
            +1
            За исключением картинки статья никак не связана с БПФ. Более того в ней делаются странные утверждения. И еще в БПФ количество точек не обязательно должно быть кратно степени двойки, достаточно что бы раскладываться на множители (и желательно малые 2,3,5...)
              0
              Мне кажется, сложность FFT можно сильно уменьшить в случае, когда предполагается представление результатов с логарифмической осью частот. Я попытался загуглить, но не смог найти по этой теме никакой информации.

              Как правильно замечено в статье, на частотах порядка 10Гц и 10кГц совершенно разные минимальные требования по частотному разрешению, в то время как FFT дает одинаковое разрешение во всем диапазоне частот. И в итоге, после применения FFT большого размера, в высоких частотах происходит обычное сглаживание — т.е., мы по сути выкидываем ощутимую часть результатов.

              Моя идея: уменьшить размер настолько, чтобы его было достаточно для высоких частот. А низкие — считать отдельно.
              Ведь от чего зависит частотное разрешение FFT? От размера и от частоты сэмплирования. И увеличить разрешение можно как увеличение размера, так и уменьшением частоты сэмплирования.
              При увеличении размера мы получаем возрастание сложности алгоритма и ухудшение временного разрешения. При уменьшении частоты сэмплирования мы получаем более узкую полосу частот для анализа — а это нам вполне подходит, если мы отдельно анализируем высокие и низкие частоты!

              То есть: вместо того, чтобы, например, сделать FFT размера 4096, мы можем сделать FFT размера 1024 3 раза: один раз с изначальными данными, второй раз с вдвое меньшей частотой сэмплирования, третий раз в четверо меньшей частотой относительно начальной (или с вдвое меньшей относительно второго преобразования). В итоге получатся результаты для участков частот [0, Fs/2], [0, Fs/4], [0, Fs/8], и нам потребуется только правильно смешать их, чтобы получить примерно одинаковую точность на всех полосах. Например, от 10 до 20 Гц брать третье преобразование, от 20 до 40 брать второе, от 40 до 80 брать первое.
              Эффективное разрешение по частоте получается «size * (2 ^ iters_count)», т.к. каждая итерация увеличивает точность в два раза для всех частот вплоть до уменьшенной вдвое (относительно предыдущей итерации) частоты Найквиста.

              Этот алгоритм имеет целых два плюса:
              1. Временное разрешение для высоких частот улучшается.
              2. Сложность падает, т.к. для того же разрешения в низких частотах достаточно FFT размера в 2^iters_count меньше, чем при обычном преобразовании. При этом сложностью за счет нескольких преобразований можно пренебречь, т.к. сумма геометрической прогрессии 1/2^n, т.е. сложность от количества преобразований увеличивается не более, чем в два раза (ровно в 2 при бесконечном числе итераций).
              Из минусов — необходимо хранить сигнал iters_count раз, что увеличивает потребление памяти.
              Для получения сигнала с меньшей частотой сэмплирования достаточно применить фильтр низких частот вида «среднее двух значений», и его можно применять на данных предыдущей итерации, так что это почти не замедляет работу.

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

              P.S. Никто не подскажет, какая нормировка результата FFT правильная?
              Если я делаю преобразование для белого шума, то для разных размеров сумма и максимум магнитуды должны оставаться примерно одинаковыми?
              А для преобразования чистого синуса — сумма и максимум магнитуды должны оставаться одинаковыми?
                0
                Описанная вами идея действительно имеет место быть и распространена. Есть два основных подхода: изменять размер временного окна или частоту дискретизации. Встречается данный алгоритм под названием multitimewindow (MTW FFT).
                  0
                  Мне кажется, сложность FFT можно сильно уменьшить в случае, когда предполагается представление результатов с логарифмической осью частот. Я попытался загуглить, но не смог найти по этой теме никакой информации.
                  Попробуйте погуглить "Constant-Q transform".
                    0
                    Нормировка может быть разной. Но если хочется симметрии, что бы нормы сигнала совпадали то нормируют на 1/sqrt(n)
                    y[k]=C*sum(x[j]*Wn^(j*k),j=0..n-1)
                    x[j]=C*sum(y[k]*Wn^(-k*j),k=0..n-1)
                    C=1/sqrt(n)
                    Wn=exp(-2*pi*i/n)

                      0

                      Один из способов вычисления FFT разделяет семплы на четные и нечетные, рекурсивно вычисляет для них FFT и объединяет результаты.
                      Поэтому ваше решение может не сработать.

                    Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.

                    Самое читаемое