Элементы спектрального анализа (Фурье, Хартман etc.)

Материал из DRKB


Рисование спектра звукового файла[править | править код]

Image 2021-10-02 181057.png

{$A+,B-,C+,D+,E-,F-,G+,H+,I+,J+,K-,L+,M-,N+,O-,P+,Q-,R-,S-,T-,U-,V+,W-,X+,Y+,Z1}
{$MINSTACKSIZE $00004000}
{$MAXSTACKSIZE $00100000}
{$IMAGEBASE $00400000}
{$APPTYPE GUI}

unit Main;

interface

uses
  Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs,
  Buttons, ExtCtrls, ComCtrls, Menus;

type
  TfmMain = class(TForm)
    MainMenu1: TMainMenu;
    N1: TMenuItem;
    N2: TMenuItem;
    StatusBar1: TStatusBar;
    N3: TMenuItem;
    imgInfo: TImage;
    Panel1: TPanel;
    btnStart: TSpeedButton;
    procedure btnStartClick(Sender: TObject);
    procedure FormCreate(Sender: TObject);
    procedure FormClose(Sender: TObject; var Action: TCloseAction);
  end;

var
  fmMain: TfmMain;

implementation

uses PFiles;

{$R *.DFM}

function Power2(lPower: Byte): LongInt;
begin
  Result := 1 shl lPower;
end;

procedure ClassicDirect(var aSignal, aSpR, aSpI: array of Double; N: LongInt);
var
  lSrch, lGarm: LongInt;
  dSumR, dSumI: Double;
begin
  for lGarm := 0 to N div 2 - 1 do
  begin
    dSumR := 0;
    dSumI := 0;
    for lSrch := 0 to N - 1 do
    begin
      dSumR := dSumR + aSignal[lSrch] * Cos(lGarm * lSrch / N * 2 * PI);
      dSumI := dSumI + aSignal[lSrch] * Sin(lGarm * lSrch / N * 2 * PI);
    end;
    aSpR[lGarm] := dSumR;
    aSpI[lGarm] := dSumI;
  end;
end;

procedure ClassicInverce(var aSpR, aSpI, aSignal: array of Double; N: LongInt);
var
  lSrch, lGarm: LongInt;
  dSum: Double;
begin
  for lSrch := 0 to N - 1 do
  begin
    dSum := 0;
    for lGarm := 0 to N div 2 - 1 do
    begin
      dSum := dSum
        + aSpR[lGarm] * Cos(lSrch * lGarm * 2 * Pi / N)
        + aSpI[lGarm] * Sin(lSrch * lGarm * 2 * Pi / N);
    end;
    aSignal[lSrch] := dSum * 2;
  end;
end;

function InvertBits(BF, DataSize, Power: Word): Word;
var
  BR, NN, L: Word;
begin
  br := 0;
  nn := DataSize;
  for l := 1 to Power do
  begin
    NN := NN div 2;
    if (BF >= NN) then
    begin
      BR := BR + Power2(l - 1);
      BF := BF - NN
    end;
  end;
  InvertBits := BR;
end;

procedure FourierDirect(var RealData, VirtData, ResultR, ResultV: array of Double; DataSize: LongInt);
var
  A1, A2, B1, B2: Real;
  D2, C2, C1, D1: Word;
  I, J, K: Word;
  Cosin, Sinus: Real;
  wIndex, Power: Word;
begin
  C1 := DataSize shr 1;
  C2 := 1;
  // hope it will be faster then round(ln(DataSize) / ln(2))
  Power := 0;
  while Power <= 15 do
  begin
    if Power2(Power) = DataSize then Break;
    Inc(Power);
  end;
  for I := 1 to Power do
  begin
    D1 := 0;
    D2 := C1;
    for J := 1 to C2 do
    begin
      wIndex := InvertBits(D1 div C1, DataSize, Power);
      Cosin := +(Cos((2 * Pi / DataSize) * wIndex));
      Sinus := -(Sin((2 * Pi / DataSize) * wIndex));
      for K := D1 to D2 - 1 do
      begin
        A1 := RealData[K];
        A2 := VirtData[K];
        B1 := ((Cosin * RealData[K + C1] - Sinus * VirtData[K + C1]));
        B2 := ((Sinus * RealData[K + C1] + Cosin * VirtData[K + C1]));
        RealData[K] := A1 + B1;
        VirtData[K] := A2 + B2;
        RealData[K + C1] := A1 - B1;
        VirtData[K + C1] := A2 - B2;
      end;
      Inc(D1, C1 * 2);
      Inc(D2, C1 * 2);
    end;
    C1 := C1 div 2;
    C2 := C2 * 2;
  end;
  for I := 0 to DataSize div 2 - 1 do
  begin
    ResultR[I] := +RealData[InvertBits(I, DataSize, Power)];
    ResultV[I] := -VirtData[InvertBits(I, DataSize, Power)];
  end;
end;

procedure Hartley(iSize: LongInt; var aData: array of Double);
type
  taDouble = array[0..MaxLongInt div SizeOf(Double) - 1] of Double;
var
  prFI, prFN, prGI: ^taDouble;
  rCos, rSin: Double;
  rA, rB, rTemp: Double;
  rC1, rC2, rC3, rC4: Double;
  rS1, rS2, rS3, rS4: Double;
  rF0, rF1, rF2, rF3: Double;
  rG0, rG1, rG2, rG3: Double;
  iK1, iK2, iK3, iK4: LongInt;
  iSrch, iK, iKX: LongInt;
begin
  iK2 := 0;
  for iK1 := 1 to iSize - 1 do
  begin
    iK := iSize shr 1;
    repeat
      iK2 := iK2 xor iK;
      if (iK2 and iK) <> 0 then Break;
      iK := iK shr 1;
    until False;
    if iK1 > iK2 then
    begin
      rTemp := aData[iK1];
      aData[iK1] := aData[iK2];
      aData[iK2] := rTemp;
    end;
  end;
  iK := 0;
  while (1 shl iK) < iSize do
    Inc(iK);
  iK := iK and 1;
  if iK = 0 then
  begin
    prFI := @aData;
    prFN := @aData;
    prFN := @prFN^[iSize];
    while Word(prFI) < Word(prFN) do
    begin
      rF1 := prFI^[0] - prFI^[1];
      rF0 := prFI^[0] + prFI^[1];
      rF3 := prFI^[2] - prFI^[3];
      rF2 := prFI^[2] + prFI^[3];
      prFI^[2] := rF0 - rF2;
      prFI^[0] := rF0 + rF2;
      prFI^[3] := rF1 - rF3;
      prFI^[1] := rF1 + rF3;
      prFI := @prFI^[4];
    end;
  end
  else
  begin
    prFI := @aData;
    prFN := @aData;
    prFN := @prFN^[iSize];
    prGI := prFI;
    prGI := @prGI^[1];
    while Word(prFI) < Word(prFN) do
    begin
      rC1 := prFI^[0] - prGI^[0];
      rS1 := prFI^[0] + prGI^[0];
      rC2 := prFI^[2] - prGI^[2];
      rS2 := prFI^[2] + prGI^[2];
      rC3 := prFI^[4] - prGI^[4];
      rS3 := prFI^[4] + prGI^[4];
      rC4 := prFI^[6] - prGI^[6];
      rS4 := prFI^[6] + prGI^[6];
      rF1 := rS1 - rS2;
      rF0 := rS1 + rS2;
      rG1 := rC1 - rC2;
      rG0 := rC1 + rC2;
      rF3 := rS3 - rS4;
      rF2 := rS3 + rS4;
      rG3 := Sqrt(2) * rC4;
      rG2 := Sqrt(2) * rC3;
      prFI^[4] := rF0 - rF2;
      prFI^[0] := rF0 + rF2;
      prFI^[6] := rF1 - rF3;
      prFI^[2] := rF1 + rF3;
      prGI^[4] := rG0 - rG2;
      prGI^[0] := rG0 + rG2;
      prGI^[6] := rG1 - rG3;
      prGI^[2] := rG1 + rG3;
      prFI := @prFI^[8];
      prGI := @prGI^[8];
    end;
  end;
  if iSize < 16 then Exit;
  repeat
    Inc(iK, 2);
    iK1 := 1 shl iK;
    iK2 := iK1 shl 1;
    iK4 := iK2 shl 1;
    iK3 := iK2 + iK1;
    iKX := iK1 shr 1;
    prFI := @aData;
    prGI := prFI;
    prGI := @prGI^[iKX];
    prFN := @aData;
    prFN := @prFN^[iSize];
    repeat
      rF1 := prFI^[000] - prFI^[iK1];
      rF0 := prFI^[000] + prFI^[iK1];
      rF3 := prFI^[iK2] - prFI^[iK3];
      rF2 := prFI^[iK2] + prFI^[iK3];
      prFI^[iK2] := rF0 - rF2;
      prFI^[000] := rF0 + rF2;
      prFI^[iK3] := rF1 - rF3;
      prFI^[iK1] := rF1 + rF3;
      rG1 := prGI^[0] - prGI^[iK1];
      rG0 := prGI^[0] + prGI^[iK1];
      rG3 := Sqrt(2) * prGI^[iK3];
      rG2 := Sqrt(2) * prGI^[iK2];
      prGI^[iK2] := rG0 - rG2;
      prGI^[000] := rG0 + rG2;
      prGI^[iK3] := rG1 - rG3;
      prGI^[iK1] := rG1 + rG3;
      prGI := @prGI^[iK4];
      prFI := @prFI^[iK4];
    until not (Word(prFI) < Word(prFN));
    rCos := Cos(Pi / 2 / Power2(iK));
    rSin := Sin(Pi / 2 / Power2(iK));
    rC1 := 1;
    rS1 := 0;
    for iSrch := 1 to iKX - 1 do
    begin
      rTemp := rC1;
      rC1 := (rTemp * rCos - rS1 * rSin);
      rS1 := (rTemp * rSin + rS1 * rCos);
      rC2 := (rC1 * rC1 - rS1 * rS1);
      rS2 := (2 * (rC1 * rS1));
      prFN := @aData;
      prFN := @prFN^[iSize];
      prFI := @aData;
      prFI := @prFI^[iSrch];
      prGI := @aData;
      prGI := @prGI^[iK1 - iSrch];
      repeat
        rB := (rS2 * prFI^[iK1] - rC2 * prGI^[iK1]);
        rA := (rC2 * prFI^[iK1] + rS2 * prGI^[iK1]);
        rF1 := prFI^[0] - rA;
        rF0 := prFI^[0] + rA;
        rG1 := prGI^[0] - rB;
        rG0 := prGI^[0] + rB;
        rB := (rS2 * prFI^[iK3] - rC2 * prGI^[iK3]);
        rA := (rC2 * prFI^[iK3] + rS2 * prGI^[iK3]);
        rF3 := prFI^[iK2] - rA;
        rF2 := prFI^[iK2] + rA;
        rG3 := prGI^[iK2] - rB;
        rG2 := prGI^[iK2] + rB;
        rB := (rS1 * rF2 - rC1 * rG3);
        rA := (rC1 * rF2 + rS1 * rG3);
        prFI^[iK2] := rF0 - rA;
        prFI^[0] := rF0 + rA;
        prGI^[iK3] := rG1 - rB;
        prGI^[iK1] := rG1 + rB;
        rB := (rC1 * rG2 - rS1 * rF3);
        rA := (rS1 * rG2 + rC1 * rF3);
        prGI^[iK2] := rG0 - rA;
        prGI^[0] := rG0 + rA;
        prFI^[iK3] := rF1 - rB;
        prFI^[iK1] := rF1 + rB;
        prGI := @prGI^[iK4];
        prFI := @prFI^[iK4];
      until not (LongInt(prFI) < LongInt(prFN));
    end;
  until not (iK4 < iSize);
end;

procedure HartleyDirect(var aData: array of Double; iSize: LongInt);
var
  rA, rB: Double;
  iI, iJ, iK: LongInt;
begin
  Hartley(iSize, aData);
  iJ := iSize - 1;
  iK := iSize div 2;
  for iI := 1 to iK - 1 do
  begin
    rA := aData[ii];
    rB := aData[ij];
    aData[iJ] := (rA - rB) / 2;
    aData[iI] := (rA + rB) / 2;
    Dec(iJ);
  end;
end;

procedure HartleyInverce(var aData: array of Double; iSize: LongInt);
var
  rA, rB: Double;
  iI, iJ, iK: LongInt;
begin
  iJ := iSize - 1;
  iK := iSize div 2;
  for iI := 1 to iK - 1 do
  begin
    rA := aData[iI];
    rB := aData[iJ];
    aData[iJ] := rA - rB;
    aData[iI] := rA + rB;
    Dec(iJ);
  end;
  Hartley(iSize, aData);
end;

// not tested
procedure HartleyDirectComplex(real, imag: array of Double; n: LongInt);
var
  a, b, c, d: Double;
  q, r, s, t: Double;
  i, j, k: LongInt;
begin
  j := n - 1;
  k := n div 2;
  for i := 1 to k - 1 do
  begin
    a := real[i]; b := real[j]; q := a + b; r := a - b;
    c := imag[i]; d := imag[j]; s := c + d; t := c - d;
    real[i] := (q + t) * 0.5; real[j] := (q - t) * 0.5;
    imag[i] := (s - r) * 0.5; imag[j] := (s + r) * 0.5;
    dec(j);
  end;
  Hartley(N, Real);
  Hartley(N, Imag);
end;

// not tested
procedure HartleyInverceComplex(real, imag: array of Double; N: LongInt);
var
  a, b, c, d: Double;
  q, r, s, t: Double;
  i, j, k: LongInt;
begin
  Hartley(N, real);
  Hartley(N, imag);
  j := n - 1;
  k := n div 2;
  for i := 1 to k - 1 do
  begin
    a := real[i]; b := real[j]; q := a + b; r := a - b;
    c := imag[i]; d := imag[j]; s := c + d; t := c - d;
    imag[i] := (s + r) * 0.5; imag[j] := (s - r) * 0.5;
    real[i] := (q - t) * 0.5; real[j] := (q + t) * 0.5;
    dec(j);
  end;
end;

procedure DrawSignal(var aSignal: array of Double; N, lColor: LongInt);
var
  lSrch: LongInt;
  lHalfHeight: LongInt;
begin
  with fmMain do
  begin
    lHalfHeight := imgInfo.Height div 2;
    imgInfo.Canvas.MoveTo(0, lHalfHeight);
    imgInfo.Canvas.Pen.Color := lColor;
    for lSrch := 0 to N - 1 do
    begin
      imgInfo.Canvas.LineTo(lSrch, Round(aSignal[lSrch]) + lHalfHeight);
    end;
    imgInfo.Repaint;
  end;
end;

procedure DrawSpector(var aSpR, aSpI: array of Double; N, lColR, lColI: LongInt);
var
  lSrch: LongInt;
  lHalfHeight: LongInt;
begin
  with fmMain do
  begin
    lHalfHeight := imgInfo.Height div 2;
    for lSrch := 0 to N div 2 do
    begin
      imgInfo.Canvas.Pixels[lSrch, Round(aSpR[lSrch] / N) + lHalfHeight] := lColR;
      imgInfo.Canvas.Pixels[lSrch + N div 2, Round(aSpI[lSrch] / N) + lHalfHeight] := lColI;
    end;
    imgInfo.Repaint;
  end;
end;

const
  N = 512;
var
  aSignalR: array[0..N - 1] of Double; //
  aSignalI: array[0..N - 1] of Double; //
  aSpR, aSpI: array[0..N div 2 - 1] of Double; //
  lFH: LongInt;

procedure TfmMain.btnStartClick(Sender: TObject);
const
  Epsilon = 0.00001;
var
  lSrch: LongInt;
  aBuff: array[0..N - 1] of ShortInt;
begin
  if lFH > 0 then
  begin
//  Repeat
    if F.Read(lFH, @aBuff, N) <> N then
    begin
      Exit;
    end;
    for lSrch := 0 to N - 1 do
    begin
      aSignalR[lSrch] := ShortInt(aBuff[lSrch] + $80);
      aSignalI[lSrch] := 0;
    end;

    imgInfo.Canvas.Rectangle(0, 0, imgInfo.Width, imgInfo.Height);
    DrawSignal(aSignalR, N, $D0D0D0);

//  ClassicDirect(aSignalR, aSpR, aSpI, N);                 // result in aSpR & aSpI, aSignal unchanged
//  FourierDirect(aSignalR, aSignalI, aSpR, aSpI, N);       // result in aSpR & aSpI, aSiggnalR & aSignalI modified

    HartleyDirect(aSignalR, N); // result in source aSignal ;-)

    DrawSpector(aSignalR, aSignalR[N div 2 - 1], N, $80, $8000);
    DrawSpector(aSpR, aSpI, N, $80, $8000);

{    for lSrch := 0 to N div 2 -1 do begin                    // comparing classic & Hartley

if (Abs(aSpR[lSrch] - aSignal[lSrch]) > Epsilon)
or ((lSrch > 0) And (Abs(aSpI[lSrch] - aSignal[N - lSrch]) > Epsilon))
then MessageDlg('Error comparing',mtError,[mbOK],-1);
end;}

    HartleyInverce(aSignalR, N); // to restore original signal with HartleyDirect
//  ClassicInverce(aSpR, aSpI, aSignalR, N);                // to restore original signal with ClassicDirect or FourierDirect

    for lSrch := 0 to N - 1 do
      aSignalR[lSrch] := aSignalR[lSrch] / N; // scaling

    DrawSignal(aSignalR, N, $D00000);
    Application.ProcessMessages;
//   Until False;

  end;
end;

procedure TfmMain.FormCreate(Sender: TObject);
begin
  lFH := F.Open('input.pcm', ForRead);
end;

procedure TfmMain.FormClose(Sender: TObject; var Action: TCloseAction);
begin
  F.Close(lFH);
end;

end.


Author: Denis Furman [000705
Source: Взято из Советов по Delphi от Валентина Озерова



Алгоритм FFT (Fast Fourier Transform)[править | править код]

Привожу 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: если это не работает у вас - значит я что-то выбросил вместе со своими глупыми комментариями :-) Итак, сложная версия:


cplx.pas[править | править код]

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
  Result.r := x;
  Result.i := y;
end;

function Sum(x, y: TComplex): TComplex;
begin
  Result.r := x.r + y.r;
  Result.i := x.i + y.i;
end;

function Difference(x, y: TComplex): TComplex;
begin
  Result.r := x.r - y.r;
  Result.i := x.i - y.i;
end;

function EiT(t: TReal): TComplex;
begin
  Result.r := cos(t);
  Result.i := sin(t);
end;

function Product(x, y: TComplex): TComplex;
begin
  Result.r := x.r * y.r - x.i * y.i;
  Result.i := x.r * y.i + x.i * y.r;
end;

function TimesReal(x: TComplex; y: TReal): TComplex;
begin
  Result.r := x.r * y;
  Result.i := x.i * y;
end;

function PlusReal(x: TComplex; y: TReal): TComplex;
begin
  Result.r := x.r + y;
  Result.i := x.i;
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.


cplxfft1.pas[править | править код]

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);

{ Перед вызовом Src и Dest ТРЕБУЕТСЯ распределение
  (Integer(1) shl Depth) * SizeOf(TScalar) байт памяти!}
procedure FFT(Depth: Word; Src: PScalars; Dest: PScalars);

implementation

{ рекурсивная часть, вызываемая при готовности FFT }
procedure DoFFT(Depth: Word;
  Src: PScalars;
  SrcSpacing: Word;
  Dest: PScalars);
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.


DemoForm.pas[править | править код]

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.


cplxfft2.pas[править | править код]

        • Версия для работы с реальными числами:
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);

{ Перед вызовом Src и Dest ТРЕБУЕТСЯ распределение
  (Integer(1) shl Depth) * SizeOf(TScalar) байт памяти! }
procedure FFT(Depth: Word; SrcR, SrcI: PScalars; DestR, DestI: PScalars);

implementation

{ рекурсивная часть, вызываемая при готовности FFT }
procedure DoFFT(Depth: Word;
  SrcR, SrcI: PScalars;
  SrcSpacing: Word;
  DestR, DestI: PScalars);
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;

finalization

  ReAllocMem(CosTable, 0);
  ReAllocMem(SinTable, 0);

end.


demofrm.pas[править | править код]

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.


Source: http://delphiworld.narod.ru
ID: 04234