Полезное для программистов:

Фриланс
Новости
Статьи
   
Рубрики:


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

Поиск:
Код

{$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

--------------------------------------------------------------------------------

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






Просмотров: 1709

 

 

Новые статьи:


Популярные:
  1. Как сделать цикличным проигрывание MIDI-файла?
  2. Создание AVI файла из рисунков
  3. Как устройство "отключить в данной конфигурации"?
  4. Kто в данный момент присоединен через Сеть?
  5. Как узнать количество доступной памяти?
  6. Как реализовать в RichEdit разноцветный текст?
  7. Как скрыть свое приложение от ProcessViewer
  8. Как программно нажать/скрыть/показ кнопку "Start"?
  9. Модуль работы с ресурсами в PE файлах
10. Функции вызова диалоговых окон выбора
11. Проверка граматики средствами Word'а из Delphi.
12. Модуль для упрощенного вызова сообщений
13. Функции для записи и чтение своих данных в, ЕХЕ- файле
14. Рекурсивный просмотр директорий
15. Network Traffic Monitor
16. Разные модули
17. Универсальная функция для обращения к любым экспортируем функциям DLL
18. Библиотека от VladS
19. Протектор для UPX'а
20. Еще об ICQ, сообщения по контакт листу?
21. Использование открытых интерфейсов
22. Теория и практика использования RTTI
23. Работа с TApplication
24. Примеры использования Drag and Drop для различных визуальных компонентов
25. Что такое порт? Правила для работы с портами
26. Симфония на клавиатуре
27. Загрузка DLL
28. Исправление автоинкремента
29. Взаимодействие с чужими окнами
30. Проверить дубляжи в столбце


 

 

 
 
На главную