Delphi World - Преобразование сигнала в спекр и обратно (методы Хартли, Фурье и классический)
Delphi World - это проект, являющийся сборником статей и малодокументированных возможностей  по программированию в среде Delphi. Здесь вы найдёте работы по следующим категориям: delphi, delfi, borland, bds, дельфи, делфи, дэльфи, дэлфи, programming, example, программирование, исходные коды, code, исходники, source, sources, сорцы, сорсы, soft, programs, программы, and, how, delphiworld, базы данных, графика, игры, интернет, сети, компоненты, классы, мультимедиа, ос, железо, программа, интерфейс, рабочий стол, синтаксис, технологии, файловая система...
Преобразование сигнала в спекр и обратно (методы Хартли, Фурье и классический)

Автор: Denis Furman


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

Проект Delphi World © Выпуск 2002 - 2017
Автор проекта: Эксклюзивные курсы программирования