{ > 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.