Быстрое преобразование сигнала в спекр и обратно (методы Хартли, Фурье и классический)

Written on 15 Октября 2007. Posted in [Delphi] Алгоритмы

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

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;

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;

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;

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;

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;

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;

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;

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;

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;

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;

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.

delphi форум