Привожу 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,
происходит вызов
, далее заполняем 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. | |
