позволяющий оперировать 256 точками данных

Привожу 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])

<
/p> , где 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 = recordr : 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 dobeginr:=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.

<

/p>

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
beginDest^[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 thenraise exception.Create('глубина рекурсии должны быть положительным целым числом');

N:= integer(1) shl d ;
GetMem(Src, N*Sizeof(TScalar));GetMem(Dest, N*SizeOf(TScalar));
for j:=0 to N-1 dobeginsrc^[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.

<

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

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
beginDestR^[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
beginDestR^[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
beginCosTable^[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.

<

/p>

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 thenraise 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]

Содержание раздела