Главная страницаОбратная связьКарта сайта

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.


Обсудить статью на форуме


Если Вас заинтересовала или понравилась информация по разработке на Delph - "FFT аглоритм для Delphi2", Вы можете поставить закладку в социальной сети или в своём блоге на данную страницу:

Так же Вы можете задать вопрос по работе этого модуля или примера через форму обратной связи, в сообщение обязательно указывайте название или ссылку на статью!
   


Copyright © 2008 - 2024 Дискета.info