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

Previous  Top  Next

    
 

 

 

Code:

{$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: LongInt;

var lGarm: LongInt;

var dSumR: Double;

var 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: LongInt;

var lGarm: LongInt;

var dSum: Double;

begin

for lSrch := 0 to N - 1 do

   begin

     dSum := 0;

     for lGarm := 0 to N div 2 - 1 do

       dSum := dSum

         + aSpR[lGarm] * Cos(lSrch * lGarm * 2 * Pi / N)

         + aSpI[lGarm] * Sin(lSrch * lGarm * 2 * Pi / N);

     aSignal[lSrch] := dSum * 2;

   end;

end;

 

function InvertBits(BF, DataSize, Power: Word): Word;

 

var BR: Word;

var NN: Word;

var 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: Real;

var A2: Real;

var B1: Real;

var B2: Real;

var D2: Word;

var C2: Word;

var C1: Word;

var D1: Word;

var I: Word;

var J: Word;

var K: Word;

var Cosin: Real;

var Sinus: Real;

var wIndex: Word;

var Power: Word;

begin

C1 := DataSize shr 1;

C2 := 1;

for Power := 0 to 15 //hope it will be faster then

   round(ln(DataSize) / ln(2)) do

   if Power2(Power) = DataSize then Break;

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;

var rCos, rSin: Double;

var rA, rB, rTemp: Double;

var rC1, rC2, rC3, rC4: Double;

var rS1, rS2, rS3, rS4: Double;

var rF0, rF1, rF2, rF3: Double;

var rG0, rG1, rG2, rG3: Double;

var iK1, iK2, iK3, iK4: LongInt;

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

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

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

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

var 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; //

var aSignalI: array[0..N - 1] of Double; //

var aSpR, aSpI: array[0..N div 2 - 1] of Double; //

var lFH: LongInt;

 

procedure TfmMain.btnStartClick(Sender: TObject);

 

const Epsilon = 0.00001;

var lSrch: LongInt;

var 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.

 

 

Denis Furman [000705

Взято из Советов по Delphi от Валентина Озерова

Сборник Kuliba

 

 


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

 

Code:

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.

 

 

 

 

Code:

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.

 

 

 

 

 

 

Code:

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.

 

 

 

 

 

 

 

 

 

**** Версия для работы с реальными числами:

 

Code:

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.

 

Code:

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.

 

 

 

 

 

©Drkb::04234

 

 

 

       

Взято с http://delphiworld.narod.ru