11.17. Вычисление быстрого преобразования Фурье

11.17. Вычисление быстрого преобразования Фурье

Проблема

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

Решение

Программный код примера 11.33 обеспечивает базовую реализацию БПФ.

Пример 11.33. Реализация БПФ

#include <iostream>

#include <complex>

#include <cmath>

#include <iterator>

using namespace std;

unsigned int bitReverse(unsigned int x, int log2n) {

 int n = 0;

 int mask = 0x1;

 for (int i=0; i < log2n; i++) {

  n <<= 1;

  n |= (x & 1);

  x >>= 1;

 }

 return n;

}

const double PI = 3.1415926536;

template<class Iter_T>

void fft(Iter_r a, Iter_r b, int log2n) {

 typedef typename iterator_traits<Iter T>::value_type complex;

 const complex J(0, 1);

 int n = 1 << log2n;

 for (unsigned int i=0; i < n; ++i) {

  b[bitReverse(i, log2n)] = a[i];

 }

 for (int s = 1; s <= log2n; ++s) {

  int m = 1 << s;

  int m2 = m >> 1;

  complex w(1, 0);

  complex wm = exp(-J * (PI / m2));

  for (int j=0; j < m2; ++j) {

   for (int k=j; k < n; k += m) {

    complex t = w * b[k + m2];

    complex u = b[k];

    b[k] = u + t;

    b[k + m2] = u - t;

   }

   w *= wm;

  }

 }

}

int main() {

 typedef complex<double> cx;

 cx a[] = { cx(0, 0), cx(1, 1), cx(3, 3), cx(4, 4),

  cx(4, 4), cx(3, 3), cx(1, 1), cx(0, 0) };

 cx b[8];

 fft(a, b, 3);

 for (int i=0; i<8; ++i) cout << b[i] << " ";

}

Программа примера 11.33 выдает следующий результат.

(16,16)

(-4.82843,-11.6569)

(0,0)

(-0.343146,0.828427)

(0.0)

(0.828427,-0.343146)

(0,0)

(-11.6569,-4.82843)

Обсуждение

Преобразование Фурье играет важную роль в спектральном анализе и часто используется в технических и научных приложениях. БПФ — это алгоритм вычисления ДПФ, который имеет сложность порядка N log₂(N) в отличие от ожидаемой сложности N? для простой реализации ДПФ. Такое впечатляющее ускорение достигается в БПФ благодаря устранению избыточных вычислений.

Очень не просто найти хорошую реализацию БПФ, написанную на «правильном» C++ (т. е. когда программа на C++ не является механическим переложением алгоритмов, написанных на Фортране или С) и которая не была бы защищена сильно ограничивающей лицензией. Представленный в примере 11.33 программный код основан на открытом коде, который можно найти в сетевой конференции Usenet, посвященной цифровой обработке сигналов (comp.dsp). Большим преимуществом реализации БПФ на правильном C++ по сравнению с более распространенным решением в стиле С является то, что стандартная библиотека содержит шаблон complex, который позволяет существенно снизить объем необходимого программного кода. В представленной в примере 11.33 функции fft() основное внимание уделялось простоте, а не эффективности.