Contributor: JACK BAKKER

{
> I'm Looking for a source for a FFT of the Soundblaster
> Input. Could somebody help me?

Do you know the basics of Fourier Transformation?
Not that you need to, for this routine, but it might be handy to know.
Anyhow, here's something I still didn't test:
}

Program FOURIER;

{ Real FFT with single sine look-up per pass }

Const
  Asize = 4096;           { Size of array goes here }
  PI2   = 1.570796327;    { PI/2 }
  F     = 16;             { Format constants }
  D     = 6;              {   "      " }

Type
  Xform = (Fwd,Inverse);  { Transform types }
  Xary  = Array[1..Asize] of Real;

Var
  I,j,N,modulo   : Integer;
  F1    : Text;  { File of Real;}
  X     : Xary;           { Array to transform }
  Inv   : Xform;          { Transform type - Forward or Inverse }
  w     : Real;           { Frequency of wave }

{*****************************************************************************}

Procedure Debug;          { Used to print intermediate results }

Var I3 : Integer;

Begin
  For I3 := 1 to N Do
    Writeln((I3-1):3,X[I3]:F:D);
End;    { Debug }

{*****************************************************************************}

Function Ibitr(j,nu : Integer) : Integer;
  { Function to bit invert the number j by nu bits }

Var
  i,j2,ib    : Integer;

Begin
  ib := 0;   { Default return integer }
  For i := 1 to nu do
  Begin
    j2 := j Div 2;          { Divide by 2 and compare lowest bits }
    { ib is doubled and bit 0 set to 1 if j is odd }
    ib := ib*2 + (j - 2*j2);
    j := j2;                { For next pass }
  End;  {For}
  ibitr := ib;              { Return bit inverted value }
End;    { Ibitr }

{*****************************************************************************}

Procedure Expand;

Var
  i,nn2,nx2    : Integer;

Begin
  nn2 := n div 2;
  nx2 := n + n;
  For i := 1 to nn2 do x[i+n] := x[i+nn2]; { Copy IM to 1st half IM position }
  For i := 2 to nn2 do
  Begin
    x[nx2-i+2] := -x[i+n];    { Replicate 2nd half Imag as negative }
    x[n-i+2] := x[i];         { Replicate 2nd half Real as mirror image }
  End;
  n := nx2;       { We have doubled number of points }
End;

{*****************************************************************************}

Procedure Post(inv : Xform);
  { Post processing for forward real transforms and pre-processing for inverse
    real transforms, depending on state of the variable inv. }

Var
  nn2,nn4,l,i,m,ipn2,mpn2   : Integer;
  arg,rmsin,rmcos,ipcos,ipsin,ic,is1,rp,rm,ip,im   : Real;

Begin
  nn2 := n div 2;   { n is global }
  nn4 := n div 4;
  { imax represents PI/2 }
  For l := 1 to nn4 Do
  { Start at ends of array and work towards middle }
  Begin
    i := l+1;       { Point near beginning }
    m := nn2-i+2;   { Point near end }
    ipn2 := i+nn2;  { Avoids recalculation each time }
    mpn2 := m+nn2;  { Calcs ptrs to imaginary part }
    rp := x[i]+x[m];
    rm := x[i]-x[m];
    ip := x[ipn2]+x[mpn2];
    im := x[ipn2]-x[mpn2];
    { Take cosine of PI/2 }
    arg := (Pi2/nn4)*(i-1);
    ic := Cos(arg);
    { Cosine term is minus if inverse }
    If inv = Inverse Then ic := -ic;
    Is1 := Sin(arg);
    ipcos := ip*ic;  { Avoid remultiplication below }
    ipsin := ip*is1;
    rmsin := rm*is1;
    rmcos := rm*ic;
    x[i] := (rp + ipcos - rmsin)/2.0;    {* Combine for real-1st pt }
    x[ipn2] := (im - ipsin - rmcos)/2.0; {* Imag-1st point }
    x[m] := (rp - ipcos + rmsin)/2.0;    {* Real - last pt }
    x[mpn2] := -(im +ipsin +rmcos)/2.0;  {* Imag, last pt }
  End;  { For }
  x[1] := x[1]+x[nn2+1];    {** For first complex pair}
  x[nn2+1] := 0.0;          {**}
End;    { Post }

{*****************************************************************************}

Procedure Shuffl(inv : Xform);
  { This procedure shuffels points from alternate real-imaginary to
    1st-half real, 2nd-half imaginary order if inv=Fwd, and reverses the
    process if inv=Inverse.  Algorithm is much like Cooley-Tukey:
      Starts with large cells and works to smaller ones for Fwd
      Starts with small cells and increases if inverse }

Var
  i,j,k,ipcm1,celdis,celnum,parnum  : Integer;
  xtemp                             : Real;

Begin
  { Choose whether to start with large cells and go down or start with small
    cells and increase }

  Case inv Of

  Fwd: Begin
    celdis := n div 2;     { Distance between cells }
    celnum := 1;           { One cell in first pass }
    parnum := n div 4;     { n/4 pairs per cell in 1st pass }
    End;   { Fwd case }

  Inverse: Begin
    celdis := 2;           { Cells are adjacent }
    celnum := n div 4;     { n/4 cells in first pass }
    parnum := 1;
    End;  { Inverse case }

  End;  { Case }

  Repeat      { Until cells large if Fwd or small if Inverse }
    i := 2;
    For j:= 1 to celnum do
    Begin
      For k := 1 to parnum do  { Do all pairs in each cell }
      Begin
        Xtemp := x[i];
        ipcm1 := i+celdis-1;   { Index of other pts }
        x[i] := x[ipcm1];
        x[ipcm1] := xtemp;
        i := i+2;
      End;  { For k }

      { End of cell, advance to next one }
      i := i+celdis;
    End;    { For j }

    { Change values for next pass }

    Case inv Of
    Fwd:    Begin
      celdis := celdis div 2;         { Decrease cell distance }
      celnum := celnum * 2;           { More cells }
      parnum := parnum div 2;         { Less pairs per cell }
      End;   { Case Fwd }

    Inverse: Begin
      celdis := celdis * 2;           { More distance between cells }
      celnum := celnum div 2;         { Less cells }
      parnum := parnum * 2;           { More pairs per cell }
      End;   { Case Inverse }
    End;  { Case }

  Until  (((inv = Fwd) And (Celdis < 2)) Or ((inv=Inverse) And (celnum = 0)));

End;  { Shuffl }

{*****************************************************************************}

Procedure FFT(inv : Xform);
  { Fast Fourier transform procedure operating on data in 1st half real,
    2nd half imaginary order and produces a complex result }

Var
  n1,n2,nu,celnum,celdis,parnum,ipn2,kpn2,jpn2,
  i,j,k,l,i2,imax,index      : Integer;
  arg,cosy,siny,r2cosy,r2siny,i2cosy,i2siny,picons,
  y,deltay,k1,k2,tr,ti,xtemp  : Real;

Begin
  { Calculate nu = log2(n) }
  nu := 0;
  n1 := n div 2;
  n2 := n1;
  While n1 >= 2 Do
  Begin
    nu := nu + 1;            { Increment power-of-2 counter }
    n1 := n1 div 2;          { divide by 2 until zero }
  End;
  { Shuffel the data into bit-inverted order }
  For i := 1 to n2 Do
  Begin
    k := ibitr(i-1,nu)+1;  { Calc bit-inverted position in array }
    If i>k Then            { Prevent swapping twice }
    Begin
      ipn2 := i+n2;
      kpn2 := k+n2;
      tr := x[k];         { Temp storage of real }
      ti := x[kpn2];      { Temp imag }
      x[k] := x[i];
      x[kpn2] := x[ipn2];
      x[i] := tr;
      x[ipn2] := ti;
    End;  { If }
  End;  { For }

  { Do first pass specially, since it has no multiplications }
  i := 1;
  While i <= n2 Do
  Begin
    k := i+1;
    kpn2 := k+n2;
    ipn2 := i+n2;
    k1 := x[i]+x[k];        { Save this sum }
    x[k] := x[i]-x[k];      { Diff to k's }
    x[i] := k1;             { Sum to I's }
    k1 := x[ipn2]+x[kpn2];  { Sum of imag }
    x[kpn2] := x[ipn2]-x[kpn2];
    x[ipn2] := k1;
    i := i+2;
  End;  { While }

  { Set up deltay for 2nd pass, deltay=PI/2 }
    deltay := PI2;   { PI/2 }
    celnum := n2 div 4;
    parnum := 2;     { Number of pairs between cell }
    celdis := 2;     { Distance between cells }


  { Each pass after 1st starts here }
  Repeat             { Until number of cells becomes zero }

{ Writeln(Lst,'After Nth Pass:');  ### }
{ Debug; }

    { Each new cell starts here }
    index := 1;
    y := 0;      { Exponent of w }
    { Do the number of pairs in each cell }
    For i2 := 1 To parnum Do
    Begin
      If y <> 0 Then
      Begin                { Use sines and cosines if y is not zero }
        cosy := cos(y);    { Calc sine and cosine }
        siny := sin(y);
        { Negate sine terms if transform is inverse }
        If inv = Inverse then siny := -siny;
      End;   { If }
      { These are the fundamental equations of the FFT }
      For l := 0 to celnum-1 Do
      Begin
        i := (celdis*2)*l + index;
        j := i+celdis;
        ipn2 := i + n2;
        jpn2 := j + n2;
        If y = 0 Then   { Special case for y=0 -- No sine or cosine terms }
        Begin
          k1 := x[i]+x[j];
          k2 := x[ipn2]+x[jpn2];
          x[j] := x[i]-x[j];
          x[jpn2] := x[ipn2]-x[jpn2];
        End   { If-Then }
        Else
        Begin
          r2cosy := x[j]*cosy;   { Calc intermediate constants }
          r2siny := x[j]*siny;
          i2cosy := x[jpn2]*cosy;
          i2siny := x[jpn2]*siny;
          { Butterfly }
          k1 := x[i] + r2cosy + i2siny;
          k2 := x[ipn2] - r2siny + i2cosy;
          x[j] := x[i] - r2cosy - i2siny;
          x[jpn2] := x[ipn2] + r2siny - i2cosy;
        End;   { Else }
        { Replace the i terms }
        x[i] := k1;
        x[ipn2] := k2;

      { Advance angle for next pair }
      End;  { For l }

      Y := y + deltay;
      index := index + 1;
    End;  { For i2 }

    { Pass done - change cell distance and number of cells }
    celnum := celnum div 2;
    parnum := parnum * 2;
    celdis := celdis * 2;
    deltay := deltay/2;

  Until celnum = 0;

End;  { FFT }

{*****************************************************************************}

Begin    { * Main program * }
  For i := 0 To Asize-1 Do
  Begin
    x[i] := 0.0;
  End;

{  Write('Enter number of points: ');
  Readln(n);}
  n := 32;
  If (n > Asize) Then
  Begin
    Writeln('Too large, will use maximum');
    n := Round(asize/2.0);
  End;
  For i := 2 to n Do x[i] := Exp((1-i)*0.25); { Create Real array }
  x[1] := 0.5;
  Writeln('Input Array:');
  Debug;
  Shuffl(Fwd);
  FFT(Fwd);
  Post(Fwd);
  For i:= 1 to n Do x[i] := x[i]*8/n;
  Writeln('Forward FFT with real array first:');
  Debug;
End.