Быстрое преобразование Фурье

Быстрое преобразование Фурье (БПФ, FFT) — это алгоритм быстрого вычисления дискретного преобразования Фурье (ДПФ). То есть, алгоритм вычисления за количество действий, меньшее чем $O(N^2)$, требуемых для прямого (по формуле) вычисления ДПФ. Иногда под БПФ понимается один из быстрых алгоритмов, называемый алгоритмом прореживания по частоте/времени или алгоритмом по основанию 2, имеющего сложность $O(N\log(N))$.

Основной алгоритм

Покажем как выполнить дискретное преобразование Фурье за $O(N(p_1+\cdots+p_n))$ действий при $N=p_1p_2\cdots p_n$. В частности, при $N=2^n$ понадобится $O(N\log(N))$ действий.

Дискретное преобразование Фурье преобразует набор чисел $a_0, \dots, a_{n-1}$ в набор чисел $b_0, \dots, b_{n-1}$, такой, что $b_i=\sum_{j=0}^{n-1}a_j\varepsilon^{ij}$, где $\varepsilon^n=1$ и $\varepsilon^k\neq 1$ при $0<k<n$. Алгоритм быстрого преобразования Фурье применим к любым коммутативным ассоциативным кольцам с единицей. Чаще всего этот алгоритм применяют к полю комплексных чисел (c $\varepsilon=e^{2\pi i/n}$) и к кольцам вычетов.

Основной шаг алгоритма состоит в сведении задачи для $N$ чисел к задаче для $p=N/q$ числам, где $q$ — делитель $N$. Пусть мы уже умеем решать задачу для $N/q$ чисел. Применим преобразование Фурье к наборам $a_i,a_{q+i}, \dots, a_{q(p-1)+i}$ для $i=0,1,\dots,q-1$. Покажем теперь, как за $O(Np)$ действий решить исходную задачу. Заметим, что $b_i=\sum_{j=0}^{q-1} \varepsilon^{ij} (\sum_{k=0}^{p-1}a_{kq+j}\varepsilon^{kiq})$. Выражения в скобках нам уже известны — это $i\pmod p$-тое число после преобразования Фурье $j$-той группы. Таким образом, для вычисления каждого $b_i$ нужно $O(q)$ действий, а для вычисления всех $b_i$ — $O(Nq)$ действий, что и требовалось получить.

Обратное преобразование Фурье

Для обратного преобразования Фурье можно применять алгоритм прямого преобразования Фурье — нужно лишь использовать $\varepsilon^{-1}$ вместо $\varepsilon$ (или применить операцию комплексного сопряжения в начале к входным данным, а затем к результату, полученному после прямого преобразования Фурье) и окончательный результат поделить на $N$.

Общий случай

Общий случай может быть сведён к предыдущему. Пусть $4N>2^k\ge2N$. Заметим, что $b_i=\varepsilon^{-i^2/2}\sum_{j=0}^{N-1}\varepsilon^{(i+j)^2/2}\varepsilon^{-j^2/2}a_j$. Обозначим $\bar{a}_i=\varepsilon^{-i^2/2}a_i, \bar{b}_i=\varepsilon^{i^2/2}b_i, c_i=\varepsilon^{(2N-2-i)^2/2}$. Тогда $\bar{b}_i=\sum_{j=0}^{2N-2-i}\bar{a}_jc_{2N-2-i-j}$, если положить $\bar{a}_i=0$ при $i\ge N$.

Таким образом задача сведена к вычислению свёртки, но это можно сделать с помощью трёх преобразований Фурье для $2^k$ элементов. Выполняем прямое преобразование Фурье для $\{\bar{a}_i\}_{i=0}^{i=2^k-1}$ и $\{c_i\}_{i=0}^{i=2^k-1}$, перемножаем поэлементно результаты и выполняем обратное преобразование Фурье.

Вычисления всех $\bar{a}_i$ и $c_i$ требуют $O(N)$ действий, три преобразования Фурье требуют $O(N\log(N))$ действий, перемножение результатов преобразований Фурье требует $O(N)$ действий, вычисление всех $b_i$ зная значения свертки требует $O(N)$ действий. Итого для дискретного преобразования Фурье требуется $O(N\log(N))$ действий для любого $N$.

Этот алгоритм быстрого преобразования Фурье может работать над кольцом только когда известны первообразные корни из единицы степеней $2N$ и $2^k$.

Вывод преобразования из ДПФ

Дискретное преобразование Фурье для вектора $\vec x$, состоящего из N элементов, имеет вид:

$\vec X = \hat A \vec x$

элементы матрицы $\hat A$ имеют вид: $a_{N}^{mn} = \exp \left( -2\pi i \frac{mn}{N} \right)$.

Пусть N четно, тогда ДПФ можно переписать следующим образом:

$X_m=\sum_{n=0}^{N-1} x_n a_{N}^{mn} = \sum_{n=0}^{N/2-1}x_{2n} a_{N}^{2nm} + \sum_{n=0}^{N/2-1}x_{2n+1} a_{N}^{(2n+1)m}$

Коэффициенты $a_{N}^{2nm}$ и $a_{N}^{(2n+1)m}$ можно переписать следующим образом (M=N/2):

$a_{N}^{2nm} = \exp \left( -2\pi i \frac{2mn}{N} \right) = \exp \left( -2\pi i \frac{mn}{N/2} \right) = a_{M}^{nm}$

$a_{N}^{(2n+1)m} = \exp \left( -2\pi i \frac{m}{N} \right) a_{M}^{nm}$

В результате получаем:

$X_m=\sum_{n=0}^{M-1} x_{2n} a_{M}^{nm} + \exp \left( -2\pi i \frac{m}{N} \right) \sum_{n=0}^{M-1} x_{2n+1} a_{M}^{nm}$

То есть дискретное преобразование Фурье от вектора, состоящего из N отсчетов, свелось к линейной композиции двух ДПФ от $\frac N2$ отсчетов, и если для первоначальной задачи требовалось $N^2$ операций, то для полученной композиции — $\frac{N^2}{2}$. Если M является степенью двух, то это разделение можно продолжать рекурсивно до тех пор, пока не дойдем до двухточечного преобразования Фурье, которое вычисляется по следующим формулам:

$\begin{cases} X_0=x_0+x_1\\ X_1=x_0-x_1 \end{cases}$

Пример программы

Ниже приведен пример расчета комплексного БПФ, написанный на С:

C  

//_________________________________________________________________________________________
//_________________________________________________________________________________________
//
// NAME:          FFT.
// PURPOSE:       Быстрое преобразование Фурье: Комплексный сигнал в комплексный спектр и обратно.
//                В случае действительного сигнала в мнимую часть (Idat) записываются нули.
//                Количество отсчетов - кратно 2**К - т.е. 2, 4, 8, 16, ... (см. комментарии ниже).
//                (C) Sergey Aleynik.  [e-mail]
//              
// PARAMETERS:  
//
//    float *Rdat    [in, out] - Real part of Input and Output Data (Signal or Spectrum)
//    float *Idat    [in, out] - Imaginary part of Input and Output Data (Signal or Spectrum)
//    int    N       [in]      - Input and Output Data length (Number of samples in arrays)
//    int    LogN    [in]      - Logarithm2(N)
//    int    Ft_Flag [in]      - Ft_Flag = FT_ERR_DIRECT  (i.e. -1) - Direct  FFT  (Signal to Spectrum)
//                               Ft_Flag = FT_ERR_INVERSE (i.e.  1) - Inverse FFT  (Spectrum to Signal)
//
// RETURN VALUE:  false on parameter error, true on success.
//_________________________________________________________________________________________
//
// NOTE: In this algorithm N and LogN can be only:
//       N    = 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096, 8192, 16384;
//       LogN = 2, 3,  4,  5,  6,   7,   8,   9,   10,   11,   12,   13,    14;
//_________________________________________________________________________________________
//_________________________________________________________________________________________
 
 
#define  NUMBER_IS_2_POW_K(x)   ((!((x)&((x)-1)))&&((x)>1))  // x is pow(2, k), k=1,2, ...
#define  FT_DIRECT        -1    // Direct transform.
#define  FT_INVERSE        1    // Inverse transform.
 
 
bool  FFT(float *Rdat, float *Idat, int N, int LogN, int Ft_Flag)
{
  // parameters error check:
  if((Rdat == NULL) || (Idat == NULL))                  return false;
  if((N > 16384) || (N < 1))                            return false;
  if(!NUMBER_IS_2_POW_K(N))                             return false;
  if((LogN < 2) || (LogN > 14))                         return false;
  if((Ft_Flag != FT_DIRECT) && (Ft_Flag != FT_INVERSE)) return false;
 
  register int  i, j, n, k, io, ie, in, nn;
  float         ru, iu, rtp, itp, rtq, itq, rw, iw, sr;
 
  static const float Rcoef[14] =
  {  -1.0000000000000000F,  0.0000000000000000F,  0.7071067811865475F,
      0.9238795325112867F,  0.9807852804032304F,  0.9951847266721969F,
      0.9987954562051724F,  0.9996988186962042F,  0.9999247018391445F,
      0.9999811752826011F,  0.9999952938095761F,  0.9999988234517018F,
      0.9999997058628822F,  0.9999999264657178F
  };
  static const float Icoef[14] =
  {   0.0000000000000000F, -1.0000000000000000F, -0.7071067811865474F,
     -0.3826834323650897F, -0.1950903220161282F, -0.0980171403295606F,
     -0.0490676743274180F, -0.0245412285229122F, -0.0122715382857199F,
     -0.0061358846491544F, -0.0030679567629659F, -0.0015339801862847F,
     -0.0007669903187427F, -0.0003834951875714F
  };
 
  nn = N >> 1;
  ie = N;
  for(n=1; n<=LogN; n++)
  {
    rw = Rcoef[LogN - n];
    iw = Icoef[LogN - n];
    if(Ft_Flag == FT_INVERSE) iw = -iw;
    in = ie >> 1;
    ru = 1.0F;
    iu = 0.0F;
    for(j=0; j<in; j++)
    {
      for(i=j; i<N; i+=ie)
      {
        io       = i + in;
        rtp      = Rdat[i]  + Rdat[io];
        itp      = Idat[i]  + Idat[io];
        rtq      = Rdat[i]  - Rdat[io];
        itq      = Idat[i]  - Idat[io];
        Rdat[io] = rtq * ru - itq * iu;
        Idat[io] = itq * ru + rtq * iu;
        Rdat[i]  = rtp;
        Idat[i]  = itp;
      }
 
      sr = ru;
      ru = ru * rw - iu * iw;
      iu = iu * rw + sr * iw;
    }
 
    ie >>= 1;
  }
 
  for(j=i=1; i<N; i++)
  {
    if(i < j)
    {
      io       = i - 1;
      in       = j - 1;
      rtp      = Rdat[in];
      itp      = Idat[in];
      Rdat[in] = Rdat[io];
      Idat[in] = Idat[io];
      Rdat[io] = rtp;
      Idat[io] = itp;
    }
 
    k = nn;
 
    while(k < j)
    {
      j   = j - k;
      k >>= 1;
    }
 
    j = j + k;
  }
 
  if(Ft_Flag == FT_DIRECT) return true;
 
  rw = 1.0F / N;
 
  for(i=0; i<N; i++)
  {
    Rdat[i] *= rw;
    Idat[i] *= rw;
  }
 
  return true;
}
 
 
// Пример вычисления БПФ от одного периода косинусного
// действительного сигнала
 
void Test_FFT()
{
  static float Re[8];
  static float Im[8];
  float  p = 2 * 3.141592653589 / 8; // будет 8 отсчетов на период
 
  int i;
  // формируем сигнал
  for(i=0; i<8; i++)
  {
    Re[i] = cos(p * i);  // заполняем действительную часть сигнала
    Im[i] = 0.0;         // заполняем мнимую часть сигнала
  }
 
  FFT(Re, Im, 8, 3, -1); // вычисляем прямое БПФ
 
  // выводим действительную и мнимую части спектра и спектр мощности
  FILE *f = fopen("spectrum.txt", "w");
  for(i=0; i<8; i++)
  {
    fprintf(f, "%10.6f  %10.6f  %10.6f\n", Re[i], Im[i], Re[i]*Re[i]+Im[i]*Im[i]);
  }
  fclose(f);
}

Ниже приведен пример вычисления модуля спектра действительного массива чисел на основе реализации быстрого преобразования Фурье, написанный на C++ :

С++  

// AVal - массив анализируемых данных, Nvl - длина массива должна быть кратна степени 2.
// FTvl - массив полученных значений, Nft - длина массива должна быть равна Nvl.
 
const double TwoPi = 6.283185307179586;
 
void FFTAnalysis(double *AVal, double *FTvl, int Nvl, int Nft) {
  int i, j, n, m, Mmax, Istp;
  double Tmpr, Tmpi, Wtmp, Theta;
  double Wpr, Wpi, Wr, Wi;
  double *Tmvl;
 
  n = Nvl * 2; Tmvl = new double[n];
 
  for (i = 0; i < Nvl; i++) {
    j = i * 2; Tmvl[j] = 0; Tmvl[j+1] = AVal[i];
  }
 
  i = 1; j = 1;
  while (i < n) {
    if (j > i) {
      Tmpr = Tmvl[i]; Tmvl[i] = Tmvl[j]; Tmvl[j] = Tmpr;
      Tmpr = Tmvl[i+1]; Tmvl[i+1] = Tmvl[j+1]; Tmvl[j+1] = Tmpr;
    }
    i = i + 2; m = Nvl;
    while ((m >= 2) && (j > m)) {
      j = j - m; m = m >> 1;
    }
    j = j + m;
  }
 
  Mmax = 2;
  while (n > Mmax) {
    Theta = -TwoPi / Mmax; Wpi = Sin(Theta);
    Wtmp = Sin(Theta / 2); Wpr = Wtmp * Wtmp * 2;
    Istp = Mmax * 2; Wr = 1; Wi = 0; m = 1;
 
    while (m < Mmax) {
      i = m; m = m + 2; Tmpr = Wr; Tmpi = Wi;
      Wr = Wr - Tmpr * Wpr - Tmpi * Wpi;
      Wi = Wi + Tmpr * Wpi - Tmpi * Wpr;
 
      while (i < n) {
        j = i + Mmax;
        Tmpr = Wr * Tmvl[j] - Wi * Tmvl[j-1];
        Tmpi = Wi * Tmvl[j] + Wr * Tmvl[j-1];
 
        Tmvl[j] = Tmvl[i] - Tmpr; Tmvl[j-1] = Tmvl[i-1] - Tmpi;
        Tmvl[i] = Tmvl[i] + Tmpr; Tmvl[i-1] = Tmvl[i-1] + Tmpi;
        i = i + Istp;
      }
    }
 
    Mmax = Istp;
  }
 
  for (i = 0; i < Nft; i++) {
    j = i * 2; FTvl[Nft - i - 1] = Sqrt(Sqr(Tmvl[j]) + Sqr(Tmvl[j+1]));
  }
 
  delete []Tmvl;
}

Пример реализации на Delphi :

Delphi  

// AVal - массив анализируемых данных, Nvl - длина массива, должна быть кратна степени 2.
// FTvl - массив полученных значений, Nft - длина массива, должна быть равна Nvl / 2 или меньше.
 
type
  TArrayValues = array of Double;
 
const
  TwoPi = 6.283185307179586;
 
procedure FFTAnalysis(var AVal, FTvl: TArrayValues; Nvl, Nft: Integer);
var
  i, j, n, m, Mmax, Istp: Integer;
  Tmpr, Tmpi, Wtmp, Theta: Double;
  Wpr, Wpi, Wr, Wi: Double;
  Tmvl: TArrayValues;
begin
  n:= Nvl * 2; SetLength(Tmvl, n);
 
  for i:= 0 to Nvl-1 do begin
    j:= i * 2; Tmvl[j]:= 0; Tmvl[j+1]:= AVal[i];
  end;
 
  i:= 1; j:= 1;
  while i < n do begin
    if j > i then begin
      Tmpr:= Tmvl[i]; Tmvl[i]:= Tmvl[j]; Tmvl[j]:= Tmpr;
      Tmpr:= Tmvl[i+1]; Tmvl[i+1]:= Tmvl[j+1]; Tmvl[j+1]:= Tmpr;
    end;
    i:= i + 2; m:= Nvl;
    while (m >= 2) and (j > m) do begin
      j:= j - m; m:= m div 2;
    end;
    j:= j + m;
  end;
 
  Mmax:= 2;
  while n > Mmax do begin
    Theta:= -TwoPi / Mmax; Wpi:= Sin(Theta);
    Wtmp:= Sin(Theta / 2); Wpr:= Wtmp * Wtmp * 2;
    Istp:= Mmax * 2; Wr:= 1; Wi:= 0; m:= 1;
 
    while m < Mmax do begin
      i:= m; m:= m + 2; Tmpr:= Wr; Tmpi:= Wi;
      Wr:= Wr - Tmpr * Wpr - Tmpi * Wpi;
      Wi:= Wi + Tmpr * Wpi - Tmpi * Wpr;
 
      while i < n do begin
        j:= i + Mmax;
        Tmpr:= Wr * Tmvl[j] - Wi * Tmvl[j-1];
        Tmpi:= Wi * Tmvl[j] + Wr * Tmvl[j-1];
 
        Tmvl[j]:= Tmvl[i] - Tmpr; Tmvl[j-1]:= Tmvl[i-1] - Tmpi;
        Tmvl[i]:= Tmvl[i] + Tmpr; Tmvl[i-1]:= Tmvl[i-1] + Tmpi;
        i:= i + Istp;
      end;
    end;
 
    Mmax:= Istp;
  end;
 
  for i:= 1 to Nft-1 do begin
    j:= i * 2; FTvl[Nft - i - 1]:= Sqrt(Sqr(Tmvl[j]) + Sqr(Tmvl[j+1]));
  end;
 
  SetLength(Tmvl, 0);
end;

Пример реализации на C#

С#  

using System;
using System.Numerics;
 
namespace FFT
{
    public class FFT
    {
        /// <summary>
        /// Вычисление поворачивающего модуля e^(-i*2*PI*k/N)
        /// </summary>
        /// <param name="k"></param>
        /// <param name="N"></param>
        /// <returns></returns>
        private static Complex w(int k, int N)
        {
            if (k % N == 0) return 1;
            double arg = -2 * Math.PI * k / N;
            return new Complex(Math.Cos(arg), Math.Sin(arg));
        }
        /// <summary>
        /// Возвращает спектр сигнала
        /// </summary>
        /// <param name="x">Массив значений сигнала. Количество значений должно быть степенью 2</param>
        /// <returns>Массив со значениями спектра сигнала</returns>
        public static Complex[] fft(Complex[] x)
        {
            Complex[] X;
            int N = x.Length;
            if (N == 2)
            {
                X = new Complex[2];
                X[0] = x[0] + x[1];
                X[1] = x[0] - x[1];
            }
            else
            {
                Complex[] x_even = new Complex[N / 2];
                Complex[] x_odd = new Complex[N / 2];
                for (int i = 0; i < N / 2; i++)
                {
                    x_even[i] = x[2 * i];
                    x_odd[i] = x[2 * i + 1];
                }
                Complex[] X_even = fft(x_even);
                Complex[] X_odd = fft(x_odd);
                X = new Complex[N];
                for (int i = 0; i < N / 2; i++)
                {
                    X[i] = X_even[i] + w(i, N) * X_odd[i];
                    X[i + N / 2] = X_even[i] - w(i, N) * X_odd[i];
                }
            }
            return X;
        }
        /// <summary>
        /// Центровка массива значений полученных в fft (спектральная составляющая при нулевой частоте будет в центре массива)
        /// </summary>
        /// <param name="X">Массив значений полученный в fft</param>
        /// <returns></returns>
        public static Complex[] nfft(Complex[] X)
        {
            int N = X.Length;
            Complex[] X_n = new Complex[N];
            for (int i = 0; i < N / 2; i++)
            {
                X_n[i] = X[N / 2 + i];
                X_n[N / 2 + i] = X[i];
            }
            return X_n;
        }
    }
}

Графическое представление работы вышеприведенного алгоритма.

Ссылки

• Быстрое преобразование Фурье и его приложения — Описан сам метод, его основные свойства. Приведены реализации алгоритма для проведения этого преобразования на языках программирования C#, C++, Delphi, Visual Basic 6, Zonnon.

• Преобразование Фурье — Суть метода, теоремы, физические эффекты при преобразовании реальных сигналов, примеры оптимизированных программ на C++.

• Описание и реализация Быстрого преобразования Фурье на сайте e-maxx.ru

• Реализация преобразования Фурье на Delphi — Проект приложения Delphi с исходными кодами, реализующего преобразование сигнала, представленного в виде квадратурных составляющих I и Q.

• БПФ по основанию 2 с прореживанием по времени  (рус.). Архивировано из первоисточника 5 февраля 2012. Проверено 15 ноября 2010.

• БПФ по основанию 2 с прореживанием по частоте  (рус.). Архивировано из первоисточника 5 февраля 2012. Проверено 15 ноября 2010.

• Полифазное БПФ  (рус.). Архивировано из первоисточника 5 февраля 2012. Проверено 15 ноября 2010.

• FFTW — свободно распространяемая библиотека, написанная на С, поддерживает разложение n-мерного сигнала, представленного как вещественными, так и комплексными числами.  (англ.)