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

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