annotate quickspec.cpp @ 13:de3961f74f30 tip

Add Linux/gcc Makefile; build fix

/*
    Harmonic sinusoidal modelling and tools

    C++ code package for harmonic sinusoidal modelling and relevant signal processing.
    Centre for Digital Music, Queen Mary, University of London.
    This file copyright 2011 Wen Xue.

    This program is free software; you can redistribute it and/or
    modify it under the terms of the GNU General Public License as
    published by the Free Software Foundation; either version 2 of the
    License, or (at your option) any later version.
*/
//---------------------------------------------------------------------------
#include <math.h>
#include <memory.h>
#include "quickspec.h"

/** \file quickspec.h */

//---------------------------------------------------------------------------
/**
  method TQuickSpectrogram::TQuickSpectrogram:

  In: AParent: pointer argument for calling G, if G is specified
      AnId: integer argument for calling G, if G is specified
      G: pointer to a function that supplies buffers used for FFT, 0 by default
      Ausex: set if complete complex spectrogram is to be buffered and accessible
      Auseph: set if phase spectrogram is to be buffered and accessible
*/
__fastcall TQuickSpectrogram::TQuickSpectrogram(void* AParent, int AnId, bool Ausex, bool Auseph, GetBuf G)
{
  memset(this, 0, sizeof(TQuickSpectrogram));
  Parent=AParent;
  Id=AnId;
  usex=Ausex;
  useph=Auseph;
  GetFFTBuffers=G;
  BufferSize=QSpec_BufferSize;
  fwt=wtRectangle;
  Wid=1024;
  Offst=512;
}//TQuickSpectrogram

//TQuickSpectrogram::~TQuickSpectrogram
__fastcall TQuickSpectrogram::~TQuickSpectrogram()
{
  FreeBuffers();
  free8(fw);
	free8(fwin);
  free(fhbi);
}//~TQuickSpectrogram

//---------------------------------------------------------------------------
/**
  method TQuickSpectrogram::A: accesses amplitude spectrogram by frame

  In: fr: frame index, 0-based

  Returns pointer to amplitude spectrum of the fr'th frame, NULL if N/A
*/
QSPEC_FORMAT* __fastcall TQuickSpectrogram::A(int fr)
{
  if (Capacity==0) SetFrCapacity((DataLength-Wid)/Offst+2);
  if (fr<0 || fr>=Capacity) return NULL;
  if (Frame[fr]<0 || !Valid[fr]) CalculateSpectrum(fr);
  return fA[Frame[fr]];
}//A

//---------------------------------------------------------------------------
/**
  method TQuickSpectrogram::AddBuffer: increases internal buffer by BufferSize frames. Allocated buffers
  beyond Capacity frames will not be indexed or used by TQuickSpectrogram.
*/
void TQuickSpectrogram::AddBuffer()
{
  int base=1, Dim=Wid/2+1;
  if (usex) base+=2;
  if (useph) base+=1;
  QSPEC_FORMAT* newbuffer=(QSPEC_FORMAT*)malloc(sizeof(QSPEC_FORMAT)*Dim*BufferSize*base);
  int fr0=BufferCount*BufferSize;
  for (int i=0; i<BufferSize; i++)
  {
    int fr=fr0+i;
    if (fr<Capacity)
    {
      fA[fr]=&newbuffer[i*Dim];
      int base=1;
      if (usex) fSpec[fr]=(cmplx<QSPEC_FORMAT>*)&newbuffer[(BufferSize+i*2)*Dim], base+=2;
      if (useph) fPh[fr]=&newbuffer[(BufferSize*base+i)*Dim];
    }
    else break;
  }
  BufferCount++;
}//AddBuffer

/**
  method TQuickSpectrogram::AddBuffer: increase internal buffer by a multiple of BufferSize so that
  it will be enough to host another AddFrCount frames.
*/
void TQuickSpectrogram::AddBuffer(int AddFrCount)
{
  while (FrCount+AddFrCount>BufferSize*BufferCount) AddBuffer();
}//AddBuffer

//---------------------------------------------------------------------------
/**
  function IntToDouble: copy content of integer array to double array

  In: in: pointer to integer array
      BytesPerSample: number of bytes each integer takes
      Count: size of integer array, in integers
  Out: vector out[Count].

  No return value.
*/
void IntToDouble(double* out, void* in, int BytesPerSample, int Count)
{
  if (BytesPerSample==1){unsigned char* in8=(unsigned char*)in; for (int k=0; k<Count; k++) *(out++)=*(in8++)-128.0;}
  else if (BytesPerSample==2) {__int16* in16=(__int16*)in; for (int k=0; k<Count; k++) *(out++)=*(in16++);}
  else {__pint24 in24=(__pint24)in; for (int k=0; k<Count; k++) *(out++)=*(in24++);}
}//IntToDouble

/**
  function CalculateSpectrum: calculate spectrum of a signal in integer format

  In: Data[Wid]: integer array hosting waveform data
      BytesPerSample: number of bytes each integer in Data[] takes
      win[Wid]: window function used for computing spectrum
      w[Wid/2], x[Wid], hbi[Wid/2]: FFT buffers
  Out:  x[Wid]: complex spectrum
        Amp[Wid/2+1]: amplitude spectrum
        Arg[Wid/2+1]: phase spectrum, optional

  No return value.
*/
void CalculateSpectrum(void* Data, int BytesPerSample, double* win, QSPEC_FORMAT* Amp, QSPEC_FORMAT* Arg, int Wid, cdouble* w, cdouble* x, int* hbi)
{
  if (BytesPerSample==2) RFFTCW((__int16*)Data, win, 0, 0, Log2(Wid), w, x, hbi);
  else {IntToDouble((double*)x, Data, BytesPerSample, Wid); RFFTCW((double*)x, win, 0, 0, Log2(Wid), w, x, hbi);}
  for (int j=0; j<=Wid/2; j++)
  {
    Amp[j]=sqrt(x[j].x*x[j].x+x[j].y*x[j].y);
    if (Arg) Arg[j]=(x[j].y==0 && x[j].x==0)?0:atan2(x[j].y, x[j].x);
  }
}//CalculateSpectrum

/**
  function CalculateSpectrum: calculate spectrum of a signal in integer format, allowing the signal
  length $eff be shorter than the DFT size Wid.

  In: Data[eff]: integer array hosting waveform data
      BytesPerSample: number of bytes each integer in Data[] takes
      win[Wid]: window function used for computing spectrum
      w[Wid/2], x[Wid], hbi[Wid/2]: FFT buffers
  Out:  x[Wid]: complex spectrum
        Amp[Wid/2+1]: amplitude spectrum
        Arg[Wid/2+1]: phase spectrum, optional

  No return value.
*/
void CalculateSpectrum(void* Data, int BytesPerSample, double* win, QSPEC_FORMAT* Amp, QSPEC_FORMAT* Arg, int Wid, int eff, cdouble* w, cdouble* x, int* hbi)
{
  if (eff<=0)
  {
    memset(Amp, 0, sizeof(double)*(Wid/2+1));
    if (Arg) memset(Arg, 0, sizeof(double)*(Wid/2+1));
  }
  else if (eff<Wid)
  {
    double* doublex=(double*)x;
    IntToDouble(doublex, Data, BytesPerSample, eff); memset(&doublex[eff], 0, sizeof(double)*(Wid-eff));
    RFFTCW(doublex, win, 0, 0, Log2(Wid), w, x, hbi);
    for (int j=0; j<=Wid/2; j++)
    {
      Amp[j]=sqrt(x[j].x*x[j].x+x[j].y*x[j].y);
      if (Arg) Arg[j]=(x[j].y==0 && x[j].x==0)?0:atan2(x[j].y, x[j].x);
    }
  }
  else
    CalculateSpectrum(Data, BytesPerSample, win, Amp, Arg, Wid, w, x, hbi);
}//CalculateSpectrum

/**
  method TQuickSpectrogram::CalculateSpectrum: computes spectrogram at fr'th frame.

  In: fr: index to the frame whose spectrum is to be computed. fr must be between 0 and Capacity-1.
*/
void __fastcall TQuickSpectrogram::CalculateSpectrum(int fr)
{
  cdouble *w, *x;
  double* win;
  int* hbi;

  //obtain FFT buffers win (window function), w (twiddle factors), x (data buffer),
  //hbi (half-size bit-inversed integer table)
  if (GetFFTBuffers)  //then use external buffers provided through GetFFTBuffers
      GetFFTBuffers(Id, w, x, win, hbi, Parent);
  else  //then use internal buffers
  {
    if (Wid!=fWid)
    { //then update internal buffers to the new window size
      free8(fw); free(fhbi);
      fw=(cdouble*)malloc8(sizeof(cdouble)*Wid*1.5); SetTwiddleFactors(Wid, fw);
      fx=&fw[Wid/2];
      fhbi=CreateBitInvTable(Log2(Wid)-1);
    }
    if (Wid!=fWid || WinType!=fwt || WinParam!=fwdp)
    { //then update internal window function to the new window type
			fwin=NewWindow8(WinType, Wid, 0, &WinParam);
      fwt=WinType; fwdp=WinParam;
    }
    fWid=Wid;

    //pick up the internal buffers
		w=fw, x=fx, win=fwin, hbi=fhbi;
  }

  //obtain the index of this frame in internal storage
  if (Frame[fr]<0) {AddBuffer(1); Frame[fr]=FrCount; FrCount++;}
  int realfr=Frame[fr];

  //obtain the pointer to this frame's phase spectrum
  QSPEC_FORMAT *lph=useph?fPh[realfr]:NULL;
  //ontain the pointer to this frame's complex spectrum
  cmplx<QSPEC_FORMAT>* lX=usex?fSpec[realfr]:NULL;
  //choose the buffer actually used for FFT - use lX if it is specified as complex double array
  //because it saves unnecessary data copying operations
  cdouble *lx=(usex && sizeof(QSPEC_FORMAT)==sizeof(double))?(cdouble*)lX:x;

  //Actual FFT
  if (fr*Offst+Wid<=DataLength)
    ::CalculateSpectrum(&((char*)Data)[fr*Offst*BytesPerSample], BytesPerSample, win, fA[realfr], lph, Wid, w, lx, hbi);
  else
    ::CalculateSpectrum(&((char*)Data)[fr*Offst*BytesPerSample], BytesPerSample, win, fA[realfr], lph, Wid, DataLength-fr*Offst, w, x, hbi);

  //optional data copy from x to lX
  if (usex && lx==x) for (int i=0; i<Wid/2+1; i++) lX[i]=x[i];

  //tag this frame as computed and valid
  Valid[fr]=1;
}//CalculateSpectrum

//---------------------------------------------------------------------------
/**
  method TQuickSpectrogram::FreeBuffers: discards all computed spectra and free all internal buffers.
  This returns the TQuickSpectrogram to its initial state before any frame is accessed. After calling
  FreeBuffers() all frames will be recomputed when they are accessed.
*/
void TQuickSpectrogram::FreeBuffers()
{
  if (fA)
  {
    for (int i=0; i<BufferCount; i++) free(fA[i*BufferSize]);
    FrCount=BufferCount=Capacity=0;
    free(Frame); free(Valid);
    free(fA);
    Frame=Valid=0, fA=0;
    if (useph) {free(fPh); fPh=0;}
    if (usex) {free(fSpec); fSpec=0;}
  }
}//FreeBuffers

//---------------------------------------------------------------------------
/**
  method TQuickSpectrogram::Invalidate: renders all frames that have overlap with interval [From, To],
  measured in samples, as invalid. Invalid frames are recomputed when they are accessed again.

  In: [From, To]: an interval spectrogram over which needs to be updated.

  Returns the number of allocated frames affected, no matter if they were valid.
*/
int TQuickSpectrogram::Invalidate(int From, int To)
{
  int result=0;
  if (Frame)
  {
    int fr1=ceil((From-Wid+1.0)/Offst), fr2=floor(1.0*To/Offst);
    if (fr1<0) fr1=0;
    if (fr2>=Capacity) fr2=Capacity-1;
    for (int fr=fr1; fr<=fr2; fr++) if (Frame[fr]>=0) Valid[fr]=false, result++;
  }
  return result;
}//Invalidate

//---------------------------------------------------------------------------
/**
  method TQuickSpectrogram::Ph: accesses phase spectrogram by frame

  In: fr: frame index, 0-based

  Returns pointer to phase spectrum of the fr'th frame, NULL if N/A
*/
QSPEC_FORMAT* __fastcall TQuickSpectrogram::Ph(int fr)
{
  if (Capacity==0) SetFrCapacity((DataLength-Wid)/Offst+2);
  if (fr<0 || fr>=Capacity) return NULL;
  if (Frame[fr]<0 || !Valid[fr]) CalculateSpectrum(fr);
  return fPh[Frame[fr]];
}//Ph

//---------------------------------------------------------------------------
/**
  method TQuickSpectrogram::SetFrCapacity: sets the capacity, i.e. the maximal number of frames handled
  by this TQuickSpectrogram.

  In: AnFrCapacity: the new Capacity, in frames

  This method should not be called to set Capacity to a smaller value.
*/
void TQuickSpectrogram::SetFrCapacity(int AnFrCapacity)
{
  //adjusting the size of index and validity arrays
  Frame=(int*)realloc(Frame, sizeof(int)*AnFrCapacity);
  Valid=(int*)realloc(Valid, sizeof(int)*AnFrCapacity);

  //
  fA=(QSPEC_FORMAT**)realloc(fA, sizeof(QSPEC_FORMAT*)*AnFrCapacity);
  if (usex) fSpec=(cmplx<QSPEC_FORMAT>**)realloc(fSpec, sizeof(cmplx<QSPEC_FORMAT>*)*AnFrCapacity);
  if (useph) fPh=(QSPEC_FORMAT**)realloc(fPh, sizeof(QSPEC_FORMAT*)*AnFrCapacity);
  if (AnFrCapacity>Capacity)
  {
    memset(&Frame[Capacity], 0xFF, sizeof(int)*(AnFrCapacity-Capacity));
    memset(&Valid[Capacity], 0x00, sizeof(int)*(AnFrCapacity-Capacity));

    if (Capacity<BufferCount*BufferSize)
    {
      for (int fr=Capacity; fr<AnFrCapacity; fr++)
      {
        int bufferno=fr/BufferSize;
        if (bufferno<BufferCount)
        {
          QSPEC_FORMAT* thisbuffer=fA[BufferSize*bufferno];
          int lfr=fr%BufferSize, base=1, Dim=Wid/2+1;
          fA[fr]=&thisbuffer[lfr*Dim];
          if (usex) fSpec[fr]=(cmplx<QSPEC_FORMAT>*)(&thisbuffer[(BufferSize+lfr*2)*Dim]), base+=2;
          if (useph) fPh[fr]=&thisbuffer[(BufferSize*base+lfr)*Dim];
        }
        else break;
      }
    }
  }
  Capacity=AnFrCapacity;
}//SetFrCapacity

//---------------------------------------------------------------------------
/**
  method TQuickSpectrogram::Ph: accesses complex spectrogram by frame

  In: fr: frame index, 0-based

  Returns pointer to complex spectrum of the fr'th frame, NULL if N/A
*/
cmplx<QSPEC_FORMAT>* __fastcall TQuickSpectrogram::Spec(int fr)
{
  if (Capacity==0) SetFrCapacity((DataLength-Wid)/Offst+2);
  if (fr<0 || fr>=Capacity) return NULL;
  if (Frame[fr]<0 || !Valid[fr]) CalculateSpectrum(fr);
  return fSpec[Frame[fr]];
}//Spec