aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Game/Code/lib/fft/UFFT.pas345
1 files changed, 345 insertions, 0 deletions
diff --git a/Game/Code/lib/fft/UFFT.pas b/Game/Code/lib/fft/UFFT.pas
new file mode 100644
index 00000000..e97812f6
--- /dev/null
+++ b/Game/Code/lib/fft/UFFT.pas
@@ -0,0 +1,345 @@
+{**********************************************************************
+
+ FFT.cpp
+
+ Dominic Mazzoni
+
+ September 2000
+
+***********************************************************************
+
+Fast Fourier Transform routines.
+
+ This file contains a few FFT routines, including a real-FFT
+ routine that is almost twice as fast as a normal complex FFT,
+ and a power spectrum routine when you know you don't care
+ about phase information.
+
+ Some of this code was based on a free implementation of an FFT
+ by Don Cross, available on the web at:
+
+ http://www.intersrv.com/~dcross/fft.html
+
+ The basic algorithm for his code was based on Numerican Recipes
+ in Fortran. I optimized his code further by reducing array
+ accesses, caching the bit reversal table, and eliminating
+ float-to-double conversions, and I added the routines to
+ calculate a real FFT and a real power spectrum.
+
+***********************************************************************
+
+ Salvo Ventura - November 2006
+ Added more window functions:
+ * 4: Blackman
+ * 5: Blackman-Harris
+ * 6: Welch
+ * 7: Gaussian(a=2.5)
+ * 8: Gaussian(a=3.5)
+ * 9: Gaussian(a=4.5)
+
+***********************************************************************
+
+ This file is part of Audacity 1.3.4 beta (http://audacity.sourceforge.net/)
+ Ported to Pascal by the UltraStar Deluxe Team
+}
+
+unit UFFT;
+
+{$ifdef fpc}
+ {$mode delphi}{$H+}
+{$endif}
+
+interface
+type
+ TSingleArray = array[0..0] of Single;
+ PSingleArray = ^TSingleArray;
+
+procedure PowerSpectrum(NumSamples: Integer; In_, Out_: PSingleArray);
+procedure WindowFunc(NumSamples: Integer; in_: PSingleArray); inline;
+
+implementation
+
+uses
+ SysUtils;
+
+type TIntArray = array[0..0] of Integer;
+type TIntIntArray = array[0..0] of ^TIntArray;
+var gFFTBitTable: ^TIntIntArray;
+const MaxFastBits: Integer = 16;
+
+function IsPowerOfTwo(x: Integer): Boolean;
+begin
+ if (x < 2) then
+ result := false
+ else if ((x and (x - 1)) <> 0) then { Thanks to 'byang' for this cute trick! }
+ result := false
+ else
+ result := true;
+end;
+
+function NumberOfBitsNeeded(PowerOfTwo: Integer): Integer;
+var i: Integer;
+begin
+ if (PowerOfTwo < 2) then begin
+ Writeln(ErrOutput, Format('Error: FFT called with size %d\n', [PowerOfTwo]));
+ Abort;
+ end;
+
+ i := 0;
+ while(true) do begin
+ if (PowerOfTwo and (1 shl i) <> 0) then begin
+ result := i;
+ Exit;
+ end;
+ Inc(i);
+ end;
+end;
+
+function ReverseBits(index, NumBits: Integer): Integer;
+var
+ i, rev: Integer;
+begin
+ rev := 0;
+ for i := 0 to NumBits-1 do begin
+ rev := (rev shl 1) or (index and 1);
+ index := index shr 1;
+ end;
+
+ result := rev;
+end;
+
+procedure InitFFT();
+var
+ len: Integer;
+ b, i: Integer;
+begin
+ GetMem(gFFTBitTable, MaxFastBits * sizeof(PSingle));
+
+ len := 2;
+ for b := 1 to MaxFastBits do begin
+ GetMem(gFFTBitTable^[b - 1], len * sizeof(Single));
+ for i := 0 to len-1 do
+ gFFTBitTable^[b - 1][i] := ReverseBits(i, b);
+ len := len shl 1;
+ end;
+end;
+
+function FastReverseBits(i, NumBits: Integer): Integer; inline;
+begin
+ if (NumBits <= MaxFastBits) then
+ result := gFFTBitTable[NumBits - 1][i]
+ else
+ result := ReverseBits(i, NumBits);
+end;
+
+{*
+ * Complex Fast Fourier Transform
+ *}
+procedure FFT(NumSamples: Integer;
+ InverseTransform: boolean;
+ RealIn, ImagIn, RealOut, ImagOut: PSingleArray);
+var
+ NumBits: Integer; { Number of bits needed to store indices }
+ i, j, k, n: Integer;
+ BlockSize, BlockEnd: Integer;
+ delta_angle: Double;
+ angle_numerator: Double;
+ tr, ti: Double; { temp real, temp imaginary }
+ sm2, sm1, cm2, cm1: Double;
+ w: Double;
+ ar0, ar1, ar2, ai0, ai1, ai2: Double;
+ denom: Single;
+begin
+
+ angle_numerator := 2.0 * Pi;
+
+ if (not IsPowerOfTwo(NumSamples)) then begin
+ Writeln(ErrOutput, Format('%d is not a power of two', [NumSamples]));
+ Abort;
+ end;
+
+ if (gFFTBitTable = nil) then
+ InitFFT();
+
+ if (InverseTransform) then
+ angle_numerator := -angle_numerator;
+
+ NumBits := NumberOfBitsNeeded(NumSamples);
+
+ {
+ ** Do simultaneous data copy and bit-reversal ordering into outputs...
+ }
+
+ for i := 0 to NumSamples-1 do begin
+ j := FastReverseBits(i, NumBits);
+ RealOut[j] := RealIn[i];
+ if(ImagIn = nil) then
+ ImagOut[j] := 0.0
+ else
+ ImagOut[j] := ImagIn[i];
+ end;
+
+ {
+ ** Do the FFT itself...
+ }
+
+ BlockEnd := 1;
+ BlockSize := 2;
+ while(BlockSize <= NumSamples) do
+ begin
+
+ delta_angle := angle_numerator / BlockSize;
+
+ sm2 := sin(-2 * delta_angle);
+ sm1 := sin(-delta_angle);
+ cm2 := cos(-2 * delta_angle);
+ cm1 := cos(-delta_angle);
+ w := 2 * cm1;
+
+ i := 0;
+ while(i < NumSamples) do
+ begin
+ ar2 := cm2;
+ ar1 := cm1;
+
+ ai2 := sm2;
+ ai1 := sm1;
+
+ j := i;
+ for n := 0 to BlockEnd-1 do
+ begin
+ ar0 := w * ar1 - ar2;
+ ar2 := ar1;
+ ar1 := ar0;
+
+ ai0 := w * ai1 - ai2;
+ ai2 := ai1;
+ ai1 := ai0;
+
+ k := j + BlockEnd;
+ tr := ar0 * RealOut[k] - ai0 * ImagOut[k];
+ ti := ar0 * ImagOut[k] + ai0 * RealOut[k];
+
+ RealOut[k] := RealOut[j] - tr;
+ ImagOut[k] := ImagOut[j] - ti;
+
+ RealOut[j] := RealOut[j] + tr;
+ ImagOut[j] := ImagOut[j] + ti;
+
+ Inc(j);
+ end;
+
+ Inc(i, BlockSize);
+ end;
+
+ BlockEnd := BlockSize;
+ BlockSize := BlockSize shl 1;
+ end;
+
+ {
+ ** Need to normalize if inverse transform...
+ }
+
+ if (InverseTransform) then begin
+ denom := NumSamples;
+
+ for i := 0 to NumSamples-1 do begin
+ RealOut[i] := RealOut[i] / denom;
+ ImagOut[i] := ImagOut[i] / denom;
+ end;
+ end;
+end;
+
+{*
+ * PowerSpectrum
+ *
+ * This function computes the same as RealFFT, above, but
+ * adds the squares of the real and imaginary part of each
+ * coefficient, extracting the power and throwing away the
+ * phase.
+ *
+ * For speed, it does not call RealFFT, but duplicates some
+ * of its code.
+ *}
+procedure PowerSpectrum(NumSamples: Integer; In_, Out_: PSingleArray);
+var
+ Half: Integer;
+ i: Integer;
+ theta: Single;
+ tmpReal, tmpImag, RealOut, ImagOut: PSingleArray;
+ wtemp: Single;
+ wpr, wpi, wr, wi: Single;
+ i3: Integer;
+ h1r, h1i, h2r, h2i, rt, it: Single;
+begin
+ Half := NumSamples div 2;
+
+ theta := Pi / Half;
+
+ GetMem(tmpReal, Half * sizeof(Single));
+ GetMem(tmpImag, Half * sizeof(Single));
+ GetMem(RealOut, Half * sizeof(Single));
+ GetMem(ImagOut, Half * sizeof(Single));
+
+ for i := 0 to Half-1 do begin
+ tmpReal[i] := In_[2 * i];
+ tmpImag[i] := In_[2 * i + 1];
+ end;
+
+ FFT(Half, false, tmpReal, tmpImag, RealOut, ImagOut);
+
+ wtemp := sin(0.5 * theta);
+
+ wpr := -2.0 * wtemp * wtemp;
+ wpi := sin(theta);
+ wr := 1.0 + wpr;
+ wi := wpi;
+
+ for i := 1 to (Half div 2)-1 do
+ begin
+ i3 := Half - i;
+
+ h1r := 0.5 * (RealOut[i] + RealOut[i3]);
+ h1i := 0.5 * (ImagOut[i] - ImagOut[i3]);
+ h2r := 0.5 * (ImagOut[i] + ImagOut[i3]);
+ h2i := -0.5 * (RealOut[i] - RealOut[i3]);
+
+ rt := h1r + wr * h2r - wi * h2i;
+ it := h1i + wr * h2i + wi * h2r;
+
+ Out_[i] := rt * rt + it * it;
+
+ rt := h1r - wr * h2r + wi * h2i;
+ it := -h1i + wr * h2i + wi * h2r;
+
+ Out_[i3] := rt * rt + it * it;
+
+ wtemp := wr;
+ wr := wtemp * wpr - wi * wpi + wr;
+ wi := wi * wpr + wtemp * wpi + wi;
+ end;
+
+ h1r := RealOut[0];
+ rt := h1r + ImagOut[0];
+ it := h1r - ImagOut[0];
+ Out_[0] := rt * rt + it * it;
+
+ rt := RealOut[Half div 2];
+ it := ImagOut[Half div 2];
+ Out_[Half div 2] := rt * rt + it * it;
+
+ FreeMem(tmpReal);
+ FreeMem(tmpImag);
+ FreeMem(RealOut);
+ FreeMem(ImagOut);
+end;
+
+procedure WindowFunc(NumSamples: Integer; in_: PSingleArray); inline;
+var
+ i: Integer;
+begin
+ for i := 0 to NumSamples-1 do
+ in_[i] := in_[i] * (0.50 - 0.50 * cos(2 * Pi * i / (NumSamples - 1)));
+end;
+
+end.
generated by cgit v1.2.3 (git 2.25.1) at 2024-12-16 07:13:58 +0000