Быстрое преобразование Фурье (FFT): туда и обратно

Код выполняет прямое быстрое преобразование Фурье (метод fft) и обратное (ifft), на простом тесте всё сошлось в пределах погрешности в 16-17 знаке после запятой (элементы тестового массива data выводятся "как есть").

В английской wiki есть псевдокод алгоритма, в русской - описание задачи.

Программа написана на консольном C++ и проверена в актуальной сборке Visual Studio 2019 (код единственного файла .cpp можно вставить в пустой проект). Учитывая, что комплексные числа теперь - часть стандарта, конечно же, написания отдельных классов для них не требуется.

#include <complex>
#include <iostream>
#include <valarray>
#define _USE_MATH_DEFINES
#include <math.h>

typedef std::complex <double> Complex;
typedef std::valarray <Complex> CArray;

void fft(CArray& x) {
 const size_t N = x.size();
 if (N <= 1) return;
 CArray even = x[std::slice(0, N / 2, 2)];
 CArray  odd = x[std::slice(1, N / 2, 2)];
 fft(even);
 fft(odd);
 for (size_t k = 0; k < N / 2; ++k) {
  Complex t = std::polar(1.0, -2 * M_PI * k / N) * odd[k];
  x[k] = even[k] + t;
  x[k + N / 2] = even[k] - t;
 }
}

void ifft(CArray& x) {
 x = x.apply(std::conj);
 fft(x);
 x = x.apply(std::conj);
 x /= x.size();
}

int main() {
 const Complex test[] = { 0, 1, 2, 3, 4, 5, 6, 7 };
 CArray data(test, 8);
 fft(data);
 std::cout << "fft" << std::endl;
 for (int i = 0; i < 8; ++i) {
  std::cout << data[i] << std::endl;
 }
 ifft(data);
 std::cout << std::endl << "ifft" << std::endl;
 for (int i = 0; i < 8; ++i) {
  std::cout << data[i] << std::endl;
 }
 return 0;
}

Функция fft работает по алгоритму Кули — Тьюки, ниже представлена альтернативная версия функции, теоретически более быстрая и точная, но практически менее понятная :)

void fft(CArray& x) {
 unsigned int N = x.size(), k = N, n;
 double thetaT = M_PI / N;
 Complex phiT = Complex(cos(thetaT), -sin(thetaT)), T;
 while (k > 1) {
  n = k;
  k >>= 1;
  phiT = phiT * phiT;
  T = 1.0L;
  for (unsigned int l = 0; l < k; l++) {
   for (unsigned int a = l; a < N; a += n) {
    unsigned int b = a + k;
    Complex t = x[a] - x[b];
    x[a] += x[b];
    x[b] = t * T;
   }
   T *= phiT;
  }
 }
 unsigned int m = (unsigned int)log2(N);
 for (unsigned int a = 0; a < N; a++) {
  unsigned int b = a;
  b = (((b & 0xaaaaaaaa) >> 1) | ((b & 0x55555555) << 1));
  b = (((b & 0xcccccccc) >> 2) | ((b & 0x33333333) << 2));
  b = (((b & 0xf0f0f0f0) >> 4) | ((b & 0x0f0f0f0f) << 4));
  b = (((b & 0xff00ff00) >> 8) | ((b & 0x00ff00ff) << 8));
  b = ((b >> 16) | (b << 16)) >> (32 - m);
  if (b > a) {
   Complex t = x[a];
   x[a] = x[b];
   x[b] = t;
  }
 }
}

26.12.2021, 00:33 [1760 просмотров]