unit YuleLevinson;

interface

uses
  Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs,
  Basic1, TeeProcs, TeEngine, Chart, StdCtrls, ExtCtrls, MtxVec,
  Toeplitz, Series, ComCtrls,Math387,MtxVecEdit, Buttons;

type
  TYuleLev = class(TBasicForm1)
    Series1: TFastLineSeries;
    Series2: TFastLineSeries;
    Panel3: TPanel;
    RadioGroup1: TRadioGroup;
    TrackBar1: TTrackBar;
    Label1: TLabel;
    Label2: TLabel;
    Label3: TLabel;
    Label4: TLabel;
    TrackBar2: TTrackBar;
    Label5: TLabel;
    Label6: TLabel;
    Button1: TBitBtn;
    procedure FormCreate(Sender: TObject);
    procedure FormDestroy(Sender: TObject);
    procedure Button1Click(Sender: TObject);
    procedure TrackBar1Change(Sender: TObject);
    procedure TrackBar2Change(Sender: TObject);
    procedure RadioGroup1Click(Sender: TObject);
  private
    x,y,spec,corr: TVec;
    CorrLen : Integer;
    ZeroPadding : Integer;
    LPCCoef : Integer;
    procedure bthHelpClick(Sender: TObject);
    { Private declarations }
  public
    { Public declarations }
  end;

var
  YuleLev: TYuleLev;

implementation

Uses MtxVecTee;

{$R *.DFM}

const LEVINSONTEHELP = 1700;

procedure TYuleLev.bthHelpClick(Sender: TObject);
begin
     Application.HelpContext(LEVINSONTEHELP);
end;

procedure TYuleLev.FormCreate(Sender: TObject);
begin
  inherited;
  With RichEdit1.Lines do
  begin
    Clear;
    Add('YuleWalker autoregressive spectra uses Levinson '
      + 'Durbin recursion to solve a toeplitz systems of '
      + 'linear equations taking only O(n2) operations '
      + 'instead of O(n3) as required by LUSolve. The chart '
      + 'compares FFT and YuleWalker AR. The corrLen defines '
      + 'the number of samples on which the Autocorrelation '
      + 'is performed and LPCCoef defines the number of '
      + 'computed autocorrelation coefficients. The method '
      + 'uses biased autocorrelation. FFT uses only LPCoef '
      + 'parameter to determine the number of sample to '
      + 'include. It then rounds LPCCoef to the nearest '
      + 'power of two. FFT uses no windowing.');
    Add('Zoom in on a chart (left-click and drag mouse '
      + 'over the chart) to see differences. Please note '
      + 'that it takes less then 10ms to compute a 32000 '
      + 'point FFT on P366.');
  end;

  CreateIt(x,y,spec,corr);
  RadioGroup1Click(RadioGroup1);
  TrackBar1Change(TrackBar1);
  TrackBar2Change(TrackBar2);
  with TButton(Self.Owner.FindComponent('btnHelp')) do
  begin
       OnClick := bthHelpClick;
       Enabled := true;
  end;
end;

procedure TYuleLev.FormDestroy(Sender: TObject);
begin
  FreeIt(x,y,spec,corr);
  TButton(Self.Owner.FindComponent('btnHelp')).OnClick := nil;
  inherited;
end;

procedure TYuleLev.Button1Click(Sender: TObject);
var DownY: TVec;
    Step : TSample;
begin
  Series1.Clear;
  Series2.Clear;
  Screen.Cursor := crHourGlass;
  TimeCheck := GetTickCount;
  { Levinson Yule Walker }
  y.LoadFromFile('FFTData.vec'); {Load signal}
  y.Resize(CorrLen);
  corr.AutoCorrBiased(y,LPCCoef);  {auto correlation}
  Levinson(corr, y);            {Levinson recursion}
  y.Resize(LargestExp2(y.Length*ZeroPadding),True); {zero padder}
  x.FFT(y);
  spec.Mag(x);
  spec.Inv(EPS);  {thresh bottom at 0.00001 before inverse}
  spec.Log10;
  TimeElapsed := GetTickCount - TimeCheck;
  Label5.Caption := 'Time needed for Levinson YW : '+ IntToStr(TimeElapsed)+' ms';
  if DownSize then
  begin
       CreateIt(DownY);
       try
          DownY.PixelDownSample(Chart1.ChartWidth,Spec);
          Step := 0.5*Trunc(Spec.Length/Chart1.ChartWidth);
          DrawValues(DownY,Series1,Step);
       finally
          FreeIt(DownY);
       end;
  end else DrawValues(spec,Series1);

  TimeCheck := GetTickCount;
  { "Regular" FFT }
  y.LoadFromFile('FFTData.vec'); {Load signal}
  y.Resize(LPCCoef);
  y.Resize(LargestExp2(LPCCoef)*ZeroPadding,True);
  x.FFT(y);
  spec.Mag(x);
  spec.ThreshBottom(0.00001); {-100 dB}
  spec.Log10;
  TimeElapsed := GetTickCount - TimeCheck;
  Label6.Caption := 'Time needed for FFT : '+ IntToStr(TimeElapsed)+' ms';
  if DownSize then
  begin
       CreateIt(DownY);
       try
          DownY.PixelDownSample(Chart1.ChartWidth,Spec);
          Step := 0.5*Trunc(Spec.Length/Chart1.ChartWidth);
          DrawValues(DownY,Series2,Step);
       finally
          FreeIt(DownY);
       end;
  end else DrawValues(spec,Series2);

  Chart1.Repaint;
  Screen.Cursor := crDefault;

end;

procedure TYuleLev.TrackBar1Change(Sender: TObject);
begin
  CorrLen := TTrackBar(Sender).Position;
  Label2.Caption := IntToStr(CorrLen);
  TrackBar2.Max := CorrLen;
  TrackBar2.Frequency := CorrLen div 20;
  TrackBar2.Position := CorrLen div 2;
  TrackBar2Change(TrackBar2);
end;

procedure TYuleLev.TrackBar2Change(Sender: TObject);
begin
  LPCCoef := TTrackBar(Sender).Position;
  Label4.Caption := IntToStr(LPCCoef);
end;

procedure TYuleLev.RadioGroup1Click(Sender: TObject);
begin
  case TRadioGroup(Sender).ItemIndex of
  0    : ZeroPadding := 1;
  1    : ZeroPadding := 2;
  2    : ZeroPadding := 4;
  3    : ZeroPadding := 8;
  4    : ZeroPadding := 16;
  end;
end;

initialization
   RegisterClass(TYuleLev);

end.
  