FFT_4(A);
  FFT_4(B);

  Gem     := c8 * (B[1].Re + B[1].Im);
  B[1].Im := c8 * (B[1].Im - B[1].Re);
  B[1].Re := Gem;
  Gem     := B[2].Im;
  B[2].Im :=-B[2].Re;
  B[2].Re := Gem;
  Gem     := c8 * (B[3].Im - B[3].Re);
  B[3].Im :=-c8 * (B[3].Re + B[3].Im);
  B[3].Re := Gem;

  Z[0] := ComplexAdd(A[0], B[0]); Z[4] := ComplexSub(A[0], B[0]);
  Z[1] := ComplexAdd(A[1], B[1]); Z[5] := ComplexSub(A[1], B[1]);
  Z[2] := ComplexAdd(A[2], B[2]); Z[6] := ComplexSub(A[2], B[2]);
  Z[3] := ComplexAdd(A[3], B[3]); Z[7] := ComplexSub(A[3], B[3]);
end;

procedure FFT_10(var Z: array of TComplex);
var
  A, B: array[0..4] of TComplex;
begin
   A[0] := Z[0];  B[0] := Z[5];
   A[1] := Z[2];  B[1] := Z[7];
   A[2] := Z[4];  B[2] := Z[9];
   A[3] := Z[6];  B[3] := Z[1];
   A[4] := Z[8];  B[4] := Z[3];

   FFT_5(A);
   FFT_5(B);

   Z[0] := ComplexAdd(A[0], B[0]); Z[5] := ComplexSub(A[0], B[0]);
   Z[6] := ComplexAdd(A[1], B[1]); Z[1] := ComplexSub(A[1], B[1]);
   Z[2] := ComplexAdd(A[2], B[2]); Z[7] := ComplexSub(A[2], B[2]);
   Z[8] := ComplexAdd(A[3], B[3]); Z[3] := ComplexSub(A[3], B[3]);
   Z[4] := ComplexAdd(A[4], B[4]); Z[9] := ComplexSub(A[4], B[4]);
end;

{
  Synthesize the FFT by taking the even factors and the odd factors multiplied by
  complex sinusoid
}
procedure SynthesizeFFT(Sofar, Radix, Remain: integer; var Y: array of TComplex);
var
  GroupOffset, DataOffset, Position: integer;
  GroupNo, DataNo, BlockNo, SynthNo: integer;
  Omega: double;
  S, CosSin: TComplex;
  Synth, Trig, Z: array[0..cMaxPrimeFactor - 1] of TComplex;

  // Local function
  procedure InitializeTrigonomials(Radix: integer);
  // Initialize trigonomial coefficients
  var
    i: integer;
    W: double;
    X: TComplex;
  begin
    W := 2 * pi / Radix;
    Trig[0] := Complex(1.0, 0.0);
    X := Complex(cos(W), -sin(W));
    Trig[1] := X;
    for i := 2 to Radix - 1 do
      Trig[i] := ComplexMul(X, Trig[i - 1]);
  end;

  // Local Function
  procedure FFT_Prime(Radix: integer);
  // This is the general DFT, which can't be made any faster by factoring because
  // Radix is a prime number
  var
    i, j, k, N, AMax: integer;
    Re, Im: TComplex;
    V, W: array[0..cMaxPrimeFactorDiv2 - 1] of TComplex;
  begin
    N := Radix;
    AMax := (N + 1) div 2;
    for j := 1 to AMax - 1 do begin
      V[j].Re := Z[j].Re + Z[n-j].Re;
      V[j].Im := Z[j].Im - Z[n-j].Im;
      W[j].Re := Z[j].Re - Z[n-j].Re;
      W[j].Im := Z[j].Im + Z[n-j].Im;
    end;

    for j := 1 to AMax - 1 do begin
      Z[j]   := Z[0];
      Z[N-j] := Z[0];
      k := j;
      for i := 1 to AMax - 1 do begin
        Re.Re := Trig[k].Re * V[i].Re;
        Im.Im := Trig[k].Im * V[i].Im;
        Re.im := Trig[k].Re * W[i].Im;
        Im.Re := Trig[k].Im * W[i].Re;

        Z[N-j].Re := Z[N-j].Re + Re.Re + Im.Im;
        Z[N-j].Im := Z[N-j].Im + Re.Im - Im.Re;
        Z[j].Re   := Z[j].Re   + Re.Re - Im.Im;
        Z[j].Im   := Z[j].Im   + Re.Im + Im.Re;

        k := k + j;
        if k >= N then
          k := k - N;
      end;
    end;

    for j := 1 to AMax - 1 do begin
      Z[0].Re := Z[0].Re + V[j].Re;
      Z[0].Im := Z[0].Im + W[j].Im;
    end;
  end;

// main
begin
  // Initialize trigonomial coefficients
  InitializeTrigonomials(Radix);

  Omega       := 2 * pi / (Sofar * Radix);
  CosSin      := Complex(cos(Omega), -sin(Omega));
  S           := Complex(1.0, 0.0);
  DataOffset  := 0;
  GroupOffset := 0;
  Position    := 0;

  for DataNo := 0 to Sofar - 1 do begin

    if Sofar > 1 then begin

      Synth[0] := Complex(1.0, 0.0);
      Synth[1] := S;
      for SynthNo := 2 to Radix - 1 do
        Synth[SynthNo] := ComplexMul(S, Synth[SynthNo - 1]);
      S := ComplexMul(CosSin, S);

    end;

    for GroupNo := 0 to Remain - 1 do begin

      if (Sofar > 1) AND (DataNo > 0) then begin

        Z[0]    := Y[Position];
        BlockNo := 1;
        repeat
          inc(Position, Sofar);
          Z[BlockNo] := ComplexMul(Synth[BlockNo], Y[Position]);
          inc(BlockNo);
        until BlockNo >= Radix;

      end else begin

        for BlockNo := 0 to Radix - 1 do begin
          Z[BlockNo] := Y[Position];
          inc(Position, Sofar);
        end;

      end;

      case Radix of
      2:  FFT_2(Z);
      3:  FFT_3(Z);
      4:  FFT_4(Z);
      5:  FFT_5(Z);
      8:  FFT_8(Z);
      10: FFT_10(Z);
      else
        // Any larger prime number than 5 (so 7, 11, 13, etc, up to cMaxPrimeFactor)
        FFT_Prime(Radix);
      end; //case

      Position := GroupOffset;
      for BlockNo := 0 to Radix - 1 do begin
        Y[Position] := Z[blockNo];
        Inc(Position, Sofar);
      end;
      GroupOffset := GroupOffset + Sofar * Radix;
      Position    := GroupOffset;
    end;
    inc(DataOffset);
    GroupOffset := DataOffset;
    Position    := DataOffset;
  end;
end;

procedure ForwardFFT(const Source: array of TComplex; var Dest: array of TComplex; Count: integer);
// Perform a FFT on the data in Source, put result in Dest. This routine works best
// for Count as a power of 2, but also works usually faster than DFT by factoring
// the series. Only in cases where Count is a prime number will this method be
// identical to regular DFT.
type
  PComplexArray = ^TComplexArray;
  TComplexArray = array[0..0] of TComplex;
var
  i: integer;
  FactorCount: integer;
  SofarRadix:  TIdx1FactorArray;
  ActualRadix: TIdx1FactorArray;
  RemainRadix: TIdx0FactorArray;
  TmpDest: PComplexArray;
begin
  if Count = 0 then exit;

  // Decompose the series with length Count into FactorCount factors in ActualRadix
  Factorize(Count, FactorCount, ActualRadix);

  // Check if our biggest prime factor is not too large
  if (ActualRadix[1] > cMaxPrimeFactor) then
    raise EMathError.Create(sErrPrimeTooLarge);

  // Setup Sofar and Remain tables
  RemainRadix[0] := Count;
  SofarRadix[1]  := 1;
  RemainRadix[1] := Count div ActualRadix[1];
  for i := 2 to FactorCount do begin
    SofarRadix[i]  := SofarRadix[i-1] * ActualRadix[i-1];
    RemainRadix[i] := RemainRadix[i-1] div ActualRadix[i];
  end;

  // Make temp copy if dest = source (otherwise the permute procedure will completely
  // ruin the structure
  if @Dest = @Source then begin
    GetMem(TmpDest, SizeOf(TComplex) * Count);;
    Move(Dest, TmpDest^, SizeOf(TComplex) * Count);
  end else begin
    TmpDest := @Dest;
  end;

  // Reorder the series so that the elements are already in the right place for
  // synthesis
  ReorderSeries(Count{, FactorCount}, ActualRadix, RemainRadix, Source, TmpDest^);

  // Free the temporary copy (if any)
  if @Dest = @Source then begin
    Move(TmpDest^, Dest, SizeOf(TComplex) * Count);
    FreeMem(TmpDest);
  end;

  // Synthesize each of the FFT factored series
  for i := 1 to FactorCount do
    SynthesizeFFT(SofarRadix[i], ActualRadix[i], RemainRadix[i], Dest);

end;

procedure InverseFFT(const Source: array of TComplex; var Dest: array of TComplex; Count: integer);
// Perform the inverse FFT on the Source data, and put result in Dest. It performs
// the forward FFT and then divides elements by N
var
  i: integer;
  S: TFloat;
  TmpSource: array of TComplex;
begin
  if Count = 0 then exit;

  // Since TmpSource is local, we do not have to free it again,
  // it will be freed automatically when out of scope
  SetLength(TmpSource, Count);

  // Create a copy with inverted imaginary part.
  for i := 0 to Count - 1 do
    with Source[i] do
      TmpSource[i] := Complex(Re, -Im);
  ForwardFFT(TmpSource, Dest, Count);

  // Scale by 1/Count, and inverse the imaginary part
  S := 1.0 / Count;
  for i := 0 to Count - 1 do begin
    Dest[i].Re :=   S * Dest[i].Re;
    Dest[i].Im := - S * Dest[i].Im;
  end;
end;

end.