mmfft.pas

一套及时通讯的原码

function InitRealFFT(iOrder: integer): PFFTReal;
begin
   Result := nil;
   if (iOrder > 0) and (iOrder <= 31) then
   begin
      Result := GlobalAllocMem(SizeOf(TFFTReal));
      if (Result <> nil) then
      with Result^ do
      begin
         Order       := iOrder;
         FFTLen      := 1 shl Order;
         Scale       := 1/FFTLen;

         BitRevTable := GlobalAllocMem((FFTlen div 2)*sizeOf(Integer));
         GenBitRevTable(BitRevTable,Order-1);

         SinTable1   := GlobalAllocMem((FFTLen div 4)*sizeOf(TfCplx));
         GenSinTable1(SinTable1, Order-1);

         SinTable2   := GlobalAllocMem((FFTLen div 4)*sizeOf(TfCplx));
         GenSinTable2(SinTable2, Order-1);
      end;
   end;
end;

{------------------------------------------------------------------------------}
procedure DoneRealFFT(var pfft: PFFTReal);
begin
   if (pfft <> nil) then
   begin
      GlobalFreeMem(Pointer(pfft^.BitRevTable));
      GlobalFreeMem(Pointer(pfft^.SinTable1));
      GlobalFreeMem(Pointer(pfft^.SinTable2));
      GlobalFreeMem(Pointer(pfft));
   end;
end;

{------------------------------------------------------------------------------}
procedure DoRealFFT(pfft: PFFTReal; Samps: PFLoatArray; Sign: integer);
var
   i: integer;

begin
   if (pfft <> nil) and (Samps <> nil) then
   with pfft^ do
   begin
      Sign := MinMax(Sign,-1,1);

      DoCplxFFT2(Pointer(Samps),sinTable1,FFTLen div 2,Sign);

      fvecBitRev2(Pointer(samps),BitRevTable,Order-1);

      i := 1;
      Samps^[0] := Samps^[0]+Samps^[i];
      Samps^[i] := 0;

      BuildRealOutput(Pointer(Samps),SinTable2, FFTLen div 2,Sign);

      if (Sign < 0) then
      begin
	 fvecMul1(Pointer(samps),Scale,FFTLen);
      end;
   end;
end;
{$ELSE}

{-------------------------------------------------------------------------}
procedure CalcWindow(Buffer: PIntArray; Window: TMMFFTWindow; FFTLen: integer);
const
   alpha = 5.0;   { Gaussian window parameter }
var
   i: integer;
   val: Float;

begin
   { Calculate FFT Windowing function }
   for i := 0 to FFTLen-1 do
   begin
      case ord(Window) of
            { Hamming }
	 1: val := 0.54-0.46*cos(2*M_PI*i/FFTLen);
            { Hanning }
         2: val := 0.5 - 0.5*cos(2*M_PI*i/FFTLen);
            { Blackman }
         3: val := 0.42-0.5*cos(2*M_PI*i/FFTLen)+0.08*cos(4*M_PI*i/FFTLen);
	    { Gaussian }
	 4: val := exp(-alpha/(Long(FFTLen)*FFTLen)*(2*i-FFTLen)*(2*i-FFTLen));
            { Welch }
	 5: val := 1 - ((2*i-FFTLen)/(FFTLen+1))*((2*i-FFTLen)/(FFTlen+1));
            { Parzen }
	 6: val := 1 - abs((2*i-FFTLen)/(FFTLen+1));
            { Rectangular }
       else val := 1;
      end;
      Buffer^[i] := Floor(val*32767+0.5);
   end;
end;

{-- TMMFFT ---------------------------------------------------------------}
constructor TMMFFT.Create;
begin
   inherited Create;

   FPoints := 0;
   FBitReversed := nil;
   FSinTable := nil;
   FFTLength := 128;
end;

{-- TMMFFT ---------------------------------------------------------------}
destructor TMMFFT.Destroy;
begin
   DoneFFT;
   inherited destroy;
end;

{-- TMMFFT ---------------------------------------------------------------}
procedure TMMFFT.SetFFTLen(FFTLen: integer);
begin
   DoneFFT;
   InitFFT(FFTLen);
end;

{-------------------------------------------------------------------------}
{ Initialize the Sine table and Twiddle pointers (bit-reversed pointers)
  for the FFT routine. }
{-- TMMFFT ---------------------------------------------------------------}
procedure TMMFFT.InitFFT(FFTLen: integer);
var
   i, temp, mask: integer;
   s,c: Real;

begin
   { FFT size is only half the number of data points.
     The full FFT output can be reconstructed from this FFT's output.
     (This optimization can be made since the data is real.) }
   FPoints := FFTlen;

   FSinTable := GlobalAllocMem(FPoints*sizeOf(SmallInt));
   FBitReversed := GlobalAllocMem(FPoints div 2*sizeOf(Integer));

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

   for i := 0 to (FPoints div 2)-1 do
   begin
      s := (-32768.0 * sin(2*M_PI*i/(FPoints))+0.5);
      c := (-32768.0 * cos(2*M_PI*i/(FPoints))+0.5);
      if (s > 32767.5) then s := 32767;
      if (c > 32767.5) then c := 32767;
      FSinTable^[FBitReversed^[i]]   := Trunc(s);
      FSinTable^[FBitReversed^[i]+1] := Trunc(c);
   end;
end;

{-------------------------------------------------------------------------}
{ Free up the memory allocated for Sin table and Twiddle Pointers }
{-- TMMFFT ---------------------------------------------------------------}
procedure TMMFFT.DoneFFT;
begin
   if FPoints <> 0 then
   begin
      GlobalFreeMem(Pointer(FBitReversed));
      GlobalFreeMem(Pointer(FSinTable));
      FPoints := 0;
   end;
end;

{-------------------------------------------------------------------------}
{ Actual FFT routine.  Must call InitFFT(FFTLen) first! }
{-------------------------------------------------------------------------}
{$IFDEF USEASM}
{$L MMFFT16.OBJ}
procedure CalcFFTASM(Buffer: PSmallInt);Far;{$IFDEF WIN32}pascal;{$ENDIF}external;
Var
   Points     : integer;
   SinTable   : PSmallArray;
   BitReversed: PIntArray;

{-- TMMFFT ---------------------------------------------------------------}
procedure TMMFFT.CalcFFT(Buffer: PSmallInt);
begin
   MMFFT.Points := FPoints;
   MMFFT.SinTable := FSinTable;
   MMFFT.BitReversed := FBitReversed;
   CalcFFTASM(Buffer);
end;
{$ELSE}
procedure TMMFFT.CalcFFT(Buffer: PSmallInt);
var
   A,B: PSmallInt;
   sptr: PSmallInt;
   endptr1, endptr2: PSmallInt;
   br1,br2: PInteger;
   HRplus,HRminus,HIplus,HIminus: Long;

   v1,v2: Long;
   temp1,temp2: Long;
   sin, cos: SmallInt;
   ButterfliesPerGroup: integer;

begin
   ButterfliesPerGroup := FPoints div 4;

   endptr1 := PSmallInt(PChar(Buffer) + FPoints * sizeOf(SmallInt));

   { Butterfly:
      Ain-----Aout
          \ /
          / \
      Bin-----Bout }

   while (ButterfliesPerGroup > 0) do
   begin
      A := Buffer;
      B := PSmallInt(PChar(Buffer) + 2*ButterfliesPerGroup * sizeOf(SmallInt));
      sptr := PSmallInt(FSinTable);

      while PChar(A) < PChar(endptr1) do
      begin
         sin := sptr^;
         cos := PSmallInt(PChar(sptr) + sizeOf(SmallInt))^;

         endptr2 := B;
         while PChar(A) < PChar(endptr2) do
         begin
            v1 := (Long(B^)*cos+Long(PSmallInt(PChar(B)+sizeOf(SmallInt))^)*sin) shr 15;
            v2 := (Long(B^)*sin-Long(PSmallInt(PChar(B)+sizeOf(SmallInt))^)*cos) shr 15;
            B^ := (A^ + v1) shr 1;
            A^ := (B^ - v1);
            inc(A);
            inc(B);
            B^ := (A^ - v2) shr 1;
            A^ := (B^ + v2);
            inc(A);
            inc(B);
         end;
	 A := B;
	 inc(B, ButterfliesPerGroup*2);
	 inc(sptr, 2);
      end;
      ButterfliesPerGroup := ButterfliesPerGroup shr 1;
   end;

   { Massage output to get the output for a real input sequence. }

   br1 := PInteger(PChar(FBitReversed)+sizeOf(Integer));
   br2 := PInteger(PChar(FBitReversed)+((FPoints div 2)-1)*sizeOf(Integer));

   while PChar(br1) <= PChar(br2) do
   begin
      sin := FSinTable^[br1^];
      cos := FSinTable^[br1^+1];
      A := PSmallInt(PChar(Buffer) + br1^ * sizeOf(SmallInt));
      B := PSmallInt(PChar(Buffer) + br2^ * sizeOf(SmallInt));
      HRminus := A^ - B^;
      HRplus := HRminus + (B^ shl 1);
      HIminus := PSmallInt(PChar(A)+sizeOf(SmallInt))^ - PSmallInt(PChar(B)+sizeOf(SmallInt))^;
      HIplus := HIminus + (PSmallInt(PChar(B)+sizeOf(SmallInt))^ shl 1);
      temp1  := (sin*HRminus - cos*HIplus) shr 15;
      temp2  := (cos*HRminus + sin*HIplus) shr 15;
      A^ := (HRplus + temp1) shr 1;
      B^ := A^ - temp1;
      inc(A);
      inc(B);
      A^ := (HIminus + temp2) shr 1;
      B^ := A^ - HIminus;

      inc(br1);
      dec(br2);
   end;

   { Handle DC bin separately }
   Buffer^ := PSmallInt(PChar(Buffer)+sizeOf(SmallInt))^;
   PSmallInt(PChar(Buffer)+sizeOf(SmallInt))^ := 0;
end;
{$ENDIF}
{$ENDIF}

{$IFDEF USEDLL}
var
   DLLHandle : THandle;
{$ENDIF}

{========================================================================}
procedure InitMMUtils;
begin
   {$IFDEF USEDLL}
   //DLLHandle := GetModuleHandle(MMUTILDLLKeyName);

   //if DLLHandle < HINSTANCE_ERROR then
      DLLHandle := GetModuleHandle(MMUTILDLLName);

   if DLLHandle >= HINSTANCE_ERROR then
   begin
      @fvecMul1       := GetProcAddress(DLLHandle,'_fvecMul1');
      @fvecMul2       := GetProcAddress(DLLHandle,'_fvecMul2');
      @fvecZero1      := GetProcAddress(DLLHandle,'_fvecZero1');
      @fvecZero2      := GetProcAddress(DLLHandle,'_fvecZero2');
      @fvecBitRev1    := GetProcAddress(DLLHandle,'_fvecBitRev1');
      @fvecBitRev2    := GetProcAddress(DLLHandle,'_fvecBitRev2');
      @GenSinTable1   := GetProcAddress(DLLHandle,'_GenSinTable1');
      @GenSinTable2   := GetProcAddress(DLLHandle,'_GenSinTable2');
      @DoCplxFFT1     := GetProcAddress(DLLHandle,'_DoCplxFFT1');
      @DoCplxFFT2     := GetProcAddress(DLLHandle,'_DoCplxFFT2');
      @BuildRealOutput:= GetProcAddress(DLLHandle,'_BuildRealOutput');
   end
   else Halt;

   {$ELSE}

   fvecMul1       := _fvecMul1;
   fvecMul2       := _fvecMul2;
   fvecZero1      := _fvecZero1;
   fvecZero2      := _fvecZero2;
   fvecBitRev1    := _fvecBitRev1;
   fvecBitRev2    := _fvecBitRev2;
   GenSinTable1   := _GenSinTable1;
   GenSinTable2   := _GenSinTable2;
   DoCplxFFT1     := _DoCplxFFT1;
   DoCplxFFT2     := _DoCplxFFT2;
   BuildRealOutput:= _BuildRealOutput;
   {$ENDIF}
end;

initialization
   InitMMUtils;
end.