Быстрое преобразование Фурье | это... Что такое Быстрое преобразование Фурье? (original) (raw)

Быстрое преобразование Фурье (БПФ, 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}$ ) и к кольцам вычетов.

Основной шаг алгоритма состоит в сведении задачи для чисел к задаче для p=N/q числам, где — делитель . Пусть мы уже умеем решать задачу для 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$ -тое число после преобразования Фурье -той группы. Таким образом, для вычисления каждого b_i нужно O(q) действий, а для вычисления всех b_i — O(Nq) действий, что и требовалось получить.

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

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

Общий случай

Общий случай может быть сведён к предыдущему. Пусть $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))$ действий для любого .

Этот алгоритм быстрого преобразования Фурье может работать над кольцом только когда известны первообразные корни из единицы степеней и 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}

](https://dic.academic.ru/dic.nsf/ruwiki/a82acaf7d5365d380a6d26ae18fbbbba.png)

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

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

//_________________________________________________________________________________________ //_________________________________________________________________________________________ // // NAME: FFT. // PURPOSE: Быстрое преобразование Фурье: Комплексный сигнал в комплексный спектр и обратно. // В случае действительного сигнала в мнимую часть (Idat) записываются нули. // Количество отсчетов - кратно 2**К - т.е. 2, 4, 8, 16, ... (см. комментарии ниже). // (C) Sergey Aleynik. saleynik@yandex.ru //
// 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;

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 { ///

/// Вычисление поворачивающего модуля e^(-i2PI*k/N) ///

/// /// /// 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)); } ///

/// Возвращает спектр сигнала ///

/// Массив значений сигнала. Количество значений должно быть степенью 2 /// Массив со значениями спектра сигнала 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; } ///

/// Центровка массива значений полученных в fft (спектральная составляющая при нулевой частоте будет в центре массива) ///

/// Массив значений полученный в fft /// 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-мерного сигнала, представленного как вещественными, так и комплексными числами. (англ.)