//============================================================================
//                         FFT ユニット Ver.0.10
//							original by 奥村晴彦「C言語によるアルゴリズム入門」
//              Pascal Converted by C60
//============================================================================

// メモリリーク対策にTFFT.FreeMemory関数を追加。

unit _FFT_C60_single;

interface

uses SysUtils;

type
  PInteger = ^Integer;
  PSingle = ^Single;

  EFFTError = class(Exception);

  TFFT = class(TObject)
  private
    { Private 宣言 }
    last_n : Integer;
    sintbl : PSingle;
    bitrev : PInteger;
    procedure make_sintbl(n : Integer; sintbl : PSingle);
    procedure make_bitrev(n : Integer; bitrev : PInteger);
    public
    { Public 宣言 }
    constructor Create;
    procedure FreeMemory;
    function FFT(n : Integer; var x, y : array of Single) : Integer;
  end;

implementation


//***************************************************************************
// コンストラクタ
//***************************************************************************
constructor TFFT.Create;
begin
	last_n  := 0;
	bitrev := nil;
	sintbl := nil;
end;

procedure TFFT.FreeMemory;
begin
	last_n  := 0;
  if sintbl<>nil then FreeMem(sintbl);
  if bitrev<>nil then FreeMem(bitrev);
end;


procedure TFFT.make_sintbl(n : Integer; sintbl : PSingle);
var
	i, n2, n4, n8 : Integer;
	c, s, dc, ds, t : Single;
begin
	n2 := n div 2; n4 := n div 4; n8 := n div 8;
	t := sin(PI / n);
	dc := 2 * t * t;
	ds := sqrt(dc * (2 - dc));
	t := 2 * dc;

//	sintbl[n4] := 1; c := 1;
//	sintbl[0] := 0; s := 0;
	PSingle(Integer(sintbl) + n4*SizeOf(Single))^ := 1; c := 1;
	sintbl^ := 0; s := 0;

	for i := 1 to n8-1 do begin
  	c := c - dc; dc := dc + t * c;
		s := s + ds; ds := ds - t * s;
//		sintbl[i] := s;	sintbl[n4 - i] := c;
		PSingle(Integer(sintbl) + i*SizeOf(Single))^ := s;
		PSingle(Integer(sintbl) + (n4-i)*SizeOf(single))^ := c;
	end;

//	if(n8 <> 0) then sintbl[n8] := sqrt(0.5);
  if(n8 <> 0) then Psingle(Integer(sintbl) + n8*SizeOf(single))^ := sqrt(0.5);
	for i := 0 to n4-1 do begin
//		sintbl[n2 - i] := sintbl[i];
		Psingle(Integer(sintbl) + (n2-i)*SizeOf(single))^ :=
														    Psingle(Integer(sintbl) + i*SizeOf(single))^;
	end;

	for i := 0 to n2 + n4 - 1 do begin
//		sintbl[i + n2] := -sintbl[i];
		Psingle(Integer(sintbl) + (i+n2)*SizeOf(single))^ :=
														    -Psingle(Integer(sintbl) + i*SizeOf(single))^;
	end;
end;


procedure TFFT.make_bitrev(n : Integer; bitrev : PInteger);
var
	i, j, k, n2 : Integer;
begin
	n2 := n div 2;
	i := 0; j := 0;

	while(True) do begin
//		bitrev[i] := j;
		PInteger(Integer(bitrev) + i*SizeOf(Integer))^ := j;
    Inc(i);
		if(i >= n) then Break;
		k := n2;
		while(k <= j) do begin
			Dec(j, k); k := k div 2;
		end;
		Inc(j, k);
	end;
end;


function TFFT.FFT(n : Integer; var x, y : array of single) : Integer;
var
	i, j, k, ik, h, d, k2, n4, inverse : Integer;
	t, s, c, dx, dy : single;
begin
	if(n < 0) then begin
		n := -n;
		inverse := 1;
	end else begin
		inverse := 0;
	end;

	n4 := n div 4;
	if(n <> last_n) or (n = 0) then begin
		last_n := n;
		if(sintbl <> nil) then begin
			FreeMem(sintbl);
		end;
		if(bitrev <> nil) then begin
			FreeMem(bitrev);
		end;

		if(n = 0) then begin
			Result := 0;
			Exit;
    end;

		GetMem(sintbl, (n + n4) * SizeOf(single));
		GetMem(bitrev, n * SizeOf(Integer));

		if(sintbl = nil) or (bitrev = nil) then begin
	  	raise EFFTError.Create('メモリを取得できません');
		end;

		make_sintbl(n, sintbl);
		make_bitrev(n, bitrev);
	end;

	for i := 0 to n-1 do begin
//		j := bitrev[i];
		j := PInteger(Integer(bitrev) + i*SizeOf(Integer))^;
		if(i < j) then begin
			t := x[i];	x[i] := x[j];	x[j] := t;
			t := y[i];	y[i] := y[j];	y[j] := t;
		end;
	end;

  k := 1;
  while(k < n) do begin
		h := 0; k2 := k + k; d := n div k2;
		for j := 0 to k-1 do begin
//    	c := sintbl[h+n4];
			c := Psingle(Integer(sintbl) + (h+n4)*SizeOf(single))^;
			if(inverse <> 0) then begin
//      	s := -sintbl[h];
				s := -Psingle(Integer(sintbl) + h*SizeOf(single))^;
			end else begin
//      	s :=  sintbl[h];
				s :=  Psingle(Integer(sintbl) + h*SizeOf(single))^;
      end;

      i := j;
      while(i < n) do begin
      	ik := i + k;
				dx := s * y[ik] + c * x[ik];
				dy := c * y[ik] - s * x[ik];
				x[ik] := x[i] - dx; x[i] := x[i] + dx;
				y[ik] := y[i] - dy; y[i] := y[i] + dy;
				Inc(i, k2);
        end;
      Inc(h, d);
		end;
  	k := k + k;
  end;

	if (inverse = 0) then begin
		for i := 0 to n-1 do begin
			x[i] := x[i] / n; y[i] := y[i] / n;
		end;
	end;

	Result := 0;
end;

end.
