aboutsummaryrefslogtreecommitdiffstats
path: root/Game/Code/lib/fft/UFFT.pas
diff options
context:
space:
mode:
authork-m_schindler <k-m_schindler@b956fd51-792f-4845-bead-9b4dfca2ff2c>2008-08-27 13:28:57 +0000
committerk-m_schindler <k-m_schindler@b956fd51-792f-4845-bead-9b4dfca2ff2c>2008-08-27 13:28:57 +0000
commit1ba91d5a0e1df7419a561f6dcf16a0839509a5e7 (patch)
tree3f76e96fc5a3f5b738dabce28642ff2415748ccb /Game/Code/lib/fft/UFFT.pas
parente9fd8ce40b4cbf006695fd6e56f84071407843c9 (diff)
downloadusdx-1ba91d5a0e1df7419a561f6dcf16a0839509a5e7.tar.gz
usdx-1ba91d5a0e1df7419a561f6dcf16a0839509a5e7.tar.xz
usdx-1ba91d5a0e1df7419a561f6dcf16a0839509a5e7.zip
Reordering of the directories[1]: moving Game/Code to src
git-svn-id: svn://svn.code.sf.net/p/ultrastardx/svn/trunk@1302 b956fd51-792f-4845-bead-9b4dfca2ff2c
Diffstat (limited to 'Game/Code/lib/fft/UFFT.pas')
-rw-r--r--Game/Code/lib/fft/UFFT.pas600
1 files changed, 0 insertions, 600 deletions
diff --git a/Game/Code/lib/fft/UFFT.pas b/Game/Code/lib/fft/UFFT.pas
deleted file mode 100644
index 87c981e0..00000000
--- a/Game/Code/lib/fft/UFFT.pas
+++ /dev/null
@@ -1,600 +0,0 @@
-{**********************************************************************
-
- 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+} // Use AnsiString
-{$ENDIF}
-
-interface
-type
- TSingleArray = array[0..0] of Single;
- PSingleArray = ^TSingleArray;
-
- TFFTWindowFunc = (
- fwfRectangular,
- fwfBartlett,
- fwfHamming,
- fwfHanning,
- fwfBlackman,
- fwfBlackman_Harris,
- fwfWelch,
- fwfGaussian2_5,
- fwfGaussian3_5,
- fwfGaussian4_5
- );
-
-const
- FFTWindowName: array[TFFTWindowFunc] of string = (
- 'Rectangular',
- 'Bartlett',
- 'Hamming',
- 'Hanning',
- 'Blackman',
- 'Blackman-Harris',
- 'Welch',
- 'Gaussian(a=2.5)',
- 'Gaussian(a=3.5)',
- 'Gaussian(a=4.5)'
- );
-
-(*
- * This is the function you will use the most often.
- * Given an array of floats, this will compute the power
- * spectrum by doing a Real FFT and then computing the
- * sum of the squares of the real and imaginary parts.
- * Note that the output array is half the length of the
- * input array, and that NumSamples must be a power of two.
- *)
-procedure PowerSpectrum(NumSamples: Integer; In_, Out_: PSingleArray);
-
-(*
- * Computes an FFT when the input data is real but you still
- * want complex data as output. The output arrays are half
- * the length of the input, and NumSamples must be a power of
- * two.
- *)
-procedure RealFFT(NumSamples: integer;
- RealIn, RealOut, ImagOut: PSingleArray);
-
-(*
- * Computes a FFT of complex input and returns complex output.
- * Currently this is the only function here that supports the
- * inverse transform as well.
- *)
-procedure FFT(NumSamples: Integer;
- InverseTransform: boolean;
- RealIn, ImagIn, RealOut, ImagOut: PSingleArray);
-
-(*
- * Applies a windowing function to the data in place
- *
- * 0: Rectangular (no window)
- * 1: Bartlett (triangular)
- * 2: Hamming
- * 3: Hanning
- * 4: Blackman
- * 5: Blackman-Harris
- * 6: Welch
- * 7: Gaussian(a=2.5)
- * 8: Gaussian(a=3.5)
- * 9: Gaussian(a=4.5)
- *)
-procedure WindowFunc(whichFunction: TFFTWindowFunc; NumSamples: Integer; in_: PSingleArray);
-
-(*
- * Returns the name of the windowing function (for UI display)
- *)
-function WindowFuncName(whichFunction: TFFTWindowFunc): string;
-
-(*
- * Returns the number of windowing functions supported
- *)
-function NumWindowFuncs(): integer;
-
-
-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; {$IFDEF HasInline}inline;{$ENDIF}
-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;
-
-(*
- * Real Fast Fourier Transform
- *
- * This function was based on the code in Numerical Recipes in C.
- * In Num. Rec., the inner loop is based on a single 1-based array
- * of interleaved real and imaginary numbers. Because we have two
- * separate zero-based arrays, our indices are quite different.
- * Here is the correspondence between Num. Rec. indices and our indices:
- *
- * i1 <-> real[i]
- * i2 <-> imag[i]
- * i3 <-> real[n/2-i]
- * i4 <-> imag[n/2-i]
- *)
-procedure RealFFT(NumSamples: integer; RealIn, RealOut, ImagOut: PSingleArray);
-var
- Half: Integer;
- i: Integer;
- theta: Single;
- tmpReal, tmpImag: PSingleArray;
- wtemp: Single;
- wpr, wpi, wr, wi: Single;
- i3: Integer;
- h1r, h1i, h2r, h2i: Single;
-begin
- Half := NumSamples div 2;
-
- theta := Pi / Half;
-
- GetMem(tmpReal, Half * sizeof(Single));
- GetMem(tmpImag, Half * sizeof(Single));
-
- for i := 0 to Half-1 do
- begin
- tmpReal[i] := RealIn[2 * i];
- tmpImag[i] := RealIn[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]);
-
- RealOut[i] := h1r + wr * h2r - wi * h2i;
- ImagOut[i] := h1i + wr * h2i + wi * h2r;
- RealOut[i3] := h1r - wr * h2r + wi * h2i;
- ImagOut[i3] := -h1i + wr * h2i + wi * h2r;
-
- wtemp := wr;
- wr := wtemp * wpr - wi * wpi + wr;
- wi := wi * wpr + wtemp * wpi + wi;
- end;
-
- h1r := RealOut[0];
- RealOut[0] := h1r + ImagOut[0];
- ImagOut[0] := h1r - ImagOut[0];
-
- FreeMem(tmpReal);
- FreeMem(tmpImag);
-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;
-
-(*
- * Windowing Functions
- *)
-function NumWindowFuncs(): integer;
-begin
- Result := Length(FFTWindowName);
-end;
-
-function WindowFuncName(whichFunction: TFFTWindowFunc): string;
-begin
- Result := FFTWindowName[whichFunction];
-end;
-
-procedure WindowFunc(whichFunction: TFFTWindowFunc; NumSamples: Integer; in_: PSingleArray);
-var
- i: Integer;
- A: Single;
-begin
- case whichFunction of
- fwfBartlett:
- begin
- // Bartlett (triangular) window
- for i := 0 to (NumSamples div 2)-1 do
- begin
- in_[i] := in_[i] * (i / (NumSamples / 2));
- in_[i + (NumSamples div 2)] :=
- in_[i + (NumSamples div 2)] *
- (1.0 - (i / (NumSamples / 2)));
- end;
- end;
- fwfHamming:
- begin
- // Hamming
- for i := 0 to NumSamples-1 do
- begin
- in_[i] := in_[i] * (0.54 - 0.46 * cos(2 * Pi * i / (NumSamples - 1)));
- end;
- end;
- fwfHanning:
- begin
- // Hanning
- for i := 0 to NumSamples-1 do
- begin
- in_[i] := in_[i] * (0.50 - 0.50 * cos(2 * Pi * i / (NumSamples - 1)));
- end;
- end;
- fwfBlackman:
- begin
- // Blackman
- for i := 0 to NumSamples-1 do
- begin
- in_[i] := in_[i] * (0.42 - 0.5 * cos (2 * Pi * i / (NumSamples - 1)) + 0.08 * cos (4 * Pi * i / (NumSamples - 1)));
- end;
- end;
- fwfBlackman_Harris:
- begin
- // Blackman-Harris
- for i := 0 to NumSamples-1 do
- begin
- in_[i] := in_[i] * (0.35875 - 0.48829 * cos(2 * Pi * i /(NumSamples-1)) + 0.14128 * cos(4 * Pi * i/(NumSamples-1)) - 0.01168 * cos(6 * Pi * i/(NumSamples-1)));
- end;
- end;
- fwfWelch:
- begin
- // Welch
- for i := 0 to NumSamples-1 do
- begin
- in_[i] := in_[i] * 4*i/NumSamples*(1-(i/NumSamples));
- end;
- end;
- fwfGaussian2_5:
- begin
- // Gaussian (a=2.5)
- // Precalculate some values, and simplify the fmla to try and reduce overhead
- A := -2*2.5*2.5;
-
- for i := 0 to NumSamples-1 do
- begin
- // full
- // in_[i] := in_[i] * exp(-0.5*(A*((i-NumSamples/2)/NumSamples/2))*(A*((i-NumSamples/2)/NumSamples/2)));
- // reduced
- //in_[i] := in_[i] * exp(A*(0.25 + ((i/NumSamples)*(i/NumSamples)) - (i/NumSamples)));
- end;
- end;
- fwfGaussian3_5:
- begin
- // Gaussian (a=3.5)
- A := -2*3.5*3.5;
-
- for i := 0 to NumSamples-1 do
- begin
- // reduced
- in_[i] := in_[i] * exp(A*(0.25 + ((i/NumSamples)*(i/NumSamples)) - (i/NumSamples)));
- end;
- end;
- fwfGaussian4_5:
- begin
- // Gaussian (a=4.5)
- A := -2*4.5*4.5;
-
- for i := 0 to NumSamples-1 do
- begin
- // reduced
- in_[i] := in_[i] * exp(A*(0.25 + ((i/NumSamples)*(i/NumSamples)) - (i/NumSamples)));
- end;
- end;
- end;
-end;
-
-end.