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

Delphi , Синтаксис , Математика

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.

Перевод:

Это реализация алгоритма быстрого преобразования Фурье (FFT) в Delphi, написанная пользователем, который утверждает, что она была воспроизведена из неизвестного источника примерно год назад.

Код состоит из двух модулей: cplx и cplxfft1, которые отвечают за комплексные числа и FFT соответственно. Второй модуль включает рекурсивные функциональные вызовы для расчета FFT, а также таблицы тригонометрических функций, инициализируются перед вызовом FFT.

Также есть демонстрационный форм (demofrm) с кнопкой, которая выполняет FFT на массиве случайных комплексных чисел и отображает результат в поле заметок.

Код использует рекурсию для расчета FFT. Процедура DoFFT рекурсивно рассчитывает FFT, делив проблему на более мелкие подпроблемы, пока не достигнет базового случая (когда глубина равна 0). Тригонометрические таблицы используются для расчета комплексных экспоненциальных факторов, необходимых при расчете FFT.

Однако, есть несколько проблем с этим кодом:

Он не очень эффективен для больших данных из-за рекурсивных функциональных вызовов.
Нет проверок ошибок или обработки ошибок в случае недопустимого ввода или ошибок аллокации памяти.
Демонастрационный форм только отображает результат как текст, более полезно было бы отобразить FFT графически.

Вот некоторые предложения по улучшению кода:

Используйте цикл вместо рекурсии для лучшей производительности и предотвращения потенциальных переполнений стека.
Добавьте проверки ошибок и обработку ошибок в случае недопустимого ввода или ошибок аллокации памяти.
Рассмотрите использование библиотеки,such as FastFFTW, которая предоставляет оптимизированную реализацию FFT в Delphi.

Вторая версия (cplxfft2) для вещественных чисел очень похожа на первую, но с модификациями для работы с вещественными числами вместо комплексных. Она также использует тригонометрические таблицы для расчета FFT. Однако, она имеет те же проблемы, что и первая версия.

FFT алгоритм для Delphi2: Привожу FFT-алгоритм, позволяющий оперировать 256 точками данных примерно за 0.008 секунд на P66 (с 72MB, YMMV). Создан на Delphi.

Комментарии и вопросы

Получайте свежие новости и обновления по Object Pascal, Delphi и Lazarus прямо в свой смартфон. Подпишитесь на наш Telegram-канал delphi_kansoftware и будьте в курсе последних тенденций в разработке под Linux, Windows, Android и iOS

Материалы статей собраны из открытых источников, владелец сайта не претендует на авторство. Там где авторство установить не удалось, материал подаётся без имени автора. В случае если Вы считаете, что Ваши права нарушены, пожалуйста, свяжитесь с владельцем сайта.

:: Главная :: Математика ::

Время компиляции файла: 2024-08-19 13:29:56
2024-11-21 13:17:12/0.0069599151611328/1