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.