program testapp;

{$APPTYPE CONSOLE}

uses

  SysUtils,

  fftreal in 'fftreal.pas',

  Math,

  Windows;

var

   nbr_points  : longint;

   x, f        : pflt_array;

   fft         : TFFTReal;

   i           : longint;

   PI          : double;

   areal, img  : double;

   f_abs       : double;

   buffer_size : longint;

   nbr_tests   : longint;

   time0, time1, time2 : int64;

   timereso    : int64;

   offset      : longint;

   t0, t1      : double;

   nbr_s_chn   : longint;

   tempp1, tempp2 : pflt_array;

begin

  (*______________________________________________

   *

   * Exactness test

   *______________________________________________

   *)

  WriteLn('Accuracy test:');

  WriteLn;

  nbr_points := 16;       // Power of 2

  GetMem(x, nbr_points * sizeof_flt);

  GetMem(f, nbr_points * sizeof_flt);

  fft := TFFTReal.Create(nbr_points);    // FFT object initialized here

  // Test signal

  PI := ArcTan(1) * 4;

  for i := 0 to nbr_points-1 do

  begin

    x^[i] := -1 + sin (3*2*PI*i/nbr_points)

                + cos (5*2*PI*i/nbr_points) * 2

                - sin (7*2*PI*i/nbr_points) * 3

                + cos (8*2*PI*i/nbr_points) * 5;

  end;

  // Compute FFT and IFFT

  fft.do_fft(f, x);

  fft.do_ifft(f, x);

  fft.rescale(x);

  // Display the result

  WriteLn('FFT:');

  for i := 0 to nbr_points div 2 do

  begin

    areal := f^[i];

    if (i > 0) and (i < nbr_points div 2) then

      img := f^[i + nbr_points div 2]

    else

      img := 0;

    f_abs := Sqrt(areal * areal + img * img);

    WriteLn(Format('%5d: %12.6f %12.6f (%12.6f)', [i, areal, img, f_abs]));

  end;

  WriteLn;

  WriteLn('IFFT:');

  for i := 0 to nbr_points-1 do

    WriteLn(Format('%5d: %f', [i, x^[i]]));

  WriteLn;

  FreeMem(x);

  FreeMem(f);

  fft.Free;

  (*______________________________________________

   *

   * Speed test

   *______________________________________________

   *)

  WriteLn('Speed test:');

  WriteLn('Please wait...');

  WriteLn;

  nbr_points := 1024;	          // Power of 2

  buffer_size := 256*nbr_points;  // Number of flt_t (float or double)

  nbr_tests := 10000;

  assert(nbr_points <= buffer_size);

  GetMem(x, buffer_size * sizeof_flt);

  GetMem(f, buffer_size * sizeof_flt);

  fft := TFFTReal.Create(nbr_points);					// FFT object initialized here

  // Test signal: noise

  for i := 0 to nbr_points-1 do

    x^[i] := Random($7fff) - ($7fff shr 1);

  // timing

  QueryPerformanceFrequency(timereso);

  QueryPerformanceCounter(time0);

  for i := 0 to nbr_tests-1 do

  begin

    offset := (i * nbr_points) and (buffer_size - 1);

    tempp1 := f;

    inc(tempp1, offset);

    tempp2 := x;

    inc(tempp2, offset);

    fft.do_fft(tempp1, tempp2);

  end;

  QueryPerformanceCounter(time1);

  for i := 0 to nbr_tests-1 do

  begin

    offset := (i * nbr_points) and (buffer_size - 1);

    tempp1 := f;

    inc(tempp1, offset);

    tempp2 := x;

    inc(tempp2, offset);

    fft.do_ifft(tempp1, tempp2);

    fft.rescale(x);

  end;

  QueryPerformanceCounter(time2);

  t0 := ((time1-time0) / timereso) / nbr_tests;

  t1 := ((time2-time1) / timereso) / nbr_tests;

  WriteLn(Format('%d-points FFT           : %.0f us.', [nbr_points, t0 * 1000000]));

  WriteLn(Format('%d-points IFFT + scaling: %.0f us.', [nbr_points, t1 * 1000000]));

  nbr_s_chn := Floor(nbr_points / ((t0 + t1) * 44100 * 2));

  WriteLn(Format('Peak performance: FFT+IFFT on %d mono channels at 44.1 KHz (with overlapping)', [nbr_s_chn]));

  WriteLn;

  FreeMem(x);

  FreeMem(f);

  fft.Free;

  WriteLn('Press [Return] key to terminate...');

  ReadLn;

end.