unit FFT;

interface

uses IOGraphics;

const N = 10;
      M = 3600;

type  {typy pro FFT a fazovou korelaci}
     Complex = record
                 Re : single;
                 Im : single;
               end;{record}

     Spectrum = array [1..N,1..N] of complex;
     BigSpectrum = array [1..M,1..M] of complex;
     FFTImage = array [1..N,1..N] of word;
     BigFFTImage = array [1..M,1..M] of word;

     pSpc     = ^Spectrum;
     pBigSpc  = ^BigSpectrum;

procedure FourierTransform(var Img : FFTImage; var Spc : Spectrum);
procedure InverseFourierTransform(var Spc : BigSpectrum; var Img : BigFFTImage);

implementation

var  b : BigSpectrum;

const FFTW_FORWARD  = -1;
      FFTW_BACKWARD = 1;
      FFTW_IN_PLACE = 8;
      FFTW_ESTIMATE = 0;
      FFTW_MEASURE  = 1;

function  fftw2d_create_plan( nx, ny, d, f : integer ) : pointer;
          cdecl external 'sfftw.dll';
procedure fftwnd_one( p, i, o : pointer );
          cdecl external 'sfftw.dll';


procedure FourierTransform(var Img : FFTImage; var Spc : Spectrum);
var sign : longint;
    i, j : integer;
    plan : pointer;
    pSpc : pointer;
begin

  sign:=1;
  for  i:=1 to N do
    for j:=1 to N do
      begin
        Spc[i,j].Re:=Img[i,j]/(N*N);
        Spc[i,j].Im:=0.0;
      end;

  pSpc := @Spc;

  plan := fftw2d_create_plan( N, N, FFTW_FORWARD, FFTW_ESTIMATE or FFTW_IN_PLACE );
  fftwnd_one( plan, pSpc, nil );
end;



procedure InverseFourierTransform(var Spc : BigSpectrum; var Img : BigFFTImage);
var sign : longint;
    i, j : integer;
    plan : pointer;
    pSpc : pointer;
begin
  b := Spc;
  pSpc := @b;
  plan := fftw2d_create_plan( M, M, FFTW_BACKWARD, FFTW_ESTIMATE or FFTW_IN_PLACE );
  fftwnd_one( plan, pSpc, nil );

  for  i:=1 to M do
    for j:=1 to M do Img[i,j] := Round(b[i,j].Re);
end;




end.
