» Главная
eXcode.ru » Статьи » Delphi/Pascal » Большой FAQ по Delphi » Синтаксис » Преобразования
» Новости
» Опросы
» Файлы
» Журнал



Пользователей: 0
Гостей: 7





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




Автор: 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.

К началу статьи





Добавил: LedWormДата публикации: 2005-11-10 19:44:40

Рейтинг статьи:0.00 [Голосов 0]Кол-во просмотров: 8390
Оцените статью:

Комментарии читателей

Всего комментариев: 1029

2017-10-18 09:21:10
Doriankiz
Всё зависит только от креативности мастера – а технология дальнейшей работы при этом существенно не меняется <a href="http://genmontage.ru/articles/teplaja-shtukaturka-dlja-vnutrennih-rabot-vidy.html">genmontage.ru</a>.

2017-10-18 04:06:37
RickyEmova
Ведущий Российский интернет магазин с огромным выбором высококачественных товаров по низким ценам http://bit.ly/2kSXFpB
[url=http://bit.ly/2kSXFpB]Конструктор-присоска Palpi[/url]

2017-10-17 22:51:26
VaidokStuts

Самые доступные заменители Виагры Большие проблемы с эрекцией в среднем и старшем возрасте испытывает огромное число мужчин, но прогрессивная медицина очень эффективно борется с данной проблемой с помощью Виагры. Виагра сама по себе — великолепный препарат с очень небольшим спектром побочных эффектов, но ее весомым минусом стала высокая стоимость. Ниже мы поговорим о некоторых, наиболее эффективных заменителях обычной Виагры:

Сиалис — это самый распространенный аналог патентованного лекарства для усиления эррекции на отечественном рынке. Купить указанное лекарство а ещё [url=http://via-sexgra.ru/levitra-kupit-moskva]дженерик левитра отзывы мужчин[/url] по наиболее низкой цене вы можете на ресурсе via-sexgra.ru

Употребляется и дозируется хороший сиалис точно так же, как и классическая Виагра.

Многие знают, что после Виагры было было разработано множество препаратов с сходным влиянием на организм. Самым сложным препаратом такого рода является Левитра, что действует на эррекцию намного быстрее, при этом меньше влияя на сердечную систему мужчины. Западные доктора, по результатам последних опросов, начинают все чаще назначать Левитру вместо простой Виагры.

Шпанская мушка в наше время является лучшим природным афродизиаком со сходным с Сиалисом эффектом, который подойдет и женщинам и мужчинам. Самой привлекательной стороной этого препарата считается минимум побочных эффектов, однако стоит он чуть больше чем другие препараты. Все описанные в материале дженерики вы сумеете заказать через интернет и заказать доставку почтой. Посылка будет на на 100 процентов анонимной, данные покупателя удаляются сразу после подтверждения доставки.

2017-10-17 03:48:18
DadgilTanda
Черная маска от прыщей

[url=http://reals-gooods.ru/blackmask-afy4/?ref=26924&lnk=1168874][img]http://s02.radikal.ru/i175/1707/62/d91f4d827e9d.jpg[/img][/url]
Черная маска ATF от проблемной кожи

[url=http://reals-gooods.ru/blackmask-afy4/?ref=26924&lnk=1168874][img]http://s019.radikal.ru/i634/1703/80/511fb7c108bc.png[/img][/url]



[url=http://bit.ly/2kSXFpB]купить mp3 диски[/url]






=xxx$$=

2017-10-16 09:07:25
Ymkmvvf
Мы платим за лайки! - Оплата ежедневно!

Наш сервис предоставляет настоящие лайки на фото заказчиков, которые готовы платить за качество.

Именно для этого мы и набираем удалённых сотрудников, которые будут выполнять работу, то есть ставить лайки и зарабатывать за это деньги.

Чтобы стать нашим удалённым сотрудником и начать ставить лайки, зарабатывая при этом 45 рублей за 1 поставленный лайк,

достаточно просто зарегистрироваться на нашем сервисе. > www.oplata-vklike.tk <

Вывод заработанных средств ежедневно в течении нескольких минут.

2017-10-16 05:01:42
RickyEmova
Ведущий Российский интернет магазин с огромным выбором высококачественных товаров по низким ценам http://bit.ly/2kSXFpB
[url=http://bit.ly/2kSXFpB]Shimaki браслет[/url]

2017-10-15 15:50:24
Cscgiaa
Пройдите 1 регистрацию, сделайте 2 действия
и получайте от 8950 рублей каждые сутки в автоматическом режиме.

Мы гарантируем:
- Первый заработок в течении 60 минут.
- Стабильный доход 24 часа в сутки.
- Поступление денег без задержек.
- Для России, стран СНГ и Европы.
- Без вложений и установки ПО.
- Специально для новичков и людей без опыта.

Более подробная информация на нашем сайте. ( www.prosto-dengi.tk )

2017-10-15 15:42:17
Snfgtkq
Пройдите 1 регистрацию, сделайте 2 действия
и получайте от 8950 рублей каждые сутки в автоматическом режиме.

Мы гарантируем:
- Первый заработок в течении 60 минут.
- Стабильный доход 24 часа в сутки.
- Поступление денег без задержек.
- Для России, стран СНГ и Европы.
- Без вложений и установки ПО.
- Специально для новичков и людей без опыта.

Более подробная информация на нашем сайте. ( prosto-dengi.tk )

2017-10-15 13:46:04
Gydrgxb
Предлагаем Вам работу без вложений, на системе автоматического приёма и обработки заказов.

Мы предоставляем:

- Наш лицезионный софт.
- документы со всей необходимой доп. информацией.
- постоянная техническая поддержка.

Оплата от 5500 в день. Выплаты ежедневно.

Более подробная информация на нашем сайте >> www.obrabotka-zakazow.tk <<

2017-10-15 13:38:19
Urskqtc
Предлагаем Вам работу без вложений, на системе автоматического приёма и обработки заказов.

Мы предоставляем:

- Наш лицезионный софт.
- документы со всей необходимой доп. информацией.
- постоянная техническая поддержка.

Оплата от 5500 в день. Выплаты ежедневно.

Более подробная информация на нашем сайте >> http://obrabotka-zakazow.tk/ <<
Ваше имя: *
Текст записи: *
Имя:

Пароль:



Регистрация

Вы купите понравившуюсь программу?
Да, если не найду Crack
19% (40)
С удовольствием если будут деньги
23% (48)
Нет, потомучто не знаю как заплатить
6% (12)
Нет, принципиально!
17% (37)
Нет, у меня никогда нет денег!
14% (29)
Да, если разработчик будет русским
9% (20)
Нет, не хочу париться с оплатой
13% (27)

Проголосовало: 213
Встpечаются два интернетчика:
- Жена меня нy совсем достала! Все вpемя шипит на меня!
- А может, она сконнектиться пытается?
Рейтинг: 9/10 (1)
Посмотреть все анекдоты