Delphi 2
FFT аглоритм для Delphi 2
Данный алгоритм я воспроизвел где-то около года назад. Вероятно он не самый оптимальный, но для повышения скорости расчета наверняка потребуются более мощное аппаратное обеспечение.
Но я не думаю что алгоритм слишком плох, в нем заложено немало математических трюков. Имеется некоторое количество рекурсий, но они занимается не копированием данных, а манипуляциями с указателями, если у нас есть массив размером 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;
begind:= 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;
beginif 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]