FFTW

program new;
{$I SRL/OSR.simba}

const
 FFTW_FORWARD        = -1;
 FFTW_BACKWARD        = 1;

 FFTW_MEASURE         = 0;
 FFTW_DESTROY_INPUT   = 1;   {1U << 0}
 FFTW_UNALIGNED       = 2;   {1U << 1}
 FFTW_CONSERVE_MEMORY = 4;   {1U << 2}
 FFTW_EXHAUSTIVE      = 8;   {1U << 3} {NO_EXHAUSTIVE is default }
 FFTW_PRESERVE_INPUT  = 16;  {1U << 4} {cancels FFTW_DESTROY_INPUT}
 FFTW_PATIENT         = 32;  {1U << 5} {IMPATIENT is default }
 FFTW_ESTIMATE        = 64;  {1U << 6}


type
  PComplex = ^Complex;
  Complex = packed record
    RE,IM: Single;
  end;
  TComplexArray = array of Complex;

  FFTW_PLAN  = type Pointer;

  LP_FFTW_MALLOC    = function(n: Int32): Pointer;
  LP_FFTW_FREE      = procedure(p: Pointer);
  LP_FFTW_PLAN_DFT1 = function(n: Int32; inData, outData: PComplex; sign: Int32; flags: UInt32): FFTW_PLAN;
  LP_FFTW_EXEC_DFT1 = procedure(plan: FFTW_PLAN; inData, outData: PComplex);
  LP_FFTW_EXEC      = procedure(plan: FFTW_PLAN);

  FN_FFTW_MALLOC    = native(LP_FFTW_MALLOC, ffi_cdecl);
  FN_FFTW_FREE      = native(LP_FFTW_FREE, ffi_cdecl);
  FN_FFTW_PLAN_DFT1 = native(LP_FFTW_PLAN_DFT1, ffi_cdecl);
  FN_FFTW_EXEC_DFT1 = native(LP_FFTW_EXEC_DFT1, ffi_cdecl);
  FN_FFTW_EXEC      = native(LP_FFTW_EXEC, ffi_cdecl);

var
  FFTW: THandle;
  n_fftw_malloc:    FN_FFTW_MALLOC;
  n_fftw_free:      FN_FFTW_FREE;
  n_fftw_plan_dft1: FN_FFTW_PLAN_DFT1;
  n_fftw_exec:      FN_FFTW_EXEC;

  fftw_malloc:    LP_FFTW_MALLOC;
  fftw_free:      LP_FFTW_FREE;
  fftw_plan_dft1: LP_FFTW_PLAN_DFT1;
  fftw_exec:      LP_FFTW_EXEC;


function ToString(x: Complex): String; override;
begin
  if Round(x.im,3) >= 0 then
    Result := ToStr(Round(x.re,3))+'+'+ToStr(Round(x.im,3))+'i'
  else
    Result := ToStr(Round(x.re,3))+ToStr(Round(x.im,3))+'i'
end;

function AsComplex(x: array of Extended): TComplexArray;
var i:Int32;
begin
  SetLength(Result, Length(x));
  for i:=0 to High(x) do
  begin
    Result[i].re := x[i];
    Result[i].im := 0;
  end;
end;


procedure LoadFFTW();
begin
  FFTW := LoadLibrary('libfftw3f-3.dll');

  n_fftw_malloc    := GetProcAddress(FFTW, 'fftwf_malloc');
  n_fftw_free      := GetProcAddress(FFTW, 'fftwf_free');
  n_fftw_plan_dft1 := GetProcAddress(FFTW, 'fftwf_plan_dft_1d');
  n_fftw_exec      := GetProcAddress(FFTW, 'fftwf_execute');

  fftw_malloc    := Lapify(n_fftw_malloc);
  fftw_free      := Lapify(n_fftw_free);
  fftw_plan_dft1 := Lapify(n_fftw_plan_dft1);
  fftw_exec      := Lapify(n_fftw_exec);
end;

var
  raw: TComplexArray;

  test, res: PComplex;
  plan: FFTW_PLAN;
  N: Int32;

  t: Double;
begin
  LoadFFTW();

  //raw  := AsComplex([1,1,1,1,0,0,0,0]);
  SetLength(raw, 1024*1024);

  N    := Length(raw);
  test := fftw_malloc(SizeOf(Complex) * N);
  res  := fftw_malloc(SizeOf(Complex) * N);

  t := PerformanceTimer();
  plan := fftw_plan_dft1(N, test, res, FFTW_FORWARD, FFTW_ESTIMATE);
  WriteLn(PerformanceTimer() - t, 'ms');

  t := PerformanceTimer();
  MemMove(raw[0], test^, SizeOf(Complex)*N);
  fftw_exec(plan);
  MemMove(res^, raw[0], SizeOf(Complex)*N);

  WriteLn(PerformanceTimer() - t, 'ms');

  FreeLibrary(FFTW);
end;