mmfft.pas

一套及时通讯的原码

{========================================================================}
{=                (c) 1995-98 SwiftSoft Ronald Dittrich                 =}
{========================================================================}
{=                          All Rights Reserved                         =}
{========================================================================}
{=  D 01099 Dresden             = Fax.: [phone]              =}
{=  Loewenstr.7a                = [e-mail]                     =}
{========================================================================}
{=  Actual versions on http://www.swiftsoft.de/mmtools.html             =}
{========================================================================}
{=  This code is for reference purposes only and may not be copied or   =}
{=  distributed in any format electronic or otherwise except one copy   =}
{=  for backup purposes.                                                =}
{=                                                                      =}
{=  No Delphi Component Kit or Component individually or in a collection=}
{=  subclassed or otherwise from the code in this unit, or associated   =}
{=  .pas, .dfm, .dcu, .asm or .obj files may be sold or distributed     =}
{=  without express permission from SwiftSoft.                          =}
{=                                                                      =}
{=  For more licence informations please refer to the associated        =}
{=  HelpFile.                                                           =}
{========================================================================}
{=  $Date: 13.04.98 - 20:21:43 $                                        =}
{========================================================================}
unit MMFFT;

{$I COMPILER.INC}

interface

uses
{$IFDEF WIN32}
    Windows,
{$ELSE}
    WinTypes,
    WinProcs,
{$ENDIF}
    SysUtils,
    MMUtils,
    MMMath,
    MMAbout;

{$IFDEF WIN32}
type
   PfCplx = ^TfCplx;
   TfCplx = record
      re: FLoat;
      im: Float;
   end;

   PfCplxArray = ^TfCplxArray;
   TfCplxArray = array[0..0] of TfCplx;

   PFFTCplx    = ^TFFTCplx;
   TFFTCplx    = record
      Order      : integer;
      FFTLen     : integer;
      Scale      : Float;
      BitRevTable: PIntArray;
      SinTable   : PfCplxArray;
   end;

   PFFTReal    = ^TFFTReal;
   TFFTReal    = record
      Order      : integer;
      FFTLen     : integer;
      Scale      : Float;
      BitRevTable: PIntArray;
      SinTable1  : PfCplxArray;
      SinTable2  : PfCplxArray;
   end;

var
   fvecMul1       : procedure(Vec: PFloatArray; val: Float; n: integer); {$IFDEF WIN32}pascal;{$ENDIF}
   fvecMul2       : procedure(SrcVec,DstVec: PfCplxArray; n: integer); {$IFDEF WIN32}pascal;{$ENDIF}
   fvecZero1      : procedure(Vec: PFloatArray; n: integer); {$IFDEF WIN32}pascal;{$ENDIF}
   fvecZero2      : procedure(Vec: PfCplxArray; n: integer); {$IFDEF WIN32}pascal;{$ENDIF}
   fvecBitRev1    : procedure(Vec: PFloatArray; BitRevTable: PIntArray; Order: integer); {$IFDEF WIN32}pascal;{$ENDIF}
   fvecBitRev2    : procedure(Vec: PfCplxArray; BitRevTable: PIntArray; Order: integer); {$IFDEF WIN32}pascal;{$ENDIF}
   GenSinTable1   : procedure(SinTable: PfCplxArray; Order: integer); {$IFDEF WIN32}pascal;{$ENDIF}
   GenSinTable2   : procedure(SinTable: PfCplxArray; Order: integer); {$IFDEF WIN32}pascal;{$ENDIF}
   DoCplxFFT1     : procedure(re, im: PFloatArray; scTable: PfCplxArray; FFTLen, sign: integer); {$IFDEF WIN32}pascal;{$ENDIF}
   DoCplxFFT2     : procedure(samps, scTable: PfCplxArray; FFTLen, sign: integer); {$IFDEF WIN32}pascal;{$ENDIF}
   BuildRealOutput: procedure(samps, SinTable: PfCplxArray; N,Sign: integer); {$IFDEF WIN32}pascal;{$ENDIF}

function  InitCplxFFT(iOrder: integer): PFFTCplx;
procedure DoneCplxFFT(var pfft: PFFTCplx);
procedure DoCplxFFTa(pfft: PFFTCplx; re, im: PFLoatArray; Sign: integer);
procedure DoCplxFFTb(pfft: PFFTCplx; samps: PfCplxArray; Sign: integer);

function  InitRealFFT(iOrder: integer): PFFTReal;
procedure DoneRealFFT(var pfft: PFFTReal);
procedure DoRealFFT(pfft: PFFTReal; Samps: PFloatArray; Sign: integer);
{$ENDIF}

function  CalcWindowFunc(Window, idx, Len: integer): Float;
procedure GenWindowTableInt(pTable: PIntArray; Window, Order: integer);
procedure GenWindowTableFloat(pTable: PFloatArray; Window, Order: integer);
procedure CalcWindowInt(Samps: PIntArray; Window,Order: integer);
procedure CalcWindowFloat(Samps: PFloatArray; Window,Order: integer);
procedure GenBitRevTable(pTable: PIntArray; Order: integer);

{$IFNDEF WIN32}

type
    {-- TMMFFT -----------------------------------------------------------}
    TMMFFT = class(TObject)
    private
      FPoints     : integer;
      FSinTable   : PSmallArray;
      FBitReversed: PIntArray;

      procedure SetFFTLen(FFTLen: integer);
      procedure InitFFT(FFTLen: integer);
      procedure DoneFFT;

    protected
      constructor Create; virtual;
      destructor  Destroy; override;
    public
      procedure CalcFFT(Buffer: PSmallInt);

      property  FFTLength: integer read FPoints write SetFFTLen default 128;
      property  BitReversed: PIntArray read FBitReversed;
      property  SinTable: PSmallArray read FSinTable;
    end;
{$ENDIF}

implementation

{$IFDEF WIN32}

{$IFNDEF USEDLL}
{$L MMFFT32.OBJ}
{$F+}
procedure _fvecMul1(Vec: PFloatArray; val: Float; n: integer); pascal; external;
procedure _fvecMul2(SrcVec,DstVec: PfCplxArray; n: integer); pascal; external;

procedure _fvecZero1(Vec: PFloatArray; n: integer); pascal; external;
procedure _fvecZero2(Vec: PfCplxArray; n: integer); pascal; external;

procedure _fvecBitRev1(Vec: PFloatArray; BitRevTable: PIntArray; Order: integer); pascal; external;
procedure _fvecBitRev2(Vec: PfCplxArray; BitRevTable: PIntArray; Order: integer); pascal; external;

procedure _GenSinTable1(SinTable: PfCplxArray; Order: integer); pascal; external;
procedure _GenSinTable2(SinTable: PfCplxArray; Order: integer); pascal; external;

procedure _DoCplxFFT1(re, im: PFloatArray; scTable: PfCplxArray; FFTLen, sign: integer); pascal; external;
procedure _DoCplxFFT2(samps, scTable: PfCplxArray; FFTLen, sign: integer); pascal; external;

procedure _BuildRealOutput(samps, SinTable: PfCplxArray; N,Sign: integer); pascal; external;
{$F-}
{$ENDIF}
{$ENDIF}

{------------------------------------------------------------------------------}
function CalcWindowFunc(Window, idx, Len: integer): Float;
const
   alpha = 5.0;   { Gaussian window parameter }
begin
   { Calculate Windowing function }
   case Window of
          { Hamming }
       1: Result := 0.54-0.46*cos(2*M_PI*idx/Len);
          { Hanning }
       2: Result := 0.5 - 0.5*cos(2*M_PI*idx/Len);
          { Blackman }
       3: Result := 0.42-0.5*cos(2*M_PI*idx/Len)+0.08*cos(4*M_PI*idx/Len);
	    { Gaussian }
       4: Result := exp(-alpha/(Len*Len)*(2*idx-Len)*(2*idx-Len));
            { Welch }
       5: Result := 1 - ((2*idx-Len)/(Len+1))*((2*idx-Len)/(Len+1));
            { Parzen }
       6: Result := 1 - abs((2*idx-Len)/(Len+1));
            { Rectangular }
     else Result := 1;
   end;
end;

{------------------------------------------------------------------------------}
procedure GenWindowTableInt(pTable: PIntArray; Window, Order: integer);
var
   i: integer;
   Len: integer;

begin
   Len := 1 shl Order;
   { Calculate Windowing function }
   for i := 0 to Len-1 do
   begin
      pTable^[i] := Floor(CalcWindowFunc(Window,i,Len)*32767+0.5);
   end;
end;

{------------------------------------------------------------------------------}
procedure GenWindowTableFloat(pTable: PFloatArray; Window, Order: integer);
var
   i: integer;
   Len: integer;

begin
   Len := 1 shl Order;
   { Calculate Windowing function }
   for i := 0 to Len-1 do
   begin
      pTable^[i] := CalcWindowFunc(Window,i,Len);
   end;
end;

{------------------------------------------------------------------------------}
procedure CalcWindowInt(Samps: PIntArray; Window, Order: integer);
var
   i: integer;
   Len: integer;

begin
   Len := 1 shl Order;
   for i := 0 to Len-1 do
   begin
      Samps^[i] := Trunc(Samps^[i]*CalcWindowFunc(Window,i,Len));
   end;
end;

{------------------------------------------------------------------------------}
procedure CalcWindowFloat(Samps: PFloatArray; Window, Order: integer);
var
   i: integer;
   Len: integer;

begin
   Len := 1 shl Order;
   for i := 0 to Len-1 do
   begin
      Samps^[i] := Samps^[i]*CalcWindowFunc(Window,i,Len);
   end;
end;

{------------------------------------------------------------------------------}
procedure GenBitRevTable(pTable: PIntArray; Order: integer);
var
   i, temp, mask: integer;
   Len: integer;

begin
   Len := 1 shl Order;

   for i := 0 to Len-1 do
   begin
      temp := 0;
      mask := Len div 2;
      while mask > 0 do
      begin
         temp := temp shr 1;
         if (i and mask > 0) then temp := temp + Len div 2;
         mask := mask shr 1;
      end;
      pTable^[i] := temp;
   end;
end;

{$IFDEF WIN32}
{------------------------------------------------------------------------------}
function InitCplxFFT(iOrder: integer): PFFTCplx;
begin
   Result := nil;
   if (iOrder > 0) and (iOrder <= 31) then
   begin
      Result := GlobalAllocMem(SizeOf(TFFTCplx));
      if (Result <> nil) then
      with Result^ do
      begin
         Order       := iOrder;
         FFTLen      := 1 shl Order;
         Scale       := 1/FFTLen;

         BitRevTable := GlobalAllocMem(FFTlen*sizeOf(Integer));
         GenBitRevTable(BitRevTable,Order);

         SinTable    := GlobalAllocMem((FFTLen div 2)*sizeOf(TfCplx));
         GenSinTable1(SinTable, Order)
      end;
   end;
end;

{------------------------------------------------------------------------------}
procedure DoneCplxFFT(var pfft: PFFTCplx);
begin
   if (pfft <> nil) then
   begin
      GlobalFreeMem(Pointer(pfft^.BitRevTable));
      GlobalFreeMem(Pointer(pfft^.SinTable));
      GlobalFreeMem(Pointer(pfft));
   end;
end;

{------------------------------------------------------------------------------}
procedure DoCplxFFTa(pfft: PFFTCplx; re, im: PFLoatArray; Sign: integer);
begin
   if (pfft <> nil) and (re <> nil) and (im <> nil) then
   with pfft^ do
   begin
      DoCplxFFT1(re,im,sinTable,FFTLen,Sign);

      fvecBitRev1(re,BitRevTable,Order);
      fvecBitRev1(im,BitRevTable,Order);

      if (Sign < 0) then
      begin
	 fvecMul1(re,Scale,FFTLen);
	 fvecMul1(im,Scale,FFTLen);
      end;
   end;
end;

{------------------------------------------------------------------------------}
procedure DoCplxFFTb(pfft: PFFTCplx; samps: PfCplxArray; Sign: integer);
var
   S: integer;
begin
   if (pfft <> nil) and (samps <> nil) then
   with pfft^ do
   begin
      if Sign = -1 then
         S := -1
      else
         S := 1;

      DoCplxFFT2(samps,sinTable,FFTLen,S);

      fvecBitRev2(samps,BitRevTable,Order);

      if (Sign < 0) then
      begin
	 fvecMul1(Pointer(samps),Scale,2*FFTLen);
      end;
   end;
end;

{==============================================================================}