unit UDFT;

interface

uses math,windows,sysutils;

type
	TArray = array of double;

	TRealValues = class
	public
		values : TArray;
		constructor Create(N : integer);
		destructor  free;
	end;

	TComplexValues = class
	public
		module : TRealValues;
    phase : TRealValues;
		im : TRealValues;
		re : TRealValues;
		constructor Create(N : integer);
		destructor  free;
		procedure computeModule;
    procedure computePhase(seuil : double);
	end;

	TDFT = class(TObject)
	private
		wcos : TRealValues;
		wsin : TRealValues;

		N : integer;
		reverseBits : array of integer;
		reverseBitsDone : array of boolean;

    seuilModule : double;
    useFFTAlgo : boolean;
    //Calcul TFD
    function		computeDFT(x : TRealValues): TComplexValues;
    //Calcul FFT
		function		computeFFT(x : TRealValues): TComplexValues;
    //Calcul FFT
		function    computeFFT2(x : TRealValues): TComplexValues;
	public
		constructor Create(N : integer; seuilModule : double; fft : boolean);
		destructor  free;
    function    compute(x : TRealValues): TComplexValues;
	end;

implementation

uses unit1;

constructor TRealValues.Create(N : integer);
begin
	inherited Create;
	setLength(values,N);
end;

destructor TRealValues.free;
begin
	SetLength(values,0);
	values := nil;

	//inherited free;
end;

constructor TComplexValues.Create(N : integer);
begin
	inherited Create;

	im := TRealValues.Create(N);
	re := TRealValues.Create(N);
	module:= TRealValues.Create(N);
  phase := TRealValues.Create(N);
end;

destructor TComplexValues.free;
begin
	im.free;
	re.free;
	module.free;
  phase.free;
  
	//inherited free;
end;

procedure TComplexValues.computeModule;
var i : integer;
begin
	for i := 0 to high(module.values)
	do module.values[i] := sqrt(re.values[i]*re.values[i] + im.values[i]*im.values[i])/(high(module.values) + 1);
end;

procedure TComplexValues.computePhase(seuil : double);
var i : integer;
begin
	for i := 0 to high(phase.values) do
  begin
    if(module.values[i] < seuil)
    then phase.values[i] := 0.0
    else begin
      if re.values[i] <> 0.0
      then phase.values[i] :=  arctan2(im.values[i],re.values[i])
      else
      begin
        if im.values[i] > 0
        then phase.values[i] := pi/2
        else phase.values[i] := -pi/2;
      end;
    end;
  end;
end;

constructor TDFT.Create(N : integer;seuilModule : double; fft:boolean);
var theta,S,C : double;
		k,l,nbbits : integer;

		function reverseBit(M : integer) : integer;
		var b,j : integer;
		begin
			nbbits := round(ln(N)/ln(2));
      b := 0;
      for j := 0 to nbbits - 1 do
      begin
        b := b + round(power(2,nbbits-1-j)*(M mod 2));
        M := M div 2;
      end;
			result := b;
		end;

begin
	inherited Create;

  Self.seuilModule := seuilModule;

  useFFTAlgo := fft;

  if fft=true then
  begin

    nbbits := trunc(ln(N)/ln(2));

    //Zero padding
    if abs(N - power(2,nbbits))>0
    then N := trunc(power(2,nbbits+1));

    Self.N := N;

    wcos := TRealValues.Create(N div 2);
    wsin := TRealValues.Create(N div 2);

    theta := 2*4*arctan(1)/N;
    S := sin(theta);
    C := 1 - 2*sin(theta/2)*sin(theta/2);

    wcos.values[0] := 1; wsin.values[0] := 0;
    for k := 0 to N div 8 - 2 do
    begin
      wcos.values[k + 1] := C*wcos.values[k] - S*wsin.values[k];
      wsin.values[k + 1] := S*wcos.values[k] + C*wsin.values[k];
    end;

    l := N div 8;
    wcos.values[l] := sqrt(2)/2; wsin.values[l] := wcos.values[l];
    for k := 1 to N div 8 - 1 do
    begin
      wcos.values[l + k] := wsin.values[l - k];
      wsin.values[l + k] := wcos.values[l - k];
    end;

    l := N div 4;
    wcos.values[l] := 0; wsin.values[l] := 1;
    for k := 1 to N div 4 - 1 do
    begin
      wcos.values[l + k] := -wcos.values[l - k];
      wsin.values[l + k] := wsin.values[l - k];
    end;

    SetLength(reverseBits,N);
    SetLength(reverseBitsDone,N);

    for k := 0 to N - 1 do
    begin
      reverseBits[k] := reverseBit(k);
      if reverseBits[k] = k
      then self.reverseBitsDone[k] := true
      else self.reverseBitsDone[k] := false;
    end;

  end else Self.N := N;
end;

destructor TDFT.free;
begin
  if useFFTAlgo = true then
  begin
    wcos.free;
    wsin.free;

    SetLength(reverseBits,0);
    SetLength(reverseBitsDone,0);
    reverseBits :=  nil;
    reverseBitsDone := nil;
  end;
	//inherited free;
end;

function TDFT.computeFFT(x : TRealValues) : TComplexValues;
var fft : TComplexValues;
		h,i,j,k,l,m,nbBits : integer;
		imTemp,reTemp : double;
		f,c1,c2 : int64; 
begin
	fft := TComplexValues.Create(N);
	for i := 0 to high(x.values) do fft.re.values[i] := x.values[i];

  QueryPerformanceFrequency(f);
	QueryPerformanceCounter(c1);

	nbBits := round(ln(N)/ln(2));

	for j := 1 to nbBits do
	begin
		i := round(N/power(2,j));
		h := round(power(2,j-1));
		for m := 1 to h do
		begin
			l := (m-1)*i*2;
			for k := 0 to i - 1 do
			begin
				reTemp := fft.re.values[k + l] - fft.re.values[k + l + i];
				imTemp := fft.im.values[k + l] - fft.im.values[k + l + i];

				fft.re.values[k + l] := fft.re.values[k+l] + fft.re.values[k + l + i];
				fft.im.values[k + l] := fft.im.values[k+l] + fft.im.values[k + l + i];

				fft.re.values[k + l + i] := reTemp*wcos.values[k*h] - imTemp*wsin.values[k*h];
				fft.im.values[k + l + i] := reTemp*wsin.values[k*h] + imTemp*wcos.values[k*h];
			end; 
		end;
	end;
	
 	for i := 0 to N div 2 - 1 do
	begin
		if not reverseBitsDone[i] then
		begin
			reTemp := fft.re.values[i];
			imTemp := fft.im.values[i];
			fft.re.values[i] := fft.re.values[reverseBits[i]];
			fft.im.values[i] := fft.im.values[reverseBits[i]];
			fft.re.values[reverseBits[i]] := reTemp;
			fft.im.values[reverseBits[i]] := imTemp;
			reTemp := fft.re.values[N-i-1];
			imTemp := fft.im.values[N-i-1];
			fft.re.values[N-i-1] := fft.re.values[reverseBits[N-i-1]];
			fft.im.values[N-i-1] := fft.im.values[reverseBits[N-i-1]];
			fft.re.values[reverseBits[N-i-1]] := reTemp;
			fft.im.values[reverseBits[N-i-1]] := imTemp;
			reverseBitsDone[reverseBits[i]] := true;
			reverseBitsDone[reverseBits[N-i-1]] := true;
		end;
	end;

	fft.computeModule;
  fft.computePhase(seuilModule);

  QueryPerformanceCounter(c2);
  form1.Caption := form1.Caption  + 'ma methode :' + floattostr(65*(c2-c1)/f);
	Result := fft;

end;

function TDFT.computeFFT2(x : TRealValues) : TComplexValues;
var fft,y : TComplexValues;
		pairsInGroup,numOfGroups,distance,i,j,k,l,jFirst,jLast,jTwiddle : integer;
		wre,wim,tempRe,tempIm : double;
		notSwitchInput : boolean;
		f,c1,c2 : int64;
begin
	fft := TComplexValues.Create(N);
	y  := TComplexValues.Create(N);
	for i := 0 to high(x.values)	do y.re.values[i] := x.values[i];

	QueryPerformanceFrequency(f);
	QueryPerformanceCounter(c1);

	pairsInGroup := N div 2;
	numOfGroups := 1;
	distance := N div 2;
	notSwitchInput := true;

	while numOfGroups < N do
	begin

		if notSwitchInput then
		begin
			l := 0;

			for k := 0 to numOfGroups - 1 do
			begin
				jFirst := 2*k*pairsInGroup;
				jLast := jFirst + pairsInGroup - 1;
				jTwiddle := k*pairsInGroup;
				wre := wcos.values[jTwiddle];
				wim := wsin.values[jTwiddle];
				for j := jFirst to jLast do
				begin
					tempRe := wre*y.re.values[j + distance] - wim*y.im.values[j + distance];
					tempIm := wre*y.im.values[j + distance] + wim*y.re.values[j + distance];

					fft.re.values[l] := y.re.values[j] + tempRe;
					fft.im.values[l] := y.im.values[j] + tempIm;
					fft.re.values[l + N div 2] := y.re.values[j] - tempRe;
					fft.im.values[l + N div 2] := y.im.values[j] - tempIm;

					inc(l);
               
				end;
			end;
			notSwitchInput := false;
		end else
		begin
      l := 0;

			for k := 0 to numOfGroups - 1 do
			begin
				jFirst := 2*k*pairsInGroup;
				jLast := jFirst + pairsInGroup - 1;
				jTwiddle := k*pairsInGroup;
				wre := wcos.values[jTwiddle];
				wim := wsin.values[jTwiddle];
				for j := jFirst to jLast do
				begin
					tempRe := wre*fft.re.values[j+ distance] - wim*fft.im.values[j + distance];
					tempIm := wre*fft.im.values[j+ distance] + wim*fft.re.values[j + distance];

					y.re.values[l] := fft.re.values[j] + tempRe;
					y.im.values[l] := fft.im.values[j] + tempIm;
					y.re.values[l + N div 2] := fft.re.values[j] - tempRe;
					y.im.values[l + N div 2] := fft.im.values[j] - tempIm;

					inc(l);
				end;
			end;
			notSwitchInput := true;
		end;
		pairsInGroup := pairsInGroup div 2;
		numOfGroups := numOfGroups*2;
		distance := distance div 2;

	end;

	fft.computeModule;
  fft.computePhase(seuilModule);

	QueryPerformanceCounter(c2);

	form1.Caption := 'DITNN : ' + floattostr(65*(c2-c1)/f) + ' - ';

	Result := fft;
end;



function TDFT.computeDFT(x : TRealValues) : TComplexValues;
var fft : TComplexValues;
    k,l : integer;
		f,c1,c2 : int64;
begin
	fft := TComplexValues.Create(N);

	QueryPerformanceFrequency(f);
	QueryPerformanceCounter(c1);

  for k := 0 to N - 1 do
  begin
    for l := 0 to N - 1 do
    begin
      fft.re.values[k] := fft.re.values[k] + x.values[l]* cos(2*pi*k*l/N);
      fft.im.values[k] := fft.im.values[k] - x.values[l]* sin(2*pi*k*l/N);
    end;
  end;

	fft.computeModule;
  fft.computePhase(seuilModule);

	QueryPerformanceCounter(c2);

	//Caption := 'DFT time : ' + floattostr(65*(c2-c1)/f) + ' - ';

	Result := fft;
end;

function TDFT.compute(x : TRealValues) : TComplexValues;
begin
  if Self.useFFTAlgo
  then result := computeFFT2(x)
  else result := computeDFT(x);
end;


end.
