Ñîâåòû ïî 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; 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]
