unadsp.pas

Voice Commnucation Components for Delphi

(*
	----------------------------------------------

	  unaDsp.pas - DSP routines and classes
	  Voice Communicator components version 2.5 Pro

	----------------------------------------------
	  This source code cannot be used without
	  proper permission granted to you as a private
	  person or an entity by the Lake of Soft, Ltd

	  Visit http://lakeofsoft.com/ for details.

	  Copyright (c) 2003 Lake of Soft, Ltd
		     All rights reserved
	----------------------------------------------

	  created by:
		Lake, 31 Oct 2003

	  modified by:
		Lake, Oct 2003

	----------------------------------------------
*)

{$I unaDef.inc}

unit
  unaDsp;

interface

uses
  unaTypes, unaUtils, unaClasses, unaWave;

type
  {DP:CLASS
    FFT implementation.
  }
  unadspFFT = class(unaObject)
  private
    f_bits: unsigned;
    f_channels: unsigned;
    //
    f_steps: unsigned;
    f_windowSize: unsigned;
    //
    f_helper: pUint32Array;
    f_cos_array: pFloatArray;
    f_sin_array: pFloatArray;
    f_dataProxy: pFloatArray;
    //
    f_dataR: pFloatArray;
    f_dataI: pFloatArray;
    //
    procedure samplesToDataProxy(samples: pointer; channel: unsigned);
    procedure dataProxyToRI();
    procedure setSteps(value: unsigned);
  public
    constructor create(windowSize: unsigned = 1024);
    //
    procedure AfterConstruction(); override;
    procedure BeforeDestruction(); override;
    //
    {DP:METHOD
      size should be power of 2.
    }
    procedure setWindowSize(size: unsigned);
    procedure setFormat(format: punaPCMFormat); overload;
    procedure setFormat(bits, channels: unsigned); overload;
    //
    {DP:METHOD
      complex DFT. Results are in dataR, dataI.
      /Works really slow/
    }
    procedure dft_complex_forward(samples: pointer; channel: unsigned = 0);
    {DP:METHOD
      complex FFT. Results are in dataR, dataI.
      /Works faster/
    }
    procedure fft_complex_forward(samples: pointer; channel: unsigned = 0);
    //
    property data: pFloatArray read f_dataProxy;
    //
    property dataR: pFloatArray read f_dataR;
    property dataI: pFloatArray read f_dataI;
    //
    property windowSize: unsigned read f_windowSize write setWindowSize;
    property steps: unsigned read f_steps write setSteps;
  end;


implementation


// -- from Math.pas --
function log2(const x: int64): float;
asm
	fld1		// ST(0) <= 1.0
	fild	x	// ST(1) <= ST(0) <= x
	fyl2x		// ST(1) <= ST(1) * log.2(ST(0))
	fwait
end;

// --  --
function bitReverse(x, steps: unsigned): unsigned;
{
	IN EAX = x
	IN ECX = steps

	OUT EAX = result
}
asm
	mov	ecx, edx
	and	ecx, $1F	// 0..31
	jecxz	@done
	//
	push	ecx
  @loop:
	ror	eax, 1
	rcl	edx, 1
	loop	@loop
	//
	pop	ecx
	//
	neg	ecx
	add	ecx, 32
	shl	edx, cl		// clear high cl bits
	shr	edx, cl
	//
  @done:
	mov	eax, edx
end;

{ unadspFFT }

// --  --
procedure unadspFFT.afterConstruction();
begin
  inherited;
  //
  setWindowSize(f_windowSize);
end;

// --  --
procedure unadspFFT.beforeDestruction();
begin
  inherited;
  //
  mrealloc(f_helper);
  mrealloc(f_cos_array);
  mrealloc(f_sin_array);
  mrealloc(f_dataProxy);
  mrealloc(f_dataR);
  mrealloc(f_dataI);
end;

// --  --
constructor unadspFFT.create(windowSize: unsigned);
begin
  inherited create();
  //
  f_windowSize := windowSize;
  //
  f_helper := nil;
  f_cos_array := nil;
  f_sin_array := nil;
  f_dataProxy := nil;
  f_dataR := nil;
  f_dataI := nil;
  //
  f_bits := 16;
  f_channels := 2;
end;

// --  --
procedure unadspFFT.dataProxyToRI();
var
  i: unsigned;
begin
  for i := 0 to f_windowSize - 1 do begin
    //
    f_dataR[i] := f_dataProxy[i];
    f_dataI[i] := 0;
  end;
end;

// --  --
procedure unadspFFT.dft_complex_forward(samples: pointer; channel: unsigned = 0);
var
  k, i: unsigned;
  sc: float;
begin
  samplesToDataProxy(samples, channel);
  //
  for i := 0 to f_windowSize - 1 do begin
    //
    // Zero REX[ ] and IMX[ ], so they can be used
    // as accumulators during the correlation
    //
    f_dataR[i] := 0;
    f_dataI[i] := 0;
  end;
  //
  for k := 0 to f_windowSize - 1 do begin	// Loop for each value in frequency domain
    //
    sc := 2 * Pi * k / f_windowSize;
    for i := 0 to f_windowSize - 1 do begin	// Correlate with the complex sinusoid, SR & SI
      //
      f_dataR[k] := f_dataR[k] + f_dataProxy[i] * cos(sc * i);
      f_dataI[k] := f_dataI[k] - f_dataProxy[i] * sin(sc * i);
    end;
  end;
end;

// --  --
procedure unadspFFT.fft_complex_forward(samples: pointer; channel: unsigned);
var
  i, j, l, le2, ip: unsigned;
  tr, ti, ur, ui, sr, si: float;
begin
  {
    THE FAST FOURIER TRANSFORM
    Upon entry, N% contains the number of points in the DFT, REX[ ] and
    IMX[ ] contain the real and imaginary parts of the input. Upon return,
    REX[ ] and IMX[ ] contain the DFT output. All signals run from 0 to N%-1.
  }
  samplesToDataProxy(samples, channel);
  dataProxyToRI();
  //
  for i := 0 to f_windowSize - 1 do begin	// Bit reversal sorting
    //
    j := f_helper[i];
    if (i < j) then begin
      //
      ti := f_dataR[i];
      f_dataR[i] := f_dataR[j];
      f_dataR[j] := ti;
    end;
  end;
  //
  for L := 1 to f_steps do begin	// Loop for each stage
    //
    LE2 := (1 shl L) shr 1;
    UR := 1;
    UI := 0;
    //
    SR := f_cos_array[L - 1];	// Calculate sine & cosine values
    SI := f_sin_array[L - 1];
    //
    for j := 1 to LE2 do begin	// Loop for each sub DFT
      //
      I := J - 1;
      while (i < f_windowSize) do begin	// Loop for each butterfly
	//
	IP := I + LE2;
	TR := f_dataR[IP] * UR - f_dataI[IP] * UI;	// Butterfly calculation
	TI := f_dataR[IP] * UI + f_dataI[IP] * UR;
	//
	f_dataR[IP] := f_dataR[I] - TR;
	f_dataI[IP] := f_dataI[I] - TI;
	f_dataR[I] := f_dataR[I] + TR;
	f_dataI[I] := f_dataI[I] + TI;
	//
	inc(i, LE2 shl 1);
      end;
      //
      TR := UR;
      UR := TR * SR - UI * SI;
      UI := TR * SI + UI * SR;
    end;
  end;
end;

// --  --
procedure unadspFFT.samplesToDataProxy(samples: pointer; channel: unsigned);
var
  i, ofs: unsigned;
begin
  ofs := channel;
  //
  for i := 0 to f_windowSize - 1 do begin
    //
    case (f_bits) of

      8:
	f_dataProxy[i] := ($7F - pArray(samples)[ofs]) shl 8;

      16:
	f_dataProxy[i] := pSmallIntArray(samples)[ofs];

      32:
	f_dataProxy[i] := pInt32Array(samples)[ofs];

    end;
    //
    inc(ofs, f_channels);
  end;
end;

// --  --
procedure unadspFFT.setFormat(format: punaPCMFormat);
begin
  setFormat(format.pcmBitsPerSample, format.pcmNumChannels);
end;

// --  --
procedure unadspFFT.setFormat(bits, channels: unsigned);
begin
  f_bits := bits;
  f_channels := channels;
end;

// --  --
procedure unadspFFT.setSteps(value: unsigned);
begin
  if (0 < (value and $1F)) then
    setWindowSize(1 shl (value and $1F));
end;

// --  --
procedure unadspFFT.setWindowSize(size: unsigned);
var
  i: unsigned;
begin
  {
	m = FFT_Step;
	n = 1 << m;
	nv = n >> 1;
  }
  i := round(log2(size));
  //
  if (0 < i) then begin
    //
    f_steps := i;
    f_windowSize := 1 shl f_steps;
    //
    mrealloc(f_helper, f_windowSize * sizeof(f_helper[0]));
    mrealloc(f_cos_array, f_steps * sizeof(f_cos_array[0]));
    mrealloc(f_sin_array, f_steps * sizeof(f_sin_array[0]));
    mrealloc(f_dataProxy, f_windowSize * sizeof(f_dataProxy[0]));
    //
    mrealloc(f_dataR, f_windowSize * sizeof(f_dataR[0]));
    mrealloc(f_dataI, f_windowSize * sizeof(f_dataI[0]));
    //
    for i := 0 to f_steps - 1 do begin
      //
      f_cos_array[i] :=  cos(Pi / (1 shl i));
      f_sin_array[i] := -sin(Pi / (1 shl i));
    end;
    //
    for i := 0 to f_windowSize - 1 do
      f_helper[i] := bitReverse(i, f_steps);
    //
  end;
end;


end.