Привожу FFT-алгоритм, позволяющий оперировать 256 точками данных примерно за 0.008 секунд на P66 (с 72MB, YMMV). Создан на Delphi.

Данный алгоритм я воспроизвел где-то около года назад. Вероятно он не самый оптимальный, но для повышения скорости расчета наверняка потребуются более мощное аппаратное обеспечение.

Но я не думаю что алгоритм слишком плох, в нем заложено немало математических трюков. Имеется некоторое количество рекурсий, но они занимается не копированием данных, а манипуляциями с указателями, если у нас есть массив размером N = 2^d, то глубина рекурсии составит всего d. Возможно имело бы смысл применить развертывающуюся рекурсию, но не пока не ясно, поможет ли ее применение в данном алгоритме. (Но вероятно мы смогли бы достаточно легко получить надежную математическую модель, развертывая в рекурсии один или два нижних слоя, то есть проще говоря:

if Depth < 2 then
{производим какие-либо действия}

вместо текущего 'if Depth = 0 then...' Это должно устранить непродуктивные вызовы функций, что несомненно хорошо в то время, пока развертывающая рекурсия работает с ресурсами.)

Имеется поиск с применением таблиц синусов и косинусов; здесь использован метод золотой середины: данный алгоритм весьма трудоемок, но дает отличные результаты при использовании малых и средних массивов.

Вероятно в машине с большим объемом оперативной памяти следует использовать VirtualAlloc(... PAGE_NOCACHE) для Src, Dest и таблиц поиска.

Если кто-либо обнаружит неверную на ваш взгляд или просто непонятную в данном совете функцию пожалуйста сообщите мне об этом.

Что делает данная технология вкратце. Имеется несколько FFT, образующих 'комплексный FT', который понимает и о котором заботится моя технология. Это означает, что если N = 2^d, Src^ и Dest^ образуют массив из N TComplexes, происходит вызов

FFT(d, Src, Dest)

, далее заполняем Dest с применением 'комплексного FT' после того, как результат вызова Dest^[j] будет равен

1/sqrt(N) * Sum(k=0.. N - 1 ; EiT(2*Pi(j*k/N)) * Src^[k])

, где EiT(t) = cos(t) + i sin(t) . То есть, стандартное преобразование Фурье.

Публикую две версии: в первой версии я использую TComplex с функциями для работы с комплексными числами. Во второй версии все числа реальные - вместо массивов Src и Dest мы используем массивы реальных чисел SrcR, SrcI, DestR, DestI (в блоке вычислений реальных чисел), и вызовы всех функций осуществляются линейно. Первая версия достаточна легка в реализации, зато вторая - значительно быстрее. (Обе версии оперируют 'комплексными FFT'.) Технология работы была опробована на алгоритме Plancherel (также известным как Parseval). Обе версии работоспособны, btw: если это не работает у вас - значит я что-то выбросил вместе со своими глупыми коментариями :-) Итак, сложная версия:

unit cplx;

interface

type
PReal = ^TReal; TReal = extended;
PComplex = ^TComplex; TComplex = record r : TReal; i : TReal; end;

function MakeComplex(x, y: TReal): TComplex;
function Sum(x, y: TComplex) : TComplex;
function Difference(x, y: TComplex) : TComplex;
function Product(x, y: TComplex): TComplex;
function TimesReal(x: TComplex; y: TReal): TComplex;
function PlusReal(x: TComplex; y: TReal): TComplex;
function EiT(t: TReal):TComplex;
function ComplexToStr(x: TComplex): string;
function AbsSquared(x: TComplex): TReal;

implementation

uses SysUtils;

function MakeComplex(x, y: TReal): TComplex;
begin
with result do begin r:=x; i:= y; end; end;

function Sum(x, y: TComplex) : TComplex;
begin
with result do
begin
r:= x.r + y.r; i:= x.i + y.i; end;
end;

function Difference(x, y: TComplex) : TComplex;
begin
with result do
begin
r:= x.r - y.r; i:= x.i - y.i; end;
end;

function EiT(t: TReal): TComplex;
begin
with result do
begin
r:= cos(t); i:= sin(t); end;
end;

function Product(x, y: TComplex): TComplex;
begin
with result do
begin
r:= x.r * y.r - x.i * y.i; i:= x.r * y.i + x.i * y.r; end;
end;

function TimesReal(x: TComplex; y: TReal): TComplex;
begin
with result do
begin
r:= x.r * y; i:= x.i * y; end;
end;

function PlusReal(x: TComplex; y: TReal): TComplex;
begin
with result do
begin
r:= x.r + y; i:= x.i; end;
end;

function ComplexToStr(x: TComplex): string;
begin
result:= FloatToStr(x.r) + ' + ' + FloatToStr(x.i) + 'i'; end;

function AbsSquared(x: TComplex): TReal;
begin
result:= x.r*x.r + x.i*x.i; end;

end.

unit cplxfft1;

interface

uses Cplx;

type
PScalar = ^TScalar; TScalar = TComplex; {Легко получаем преобразование в реальную величину}
PScalars = ^TScalars; TScalars = array[0..High(integer) div SizeOf(TScalar) - 1] of TScalar;
const

TrigTableDepth: word = 0; TrigTable : PScalars = nil;
procedure InitTrigTable(Depth: word);

procedure FFT(Depth: word;
Src: PScalars; Dest: PScalars);
{Перед вызовом Src и Dest ТРЕБУЕТСЯ распределение

(integer(1) shl Depth) * SizeOf(TScalar)

байт памяти!}

implementation

procedure DoFFT(Depth: word;
Src: PScalars; SrcSpacing: word; Dest: PScalars); {рекурсивная часть, вызываемая при готовности FFT}
var j, N: integer; Temp: TScalar; Shift: word;
begin
if Depth = 0 then
begin Dest^[0]:= Src^[0]; exit; end;
N:= integer(1) shl (Depth - 1);

DoFFT(Depth - 1, Src, SrcSpacing * 2, Dest);
DoFFT(Depth - 1, @Src^[SrcSpacing], SrcSpacing * 2, @Dest^[N] );

Shift:= TrigTableDepth - Depth;

for j:= 0 to N - 1 do
begin
Temp:= Product(TrigTable^[j shl Shift], Dest^[j + N]); Dest^[j + N]:= Difference(Dest^[j], Temp); Dest^[j] := Sum(Dest^[j], Temp); end;

end;

procedure FFT(Depth: word;
Src: PScalars; Dest: PScalars); var j, N: integer; Normalizer: extended;
begin

N:= integer(1) shl depth;

if Depth TrigTableDepth then
InitTrigTable(Depth);
DoFFT(Depth, Src, 1, Dest);

Normalizer:= 1 / sqrt(N) ;

for j:=0 to N - 1 do
Dest^[j]:= TimesReal(Dest^[j], Normalizer);
end;

procedure InitTrigTable(Depth: word);
var j, N: integer;
begin

N:= integer(1) shl depth;
ReAllocMem(TrigTable, N * SizeOf(TScalar));
for j:=0 to N - 1 do
TrigTable^[j]:= EiT(-(2*Pi)*j/N); TrigTableDepth:= Depth;

end;

initialization

;
finalization
ReAllocMem(TrigTable, 0);
end.

unit DemoForm;

interface

uses
Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs, StdCtrls;
type
TForm1 = class(TForm) Button1: TButton; Memo1: TMemo; Edit1: TEdit; Label1: TLabel; procedure Button1Click(Sender: TObject); private { Private declarations } public { Public declarations } end;
var
Form1: TForm1;
implementation

{$R *.DFM}

uses cplx, cplxfft1, MMSystem;

procedure TForm1.Button1Click(Sender: TObject);
var j: integer; s:string;
src, dest: PScalars; norm: extended; d,N,count:integer; st,et: longint; begin

d:= StrToIntDef(edit1.text, -1) ; if d <1 then raise exception.Create('глубина рекурсии должны быть положительным целым числом');

N:= integer(1) shl d ;
GetMem(Src, N*Sizeof(TScalar)); GetMem(Dest, N*SizeOf(TScalar));
for j:=0 to N-1 do begin src^[j]:= MakeComplex(random, random); end;
begin

st:= timeGetTime; FFT(d, Src, dest); et:= timeGetTime;
end;

Memo1.Lines.Add('N = ' + IntToStr(N)); Memo1.Lines.Add('норма ожидания: ' +#9+ FloatToStr(N*2/3));
norm:=0; for j:=0 to N-1 do norm:= norm + AbsSquared(src^[j]); Memo1.Lines.Add('Норма данных: '+#9+FloatToStr(norm)); norm:=0; for j:=0 to N-1 do norm:= norm + AbsSquared(dest^[j]); Memo1.Lines.Add('Норма FT: '+#9#9+FloatToStr(norm));

Memo1.Lines.Add('Время расчета FFT: '+#9 + inttostr(et - st) + ' мс.'); Memo1.Lines.Add(' ');
FreeMem(Src); FreeMem(DEst); end;

end.

**** Версия для работы с реальными числами:

unit cplxfft2;

interface

type
PScalar = ^TScalar; TScalar = extended;
PScalars = ^TScalars; TScalars = array[0..High(integer) div SizeOf(TScalar) - 1] of TScalar;
const

TrigTableDepth: word = 0; CosTable : PScalars = nil; SinTable : PScalars = nil;
procedure InitTrigTables(Depth: word);

procedure FFT(Depth: word;
SrcR, SrcI: PScalars; DestR, DestI: PScalars);
{Перед вызовом Src и Dest ТРЕБУЕТСЯ распределение

(integer(1) shl Depth) * SizeOf(TScalar)

байт памяти!}

implementation

procedure DoFFT(Depth: word;
SrcR, SrcI: PScalars; SrcSpacing: word; DestR, DestI: PScalars); {рекурсивная часть, вызываемая при готовности FFT}
var j, N: integer;
TempR, TempI: TScalar; Shift: word; c, s: extended; begin
if Depth = 0 then
begin DestR^[0]:= SrcR^[0]; DestI^[0]:= SrcI^[0]; exit; end;
N:= integer(1) shl (Depth - 1);

DoFFT(Depth - 1, SrcR, SrcI, SrcSpacing * 2, DestR, DestI);
DoFFT(Depth - 1,
@SrcR^[srcSpacing], @SrcI^[SrcSpacing], SrcSpacing * 2, @DestR^[N], @DestI^[N]);
Shift:= TrigTableDepth - Depth;

for j:= 0 to N - 1 do
begin
c:= CosTable^[j shl Shift]; s:= SinTable^[j shl Shift];
TempR:= c * DestR^[j + N] - s * DestI^[j + N]; TempI:= c * DestI^[j + N] + s * DestR^[j + N];
DestR^[j + N]:= DestR^[j] - TempR; DestI^[j + N]:= DestI^[j] - TempI;
DestR^[j]:= DestR^[j] + TempR; DestI^[j]:= DestI^[j] + TempI; end;

end;

procedure FFT(Depth: word;
SrcR, SrcI: PScalars; DestR, DestI: PScalars); var j, N: integer; Normalizer: extended;
begin

N:= integer(1) shl depth;

if Depth TrigTableDepth then
InitTrigTables(Depth);
DoFFT(Depth, SrcR, SrcI, 1, DestR, DestI);

Normalizer:= 1 / sqrt(N) ;

for j:=0 to N - 1 do
begin DestR^[j]:= DestR^[j] * Normalizer; DestI^[j]:= DestI^[j] * Normalizer; end;
end;

procedure InitTrigTables(Depth: word);
var j, N: integer;
begin

N:= integer(1) shl depth;
ReAllocMem(CosTable, N * SizeOf(TScalar));
ReAllocMem(SinTable, N * SizeOf(TScalar));
for j:=0 to N - 1 do
begin CosTable^[j]:= cos(-(2*Pi)*j/N); SinTable^[j]:= sin(-(2*Pi)*j/N); end; TrigTableDepth:= Depth;

end;

initialization

;
finalization
ReAllocMem(CosTable, 0); ReAllocMem(SinTable, 0);
end.

unit demofrm;

interface

uses
Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs, cplxfft2, StdCtrls;
type
TForm1 = class(TForm) Button1: TButton; Memo1: TMemo; Edit1: TEdit; Label1: TLabel; procedure Button1Click(Sender: TObject); private { Private declarations } public { Public declarations } end;
var
Form1: TForm1;
implementation

{$R *.DFM}

uses MMSystem;

procedure TForm1.Button1Click(Sender: TObject);
var SR, SI, DR, DI: PScalars;
j,d,N:integer;
st, et: longint;
norm: extended;
begin

d:= StrToIntDef(edit1.text, -1) ; if d <1 then raise exception.Create('глубина рекурсии должны быть положительным целым числом');
N:= integer(1) shl d;

GetMem(SR, N * SizeOf(TScalar));
GetMem(SI, N * SizeOf(TScalar));
GetMem(DR, N * SizeOf(TScalar));
GetMem(DI, N * SizeOf(TScalar));

for j:=0 to N - 1 do
begin
SR^[j]:=random; SI^[j]:=random; end;

st:= timeGetTime; FFT(d, SR, SI, DR, DI);
et:= timeGetTime;
memo1.Lines.Add('N = '+inttostr(N));
memo1.Lines.Add('норма ожидания: '+#9+FloatToStr(N*2/3));

norm:=0;
for j:=0 to N - 1 do
norm:= norm + SR^[j]*SR^[j] + SI^[j]*SI^[j]; memo1.Lines.Add('норма данных: '+#9+FloatToStr(norm));

norm:=0;
for j:=0 to N - 1 do
norm:= norm + DR^[j]*DR^[j] + DI^[j]*DI^[j]; memo1.Lines.Add('норма FT: '+#9#9+FloatToStr(norm));

memo1.Lines.Add('Время расчета FFT: '+#9+inttostr(et-st));
memo1.Lines.add('');
(*for j:=0 to N - 1 do
Memo1.Lines.Add(FloatToStr(SR^[j]) + ' + ' + FloatToStr(SI^[j]) + 'i');
for j:=0 to N - 1 do
Memo1.Lines.Add(FloatToStr(DR^[j]) + ' + ' + FloatToStr(DI^[j]) + 'i');*)
FreeMem(SR, N * SizeOf(TScalar));
FreeMem(SI, N * SizeOf(TScalar));
FreeMem(DR, N * SizeOf(TScalar));
FreeMem(DI, N * SizeOf(TScalar));
end;

end.

[000111]