{ Unit FFTs

  This unit provides a forward and inverse FFT pascal implementation
  for complex number series.

  The formal definition of the complex DFT is:
    y[k] = sum(x[m]*exp(-i*2*pi*k*m/n), m = 0..n-1), k = 0..n-1

  Copyright: Nils Haeck M.Sc. (email: n.haeck@simdesign.nl)
  For more information visit http://www.simdesign.nl
  Original date of publication: 10 Mar 2003

  This unit requires these other units:
  - Complexs: Complex number unit
  - Types:    Additional mathematical variable types
  - SysUtils: Delphi system utilities

  ****************************************************************

  The contents of this file are subject to the Mozilla Public
  License Version 1.1 (the "License"); you may not use this file
  except in compliance with the License. You may obtain a copy of
  the License at:
  http://www.mozilla.org/MPL/

  Software distributed under the License is distributed on an
  "AS IS" basis, WITHOUT WARRANTY OF ANY KIND, either express or
  implied. See the License for the specific language governing
  rights and limitations under the License.
}
unit FFTs;

interface

uses
  Complexs, Types, SysUtils;

const

  cMaxPrimeFactor     = 1021;
  cMaxPrimeFactorDiv2 = (cMaxPrimeFactor + 1) div 2;
  cMaxFactorCount     = 20;

resourcestring

  sErrPrimeTooLarge = 'Prime factor for FFT length too large. Change value for cMaxPrimeFactor in FFTs unit';

{ ForwardFFT:
  Perform a complex 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 complex DFT.

  The largest prime factor in Count should be less or equal to cMaxPrimeFactor.

  The remaining factors are handled by optimised partial FFT code, that can be
  found in the FFT_X procedures

  Inputs:
    Source: this can be any zero-based array type of TComplex
    Count: The number of elements in the array.

  Outputs:
    Dest: this can be any zero-based array type of TComplex, and will contain
      the FFT transformed data (frequency spectrum). Source may be equal to
      Dest. In this case, the original series will be overwritten with the new
      fourier-transformed series.
}
procedure ForwardFFT(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. This is based
  on the forward FFT with some additional customisation. The result of a forward
  FFT followed by an inverse FFT should yield the same data, except for rounding
  errors.
}
procedure InverseFFT(const Source: array of TComplex; var Dest: array of TComplex; Count: integer);

implementation

const
  // Some helper constants for the FFT optimisations
  c31: TFloat = -1.5000000000000E+00; //  cos(2*pi / 3) - 1;
  c32: TFloat =  8.6602540378444E-01; //  sin(2*pi / 3);

  u5:  TFloat =  1.2566370614359E+00; //  2*pi / 5;
  c51: TFloat = -1.2500000000000E+00; // (cos(u5) + cos(2*u5))/2 - 1;
  c52: TFloat =  5.5901699437495E-01; // (cos(u5) - cos(2*u5))/2;
  c53: TFloat = -9.5105651629515E-01; //- sin(u5);
  c54: TFloat = -1.5388417685876E+00; //-(sin(u5) + sin(2*u5));
  c55: TFloat =  3.6327126400268E-01; // (sin(u5) - sin(2*u5));
  c8:  TFloat =  7.0710678118655E-01; //  1 / sqrt(2);

type
  // Base 1 and Base 0 arrays
  TIdx0FactorArray = array[0..cMaxFactorCount] of integer;
  TIdx1FactorArray = array[1..cMaxFactorCount] of integer;

// Factorise the series with length Count into FactorCount factors, stored in Fact
procedure Factorize(Count: integer; var FactorCount: integer; var Fact: TIdx1FactorArray);
var
  i: integer;
  Factors: TIdx1FactorArray;
const
  // Define specific FFT lengths (radices) that we can process with optimised routines
  cRadixCount = 6;
  cRadices: array[1..6] of integer =
    (2, 3, 4, 5, 8, 10);
begin

  if Count = 1 then begin
    FactorCount := 1;
    Factors[1]  := 1;
  end else begin
    FactorCount := 0;
  end;

  // Factorise the original series length Count into known factors and rest value
  i := cRadixCount;
  while (Count > 1) AND (i > 0) do begin
    if Count mod cRadices[i] = 0 then begin
      Count := Count div cRadices[i];
      inc(FactorCount);
      Factors[FactorCount] := cRadices[i];
    end else
      dec(i);
  end;

  // substitute factors 2*8 with more optimal 4*4
  if Factors[FactorCount] = 2 then begin
    i := FactorCount - 1;
    while (i > 0) AND (Factors[i] <> 8) do
      dec(i);
    if i > 0 then begin
      Factors[FactorCount] := 4;
      Factors[i] := 4;
    end;
  end;

  // Analyse the rest value and see if it can be factored in primes
  if Count > 1 then begin
    for i := 2 to trunc(sqrt(Count)) do begin
      while Count mod i = 0 do begin
        Count := Count div i;
        inc(FactorCount);
        Factors[FactorCount] := i;
      end;
    end;

    if (Count > 1) then begin
      inc(FactorCount);
      Factors[FactorCount] := Count;
    end;
  end;

  // Reverse factors so that primes are first
  for i := 1 to FactorCount do
    Fact[i] := Factors[FactorCount - i + 1];

end;

{ Reorder the series in X to a permuted sequence in Y so that the later step can
  be done in place, and the final FFT result is in correct order.
  The series X and Y must be different series!
}
procedure ReorderSeries(Count: integer; var Factors: TIdx1FactorArray; var Remain: TIdx0FactorArray;
  const X: array of TComplex; var Y: array of TComplex);
var
  i, j, k: integer;
  Counts: TIdx1FactorArray;
begin
  FillChar(Counts, SizeOf(Counts), 0);

  k := 0;
  for i := 0 to Count - 2 do begin
    Y[i] := X[k];
    j := 1;
    k := k + Remain[j];
    Counts[1] := Counts[1] + 1;
    while Counts[j] >= Factors[j] do begin
      Counts[j] := 0;
      k := k - Remain[j - 1] + Remain[j + 1];
      inc(j);
      inc(Counts[j]);
    end;
  end;

  Y[Count - 1] := X[Count - 1];
end;

procedure FFT_2(var Z: array of TComplex);
var
  T1: TComplex;
begin
  T1   := ComplexAdd(Z[0], Z[1]);
  Z[1] := ComplexSub(Z[0], Z[1]);
  Z[0] := T1;
end;

procedure FFT_3(var Z: array of TComplex);
var
  T1, M1, M2, S1: TComplex;
begin
  T1   := ComplexAdd(Z[1], Z[2]);
  Z[0] := ComplexAdd(Z[0], T1);
  M1   := ComplexScl(c31, T1);
  M2.Re := c32 * (Z[1].Im - Z[2].Im);
  M2.Im := c32 * (Z[2].Re - Z[1].Re);
  S1   := ComplexAdd(Z[0], M1);
  Z[1] := ComplexAdd(S1, M2);
  Z[2] := ComplexSub(S1, M2);
end;

procedure FFT_4(var Z: array of TComplex);
var
  T1, T2, M2, M3: TComplex;
begin
  T1 := ComplexAdd(Z[0], Z[2]);
  T2 := ComplexAdd(Z[1], Z[3]);

  M2 := ComplexSub(Z[0], Z[2]);
  M3.Re := Z[1].Im - Z[3].Im;
  M3.Im := Z[3].Re - Z[1].Re;

  Z[0] := ComplexAdd(T1, T2);
  Z[2] := ComplexSub(T1, T2);
  Z[1] := ComplexAdd(M2, M3);
  Z[3] := ComplexSub(M2, M3);
end;

procedure FFT_5(var Z: array of TComplex);
var
  T1, T2, T3, T4, T5: TComplex;
  M1, M2, M3, M4, M5: TComplex;
  S1, S2, S3, S4, S5: TComplex;
begin
  T1 := ComplexAdd(Z[1], Z[4]);
  T2 := ComplexAdd(Z[2], Z[3]);
  T3 := ComplexSub(Z[1], Z[4]);
  T4 := ComplexSub(Z[3], Z[2]);

  T5   := ComplexAdd(T1, T2);
  Z[0] := ComplexAdd(Z[0], T5);
  M1   := ComplexScl(c51, T5);
  M2   := ComplexScl(c52, ComplexSub(T1, T2));

  M3.Re := -c53 * (T3.Im + T4.Im);
  M3.Im :=  c53 * (T3.Re + T4.Re);
  M4.Re := -c54 * T4.Im;
  M4.Im :=  c54 * T4.Re;
  M5.Re := -c55 * T3.Im;
  M5.Im :=  c55 * T3.Re;

  S3 := ComplexSub(M3, M4);
  S5 := ComplexAdd(M3, M5);;
  S1 := ComplexAdd(Z[0], M1);
  S2 := ComplexAdd(S1, M2);
  S4 := ComplexSub(S1, M2);

  Z[1] := ComplexAdd(S2, S3);
  Z[2] := ComplexAdd(S4, S5);
  Z[3] := ComplexSub(S4, S5);
  Z[4] := ComplexSub(S2, S3);
end;

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