29-09-2023
Быстрое преобразование Фурье (БПФ, FFT) — это алгоритм быстрого вычисления дискретного преобразования Фурье (ДПФ). То есть, алгоритм вычисления за количество действий, меньшее чем , требуемых для прямого (по формуле) вычисления ДПФ. Иногда под БПФ понимается один из быстрых алгоритмов, называемый алгоритмом прореживания по частоте/времени или алгоритмом по основанию 2, имеющего сложность .
Содержание |
Покажем как выполнить дискретное преобразование Фурье за действий при . В частности, при понадобится действий.
Дискретное преобразование Фурье преобразует набор чисел в набор чисел , такой, что , где и при . Алгоритм быстрого преобразования Фурье применим к любым коммутативным ассоциативным кольцам с единицей. Чаще всего этот алгоритм применяют к полю комплексных чисел (c ) и к кольцам вычетов.
Основной шаг алгоритма состоит в сведении задачи для чисел к задаче для числам, где — делитель . Пусть мы уже умеем решать задачу для чисел. Применим преобразование Фурье к наборам для . Покажем теперь, как за действий решить исходную задачу. Заметим, что . Выражения в скобках нам уже известны — это -тое число после преобразования Фурье -той группы. Таким образом, для вычисления каждого нужно действий, а для вычисления всех — действий, что и требовалось получить.
Для обратного преобразования Фурье можно применять алгоритм прямого преобразования Фурье — нужно лишь использовать вместо (или применить операцию комплексного сопряжения в начале к входным данным, а затем к результату, полученному после прямого преобразования Фурье) и окончательный результат поделить на .
Общий случай может быть сведён к предыдущему. Пусть . Заметим, что . Обозначим . Тогда , если положить при .
Таким образом задача сведена к вычислению свёртки, но это можно сделать с помощью трёх преобразований Фурье для элементов. Выполняем прямое преобразование Фурье для и , перемножаем поэлементно результаты и выполняем обратное преобразование Фурье.
Вычисления всех и требуют действий, три преобразования Фурье требуют действий, перемножение результатов преобразований Фурье требует действий, вычисление всех зная значения свертки требует действий. Итого для дискретного преобразования Фурье требуется действий для любого .
Этот алгоритм быстрого преобразования Фурье может работать над кольцом только когда известны первообразные корни из единицы степеней и .
Дискретное преобразование Фурье для вектора , состоящего из N элементов, имеет вид:
элементы матрицы имеют вид: .
Пусть N четно, тогда ДПФ можно переписать следующим образом:
Коэффициенты и можно переписать следующим образом (M=N/2):
В результате получаем:
То есть дискретное преобразование Фурье от вектора, состоящего из N отсчетов, свелось к линейной композиции двух ДПФ от отсчетов, и если для первоначальной задачи требовалось операций, то для полученной композиции — . Если M является степенью двух, то это разделение можно продолжать рекурсивно до тех пор, пока не дойдем до двухточечного преобразования Фурье, которое вычисляется по следующим формулам:
Ниже приведен пример расчета комплексного БПФ написанный на С:
//_________________________________________________________________________________________ //_________________________________________________________________________________________ // // 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; 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_INVERSE) 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 / 2 или меньше. 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+1]; 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 >> 2; } 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 = 1; i < Nft; i++) { j = i * 2; FTvl[i] = Sqrt(Sqr(Tmvl[j]) + Sqr(Tmvl[j+1])); } delete []Tmvl; }
Пример реализации на 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[i]:= Sqrt(Sqr(Tmvl[j]) + Sqr(Tmvl[j+1])); end; SetLength(Tmvl, 0); end;
Быстрое преобразование Фурье.