comparison quickspec.cpp @ 2:fc19d45615d1

* Make all file names lower-case to avoid case ambiguity
(some includes differed in case from the filenames they were
trying to include).  Also replace MinGW-specific mem.h with
string.h

+//---------------------------------------------------------------------------
+
+#include <math.h>
+#include <memory.h>
+#include "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
+
+