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

Публикую присланный читателем алгоритм:

{$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;
beginResult := 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;beginfor lGarm := 0 to N div 2 - 1 do begindSumR := 0;dSumI := 0;for lSrch := 0 to N - 1 do begindSumR := 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;beginfor lSrch := 0 to N-1 do begindSum := 0;for lGarm := 0 to N div 2 -1do 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;Beginbr:= 0;nn:= DataSize;For l:= 1 To Power DoBeginNN:= NN Div 2;If (BF >= NN) ThenBeginBR:= BR + Power2(l - 1);BF:= BF - NNEnd;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;BeginC1:= DataSize Shr 1;C2:= 1;for Power:=0 to 15 //hope it will be faster thenround(ln(DataSize)/ln(2))
do if Power2(Power)=DataSizethen Break;For I:= 1 To Power Do BeginD1:= 0;D2:= C1;For J:= 1 To C2 Do BeginwIndex:=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 BeginA1:= 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 BeginResultR[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;BeginiK2:=0;For iK1:=1 To iSize-1 Do BeginiK:=iSize Shr 1;RepeatiK2:=iK2 Xor iK;If (iK2 And iK)<>0 Then Break;iK:=iK Shr 1;Until False;If iK1>iK2 Then BeginrTemp:=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 BeginprFI:=@aData;prFN:=@aData;prFN := @prFN[iSize];While Word(prFI)<Word(prFN) Do BeginrF1:=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 BeginprFI:=@aData;prFN:=@aData;prFN := @prFN[iSize];prGI:=prFI;prGI := @prGI[1];While Word(prFI)<Word(prFN) Do beginrC1:=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;RepeatInc(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];RepeatrF1:= 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 BeginrTemp:=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];RepeatrB:=(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;BeginHartley(iSize,aData);iJ:=iSize-1;iK:=iSize Div 2;For iI:=1 To iK-1 Do BeginrA:=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;BeginiJ:=iSize-1;iK:=iSize Div 2;For iI:=1 To iK-1 Do BeginrA:=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 begina := 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;beginHartley(N,real);Hartley(N,imag);j:=n-1;k:=n div 2;for i:=1 to k-1 do begina := 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;beginwith fmMain do beginlHalfHeight := imgInfo.Height Div 2;imgInfo.Canvas.MoveTo(0,lHalfHeight);imgInfo.Canvas.Pen.Color := lColor;for lSrch := 0 to N-1 do beginimgInfo.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;beginwith fmMain do beginlHalfHeight := imgInfo.Height Div 2;for lSrch := 0 to N Div 2 do beginimgInfo.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;beginif lFH > 0 then begin// Repeat
if F.Read(lFH,@aBuff,N) <> N then beginExit;end;for lSrch := 0 to N-1 do beginaSignalR[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 withHartleyDirect
// ClassicInverce(aSpR, aSpI, aSignalR, N); //to restore original
signal with ClassicDirect or FourierDirect

for lSrch := 0 to N -1do aSignalR[lSrch]:= aSignalR[lSrch]/N; //scaling
DrawSignal(aSignalR, N, $D00000);Application.ProcessMessages;// Until False;
end;end;
procedure TfmMain.FormCreate(Sender: TObject);
beginlFH := F.Open('input.pcm', ForRead);end;
procedure TfmMain.FormClose(Sender: TObject; var Action: TCloseAction);
beginF.Close(lFH);end;
end.

[000705]

Содержание раздела