Delphi World - FFT аглоритм для Delphi2

Delphi World - это проект, являющийся сборником статей и малодокументированных возможностей по программированию в среде Delphi. Здесь вы найдёте работы по следующим категориям: delphi, delfi, borland, bds, дельфи, делфи, дэльфи, дэлфи, programming, example, программирование, исходные коды, code, исходники, source, sources, сорцы, сорсы, soft, programs, программы, and, how, delphiworld, базы данных, графика, игры, интернет, сети, компоненты, классы, мультимедиа, ос, железо, программа, интерфейс, рабочий стол, синтаксис, технологии, файловая система...

FFT аглоритм для Delphi2

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