fftspec.pas

Delphi FFT Spectrum Analyzer

{ **************************************************************************** }
{ FileName............: FFTSpec.PAS                                            }
{ Project.............: FFT                                                    }
{ Author(s)...........: M.Majoor                                               }
{ Version.............: 2.00                                                   }
{ Last revision date..: 23 January, 2001                                       }
{ ---------------------------------------------------------------------------- }
{ Calculate FFT spectrum unit.                                                 }
{                                                                              }
{ Version  Date      Comment                                                   }
{ 1.00     19980502  - Initial release                                         }
{ 2.00     20010123  - Recompiled for Delphi 3 and made it up-to-date          }
{                    - Less fixed settings (more dynamic)                      }
{ **************************************************************************** }
unit FFTSpec;

interface
type
  TFFTArray  = array[0..$FFFF] of single;                                      { Placeholder }
  PFFTArray  = ^TFFTArray;
  TFFTArray2 = array[0..$FFFF] of smallint;                                    { Placeholder }
  PFFTArray2 = ^TFFTArray2;

  SpectrumWindows = (idRectangle, idBlackman, idCos2, idGaussian, idHamming,
                     idKaiser, idTriangle);

var
  BlackmanAlpha : single;
  BlackmanBeta  : single;
  BlackmanGamma : single;
  GaussianAlpha : single;
  HammingAlpha  : single;
  KaiserAlpha   : single;

procedure Spectrum1(WindowFunction: SpectrumWindows;
                    DataSize      : word;
                    RealIn        : PFFTArray;
                    RealOut       : PFFTArray);
procedure Spectrum2(WindowFunction: SpectrumWindows;
                    DataSize      : word;
                    RealIn        : PFFTArray2;
                    RealOut       : PFFTArray);

procedure CalcWindowFunctions(DataSize: word);

var
  CoherentGain    : array[idRectangle..idTriangle] of single;

implementation
uses SysUtils;

var
   WindowingTables : array[idRectangle..idTriangle] of array[0..$FFFF] of single;
   ReversalBits    : array[0..$FFFF] of word;
   AlphaTable      : array[0..$FFFF] of single;
   BetaTable       : array[0..$FFFF] of single;


{ **************************************************************************** }
{ Params   : <X>       Value to calculate Bessel for                           }
{ Return   : <Result>  Bessel value                                            }
{                                                                              }
{ Descript : Return Bessel result for index X.                                 }
{ Notes    : Bessel is calculated using:                                       }
{            (x^0/1)+(x^2/4)+(x^4/64)+(x^6/2304)+(x^8/14745600).... -> 25 terms}
{            The 'top' numbers increase with 2 for each term or alternatively  }
{            termindex*2. Also it can be seen that the upperterm is everytime  }
{            sqr(x) larger for each term.                                      }
{            This in a way also applies to the bottom numbers. These are       }
{            calculated with:                                                  }
{              sqr( (termindex*2)*(termindex-1)*2) )                           }
{ **************************************************************************** }
function CalcBessel(X: single): single;
var
  PrevBTerm : single;                  { Previous bottom term = (termindex-1)*2 }
  BTerm     : single;                  { Bottom term }
  Terms     : integer;                 { Term loop }
  UTermMult : single;                  { Upper term multiplier }
  UTerm     : single;                  { Upper term }
  CalcResult: single;                  { Result of Bessel function }
begin
  UTermMult  := sqr(X);                { Multiplier upperterm value }
  UTerm      := UTermMult;
  PrevBTerm  := 1;
  CalcResult := 1.0;
  for Terms := 1 to 10 do              { Should be accurate enough }
  begin
    BTerm := (Terms * 2) * PrevBTerm;
    PrevBTerm := BTerm;
    BTerm := sqr(BTerm);
    CalcResult := CalcResult + (UTerm / BTerm);
    UTerm := UTerm * UTermMult;        { Upper term increments }
  end;
  Result := CalcResult;
end;


{ **************************************************************************** }
{ Params   : <Data>    Data to reverse                                         }
{            <NumBits> Number of bits to reverse                               }
{ Return   : <Result>  Reversed bit data                                       }
{                                                                              }
{ Descript : Reverse the bits from indicated bit index.                        }
{ Notes    :                                                                   }
{ **************************************************************************** }
function ReverseBits(Data, NumBits: word): word;
var
  Loop    : word;
  Reversed: word;
begin
  Reversed := 0;
  for Loop := 0 to NumBits-1 do
  begin
    Reversed := (Reversed shl 1) or (Data and $01);
    Data := Data shr 1;
  end;
  Result := Reversed;
end;


{ **************************************************************************** }
{ Params   : <WindowsFunction>  Windowing to perform                           }
{            <DataSize>         Number of data items                           }
{                               MUST be 2^x number (0, 2, 4, 8, 16, 32..)      }                    
{            <RealIn>           Data to process                                }
{ Return   : <RealOut>          Processed data                                 }
{                                                                              }
{ Descript : Calculate FFT followed by spectrum (using reals as input).        }
{ Notes    :                                                                   }
{ **************************************************************************** }
procedure Spectrum1(WindowFunction: SpectrumWindows;
                    DataSize      : word;
                    RealIn        : PFFTArray;
                    RealOut       : PFFTArray);
var
  i        : integer;
  j, k, n  : word;
  BlockSize: word;
  BlockEnd : word;
  Alpha    : single;
  Beta     : single;
  Delta_ar : single;
  tr, ti   : single;
  ar, ai   : single;
  ImagOut  : PFFTArray;                { Temporary buffer for imaginary data }
begin
  ImagOut := AllocMem(DataSize * sizeof(single));                              { Allocate and clear }
  try
    for i := 0 to DataSize-1
      do RealOut^[ReversalBits[i]] := RealIn^[i]*WindowingTables[WindowFunction][i];

    BlockEnd  := 1;
    BlockSize := 2;
    while BlockSize <= DataSize do
    begin
      Alpha := AlphaTable[BlockSize-2];
      Beta  := BetaTable [BlockSize-2];
      i := 0;
      while i < DataSize do
      begin
        ar := 1.0;                                           { cos(0) }
        ai := 0.0;                                           { sin(0) }
        j := i;
        for n := 0 to BlockEnd-1 do
        begin
          k  := j + BlockEnd;                                { K  = 0..DataSize-1 }
          tr := ar*RealOut^[k] - ai*ImagOut^[k];             { tr = appr. 0.5 * datainput in realout }
          ti := ar*ImagOut^[k] + ai*RealOut^[k];             { ti = appr. 0.5 * datainput in realout }
          RealOut^[k] := RealOut^[j] - tr;
          ImagOut^[k] := ImagOut^[j] - ti;
          RealOut^[j] := RealOut^[j] + tr;
          ImagOut^[j] := ImagOut^[j] + ti;
          delta_ar := Alpha*ar + Beta*ai;                    { delta_ar = 0..2 }
          ai := ai - (Alpha*ai - Beta*ar);                   { ai = 0..1 }
          ar := ar - Delta_ar;                               { ar = -1..1 }
          inc (j);
        end;
        inc(i, BlockSize);
      end;
      BlockEnd  := BlockSize;
      BlockSize := BlockSize shl 1;
    end;

    { At this point FFT is calculated, now calculate spectrum }
    for i := 0 to DataSize-1 do RealOut^[i] := sqrt (sqr(RealOut^[i]) + sqr(ImagOut^[i])) *
                                               CoherentGain[WindowFunction];
  finally
    FreeMem(ImagOut);
  end;
end;


{ **************************************************************************** }
{ Params   : <WindowsFunction>  Windowing to perform                           }
{            <DataSize>         Number of data items                           }
{                               MUST be 2^x number (0, 2, 4, 8, 16, 32..)      }