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

Материал из Википедии — свободной энциклопедии

Быстрое преобразование Фурье (БПФ) — это быстрый алгоритм вычисления дискретного преобразования Фурье (см. Дискретное преобразование Фурье, ДПФ). То есть алгоритм вычисления за количество действий, меньше чем $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 (N p)$ действий решить исходную задачу. Заметим, что $b_i=\sum_{j=0}^{q-1} \varepsilon^{ij} (\sum_{k=0}^{p-1}a_{kq+j}\varepsilon^{kiq})$ . Выражения в скобках нам уже известны — это $i (mod p)$ -тое число после преобразования Фурье $j$ -той группы. Таким образом, для вычисления каждого $b i$ нужно $O (q)$ действий, а для вычисления всех $b i$ — $O (N q)$ действий, что и требовалось получить.

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

Для обратного преобразования Фурье можно применять алгоритм прямого преобразования Фурье — нужно лишь использовать $\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$ .

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

[править] Пример программы

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

#include <сmath>
#include <vector>
#include <iostream>
 
using namespace std;
 
#ifndef M_PI 
#define M_PI 3.14159265358979323846
#endif
 
vector<double> FFT(const vector<int>& dIn, int nn, int beginData)
{
	int i, j, n, m, mmax, istep;
	double tempr, tempi, wtemp, theta, wpr, wpi, wr, wi;
 
	int isign = -1;
	vector<double> data(nn*2 + 1);
 
	j = 0;
	for (i = beginData; i < beginData + nn; i++)
	{
		if (i < dIn.size())
		{
			data[j*2]   = 0;
			data[j*2+1] = dIn[i];
		}
		else
		{
			data[j*2]   = 0;
			data[j*2+1] = 0;
		}
		j++;
	}
 
	n = nn << 1;
	j = 1;
	i = 1;
	while (i < n)
	{
		if (j > i)
		{
			tempr = data[i];   data[i]   = data[j];   data[j]   = tempr;
			tempr = data[i+1]; data[i+1] = data[j+1]; data[j+1] = tempr;
		}
		m = n >> 1;
		while ((m >= 2) && (j > m))
		{
			j = j - m;
			m = m >> 1;
		}
		j = j + m;
		i = i + 2;
	}
	mmax = 2;
	while (n > mmax)
	{
		istep = 2 * mmax;
		theta = 2.0*M_PI / (isign * mmax);
		wtemp = sin(0.5 * theta);
		wpr   = -2.0 * wtemp * wtemp;
		wpi   = sin(theta);
		wr    = 1.0;
		wi    = 0.0;
		m    = 1;
		while (m < mmax)
		{
			i = m;
			while (i < n)
			{
				j         = i + mmax;
				tempr     = wr * data[j] - wi * data[j+1];
				tempi     = wr * data[j+1] + wi * data[j];
				data[j]   = data[i] - tempr;
				data[j+1] = data[i+1] - tempi;
				data[i]   = data[i] + tempr;
				data[i+1] = data[i+1] + tempi;
				i         = i + istep;
			}
			wtemp = wr;
			wr    = wtemp * wpr - wi * wpi + wr;
			wi    = wi * wpr + wtemp * wpi + wi;
			m     = m + 2;
		}
		mmax = istep;
	}
	vector<double> dOut(nn / 2);
 
	for (i = 0; i < (nn / 2); i++)
	{
		dOut[i] = sqrt( data[i*2] * data[i*2] + data[i*2+1] * data[i*2+1] );
	}
 
	return dOut;
}
 
int main()
{
	vector<int> dsin;
 
	for (double x = 0; x <= 10*M_PI; x += M_PI/180.0)
	{
		dsin.push_back( sin(x)*1000 + cos(0.5*x)*1000 );
	}
 
	vector<double> dfourier = FFT(dsin, 1024, 1);
 
	for (int i = 0; i < dfourier.size(); i++)
	{
		cout << i << "\t" << dfourier[i] << endl;
	}
	return 0;
}

Реализация алгоритма быстрого преобразования Фурье на языке Pascal:

procedure FFT(var a : TSingleArray;
     nn : Integer;
     InverseFFT : Boolean);
var
    ii, jj, n, mmax, m, j, istep, i, isign : Integer;
    wtemp, wr, wpr, wpi, wi, theta, tempr, tempi : Double;
begin
    if InverseFFT then isign := -1
    else isign := 1;
    n := 2*nn; j := 1; ii:=1;
    while ii <= nn do
    begin
        i := 2*ii-1;
        if j>i then
        begin
            tempr := a[j-1];
            tempi := a[j];
            a[j-1] := a[i-1];
            a[j] := a[i];
            a[i-1] := tempr;
            a[i] := tempi;
        end;
        m := n div 2;
        while (m>=2) and (j>m) do
        begin
            j := j-m;
            m := m div 2;
        end;
        j := j+m;
        Inc(ii);
    end;
    mmax := 2;
    while n>mmax do
    begin
        istep := 2*mmax;
        theta := 2*Pi/(isign*mmax);
        wpr := -2.0*sqr(sin(0.5*theta));
        wpi := sin(theta);
        wr := 1.0;
        wi := 0.0;
        ii:=1;
        while ii<=mmax div 2 do
        begin
            m := 2*ii-1;
            jj:=0;
            while jj<=(n-m) div istep do
            begin
                i := m+jj*istep;
                j := i+mmax;
                tempr := wr*a[j-1]-wi*a[j];
                tempi := wr*a[j]+wi*a[j-1];
                a[j-1] := a[i-1]-tempr;
                a[j] := a[i]-tempi;
                a[i-1] := a[i-1]+tempr;
                a[i] := a[i]+tempi;
                Inc(jj);
            end;
            wtemp := wr;
            wr := wr*wpr-wi*wpi+wr;
            wi := wi*wpr+wtemp*wpi+wi;
            Inc(ii);
        end;
        mmax := istep;
    end;
    if InverseFFT then
    begin
        I:=1;
        while I<=2*nn do
        begin
            a[I-1] := a[I-1]/nn;
            Inc(I);
        end;
    end;
end;

[править] Ссылки

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

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

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

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

Материал из Википедии — свободной энциклопедии

Содержание

[править] Основной алгоритм

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

[править] Общий случай

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

[править] Пример программы

[править] Ссылки

Просмотры

Личные инструменты

Навигация

Поиск

Участие

Инструменты

На других языках