diff options
author | tobigun <tobigun@b956fd51-792f-4845-bead-9b4dfca2ff2c> | 2010-10-10 22:59:33 +0000 |
---|---|---|
committer | tobigun <tobigun@b956fd51-792f-4845-bead-9b4dfca2ff2c> | 2010-10-10 22:59:33 +0000 |
commit | 35187604cef84864a908972d07361a5bd57e29ca (patch) | |
tree | dc95a8b1abeabd3a466729056ab8d37aaa6e72ea /src/lib/fft | |
parent | 58c1daf3692d4c5c534750a4fda97e087b0f0cbb (diff) | |
parent | 02bd10f0798829ab69d2028b988cb2a54eae292a (diff) | |
download | usdx-35187604cef84864a908972d07361a5bd57e29ca.tar.gz usdx-35187604cef84864a908972d07361a5bd57e29ca.tar.xz usdx-35187604cef84864a908972d07361a5bd57e29ca.zip |
rename trunk to 1.1svn/1.1github/svn/1.1
git-svn-id: svn://svn.code.sf.net/p/ultrastardx/svn/branches/1.1@2662 b956fd51-792f-4845-bead-9b4dfca2ff2c
Diffstat (limited to '')
-rw-r--r-- | src/lib/fft/UFFT.pas | 602 |
1 files changed, 602 insertions, 0 deletions
diff --git a/src/lib/fft/UFFT.pas b/src/lib/fft/UFFT.pas new file mode 100644 index 00000000..5a056a8c --- /dev/null +++ b/src/lib/fft/UFFT.pas @@ -0,0 +1,602 @@ +{********************************************************************** + + 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 long strings +{$ENDIF} + +interface +type + TSingleArray = array[0 .. (MaxInt div SizeOf(Single))-1] 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 .. (MaxInt div SizeOf(Integer))-1] of Integer; +type PIntArray = ^TIntArray; +type TIntIntArray = array[0 .. (MaxInt div SizeOf(PIntArray))-1] of PIntArray; +type PIntIntArray = ^TIntIntArray; +var gFFTBitTable: PIntIntArray; +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. |