xue@11: /*
xue@11:     Harmonic sinusoidal modelling and tools
xue@11: 
xue@11:     C++ code package for harmonic sinusoidal modelling and relevant signal processing.
xue@11:     Centre for Digital Music, Queen Mary, University of London.
xue@11:     This file copyright 2011 Wen Xue.
xue@11: 
xue@11:     This program is free software; you can redistribute it and/or
xue@11:     modify it under the terms of the GNU General Public License as
xue@11:     published by the Free Software Foundation; either version 2 of the
xue@11:     License, or (at your option) any later version.
xue@11: */
xue@1: //---------------------------------------------------------------------------
xue@1: #include <math.h>
xue@1: #include <memory.h>
Chris@2: #include "quickspec.h"
xue@1: 
Chris@5: /** \file quickspec.h */
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   method TQuickSpectrogram::TQuickSpectrogram:
xue@1: 
xue@1:   In: AParent: pointer argument for calling G, if G is specified
xue@1:       AnId: integer argument for calling G, if G is specified
xue@1:       G: pointer to a function that supplies buffers used for FFT, 0 by default
xue@1:       Ausex: set if complete complex spectrogram is to be buffered and accessible
xue@1:       Auseph: set if phase spectrogram is to be buffered and accessible
xue@1: */
xue@1: __fastcall TQuickSpectrogram::TQuickSpectrogram(void* AParent, int AnId, bool Ausex, bool Auseph, GetBuf G)
xue@1: {
xue@1:   memset(this, 0, sizeof(TQuickSpectrogram));
xue@1:   Parent=AParent;
xue@1:   Id=AnId;
xue@1:   usex=Ausex;
xue@1:   useph=Auseph;
xue@1:   GetFFTBuffers=G;
xue@1:   BufferSize=QSpec_BufferSize;
xue@1:   fwt=wtRectangle;
xue@1:   Wid=1024;
xue@1:   Offst=512;
xue@1: }//TQuickSpectrogram
xue@1: 
xue@1: //TQuickSpectrogram::~TQuickSpectrogram
xue@1: __fastcall TQuickSpectrogram::~TQuickSpectrogram()
xue@1: {
xue@1:   FreeBuffers();
xue@1:   free8(fw);
xue@1: 	free8(fwin);
xue@1:   free(fhbi);
xue@1: }//~TQuickSpectrogram
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   method TQuickSpectrogram::A: accesses amplitude spectrogram by frame
xue@1: 
xue@1:   In: fr: frame index, 0-based
xue@1: 
xue@1:   Returns pointer to amplitude spectrum of the fr'th frame, NULL if N/A
xue@1: */
xue@1: QSPEC_FORMAT* __fastcall TQuickSpectrogram::A(int fr)
xue@1: {
xue@1:   if (Capacity==0) SetFrCapacity((DataLength-Wid)/Offst+2);
xue@1:   if (fr<0 || fr>=Capacity) return NULL;
xue@1:   if (Frame[fr]<0 || !Valid[fr]) CalculateSpectrum(fr);
xue@1:   return fA[Frame[fr]];
xue@1: }//A
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   method TQuickSpectrogram::AddBuffer: increases internal buffer by BufferSize frames. Allocated buffers
xue@1:   beyond Capacity frames will not be indexed or used by TQuickSpectrogram.
xue@1: */
xue@1: void TQuickSpectrogram::AddBuffer()
xue@1: {
xue@1:   int base=1, Dim=Wid/2+1;
xue@1:   if (usex) base+=2;
xue@1:   if (useph) base+=1;
xue@1:   QSPEC_FORMAT* newbuffer=(QSPEC_FORMAT*)malloc(sizeof(QSPEC_FORMAT)*Dim*BufferSize*base);
xue@1:   int fr0=BufferCount*BufferSize;
xue@1:   for (int i=0; i<BufferSize; i++)
xue@1:   {
xue@1:     int fr=fr0+i;
xue@1:     if (fr<Capacity)
xue@1:     {
xue@1:       fA[fr]=&newbuffer[i*Dim];
xue@1:       int base=1;
xue@1:       if (usex) fSpec[fr]=(cmplx<QSPEC_FORMAT>*)&newbuffer[(BufferSize+i*2)*Dim], base+=2;
xue@1:       if (useph) fPh[fr]=&newbuffer[(BufferSize*base+i)*Dim];
xue@1:     }
xue@1:     else break;
xue@1:   }
xue@1:   BufferCount++;
xue@1: }//AddBuffer
xue@1: 
Chris@5: /**
xue@1:   method TQuickSpectrogram::AddBuffer: increase internal buffer by a multiple of BufferSize so that
xue@1:   it will be enough to host another AddFrCount frames.
xue@1: */
xue@1: void TQuickSpectrogram::AddBuffer(int AddFrCount)
xue@1: {
xue@1:   while (FrCount+AddFrCount>BufferSize*BufferCount) AddBuffer();
xue@1: }//AddBuffer
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function IntToDouble: copy content of integer array to double array
xue@1: 
xue@1:   In: in: pointer to integer array
xue@1:       BytesPerSample: number of bytes each integer takes
xue@1:       Count: size of integer array, in integers
xue@1:   Out: vector out[Count].
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void IntToDouble(double* out, void* in, int BytesPerSample, int Count)
xue@1: {
xue@1:   if (BytesPerSample==1){unsigned char* in8=(unsigned char*)in; for (int k=0; k<Count; k++) *(out++)=*(in8++)-128.0;}
xue@1:   else if (BytesPerSample==2) {__int16* in16=(__int16*)in; for (int k=0; k<Count; k++) *(out++)=*(in16++);}
xue@1:   else {__pint24 in24=(__pint24)in; for (int k=0; k<Count; k++) *(out++)=*(in24++);}
xue@1: }//IntToDouble
xue@1: 
Chris@5: /**
xue@1:   function CalculateSpectrum: calculate spectrum of a signal in integer format
xue@1: 
xue@1:   In: Data[Wid]: integer array hosting waveform data
xue@1:       BytesPerSample: number of bytes each integer in Data[] takes
xue@1:       win[Wid]: window function used for computing spectrum
xue@1:       w[Wid/2], x[Wid], hbi[Wid/2]: FFT buffers
xue@1:   Out:  x[Wid]: complex spectrum
xue@1:         Amp[Wid/2+1]: amplitude spectrum
xue@1:         Arg[Wid/2+1]: phase spectrum, optional
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void CalculateSpectrum(void* Data, int BytesPerSample, double* win, QSPEC_FORMAT* Amp, QSPEC_FORMAT* Arg, int Wid, cdouble* w, cdouble* x, int* hbi)
xue@1: {
xue@11:   if (BytesPerSample==2) RFFTCW((__int16*)Data, win, 0, 0, Log2(Wid), w, x, hbi);
xue@11:   else {IntToDouble((double*)x, Data, BytesPerSample, Wid); RFFTCW((double*)x, win, 0, 0, Log2(Wid), w, x, hbi);}
xue@1:   for (int j=0; j<=Wid/2; j++)
xue@1:   {
xue@1:     Amp[j]=sqrt(x[j].x*x[j].x+x[j].y*x[j].y);
xue@1:     if (Arg) Arg[j]=(x[j].y==0 && x[j].x==0)?0:atan2(x[j].y, x[j].x);
xue@1:   }
xue@1: }//CalculateSpectrum
xue@1: 
Chris@5: /**
xue@1:   function CalculateSpectrum: calculate spectrum of a signal in integer format, allowing the signal
xue@1:   length $eff be shorter than the DFT size Wid.
xue@1: 
xue@1:   In: Data[eff]: integer array hosting waveform data
xue@1:       BytesPerSample: number of bytes each integer in Data[] takes
xue@1:       win[Wid]: window function used for computing spectrum
xue@1:       w[Wid/2], x[Wid], hbi[Wid/2]: FFT buffers
xue@1:   Out:  x[Wid]: complex spectrum
xue@1:         Amp[Wid/2+1]: amplitude spectrum
xue@1:         Arg[Wid/2+1]: phase spectrum, optional
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void CalculateSpectrum(void* Data, int BytesPerSample, double* win, QSPEC_FORMAT* Amp, QSPEC_FORMAT* Arg, int Wid, int eff, cdouble* w, cdouble* x, int* hbi)
xue@1: {
xue@1:   if (eff<=0)
xue@1:   {
xue@1:     memset(Amp, 0, sizeof(double)*(Wid/2+1));
xue@1:     if (Arg) memset(Arg, 0, sizeof(double)*(Wid/2+1));
xue@1:   }
xue@1:   else if (eff<Wid)
xue@1:   {
xue@1:     double* doublex=(double*)x;
xue@1:     IntToDouble(doublex, Data, BytesPerSample, eff); memset(&doublex[eff], 0, sizeof(double)*(Wid-eff));
xue@11:     RFFTCW(doublex, win, 0, 0, Log2(Wid), w, x, hbi);
xue@1:     for (int j=0; j<=Wid/2; j++)
xue@1:     {
xue@1:       Amp[j]=sqrt(x[j].x*x[j].x+x[j].y*x[j].y);
xue@1:       if (Arg) Arg[j]=(x[j].y==0 && x[j].x==0)?0:atan2(x[j].y, x[j].x);
xue@1:     }
xue@1:   }
xue@1:   else
xue@1:     CalculateSpectrum(Data, BytesPerSample, win, Amp, Arg, Wid, w, x, hbi);
xue@1: }//CalculateSpectrum
xue@1: 
Chris@5: /**
xue@1:   method TQuickSpectrogram::CalculateSpectrum: computes spectrogram at fr'th frame.
xue@1: 
xue@1:   In: fr: index to the frame whose spectrum is to be computed. fr must be between 0 and Capacity-1.
xue@1: */
xue@1: void __fastcall TQuickSpectrogram::CalculateSpectrum(int fr)
xue@1: {
xue@1:   cdouble *w, *x;
xue@1:   double* win;
xue@1:   int* hbi;
xue@1: 
xue@1:   //obtain FFT buffers win (window function), w (twiddle factors), x (data buffer),
xue@1:   //hbi (half-size bit-inversed integer table)
xue@1:   if (GetFFTBuffers)  //then use external buffers provided through GetFFTBuffers
xue@1:       GetFFTBuffers(Id, w, x, win, hbi, Parent);
xue@1:   else  //then use internal buffers
xue@1:   {
xue@1:     if (Wid!=fWid)
xue@1:     { //then update internal buffers to the new window size
xue@1:       free8(fw); free(fhbi);
xue@1:       fw=(cdouble*)malloc8(sizeof(cdouble)*Wid*1.5); SetTwiddleFactors(Wid, fw);
xue@1:       fx=&fw[Wid/2];
xue@11:       fhbi=CreateBitInvTable(Log2(Wid)-1);
xue@1:     }
xue@1:     if (Wid!=fWid || WinType!=fwt || WinParam!=fwdp)
xue@1:     { //then update internal window function to the new window type
xue@1: 			fwin=NewWindow8(WinType, Wid, 0, &WinParam);
xue@1:       fwt=WinType; fwdp=WinParam;
xue@1:     }
xue@1:     fWid=Wid;
xue@1: 
xue@1:     //pick up the internal buffers
xue@1: 		w=fw, x=fx, win=fwin, hbi=fhbi;
xue@1:   }
xue@1: 
xue@1:   //obtain the index of this frame in internal storage
xue@1:   if (Frame[fr]<0) {AddBuffer(1); Frame[fr]=FrCount; FrCount++;}
xue@1:   int realfr=Frame[fr];
xue@1: 
xue@1:   //obtain the pointer to this frame's phase spectrum
xue@1:   QSPEC_FORMAT *lph=useph?fPh[realfr]:NULL;
xue@1:   //ontain the pointer to this frame's complex spectrum
xue@1:   cmplx<QSPEC_FORMAT>* lX=usex?fSpec[realfr]:NULL;
xue@1:   //choose the buffer actually used for FFT - use lX if it is specified as complex double array
xue@1:   //because it saves unnecessary data copying operations
xue@1:   cdouble *lx=(usex && sizeof(QSPEC_FORMAT)==sizeof(double))?(cdouble*)lX:x;
xue@1: 
xue@1:   //Actual FFT
xue@1:   if (fr*Offst+Wid<=DataLength)
xue@1:     ::CalculateSpectrum(&((char*)Data)[fr*Offst*BytesPerSample], BytesPerSample, win, fA[realfr], lph, Wid, w, lx, hbi);
xue@1:   else
xue@1:     ::CalculateSpectrum(&((char*)Data)[fr*Offst*BytesPerSample], BytesPerSample, win, fA[realfr], lph, Wid, DataLength-fr*Offst, w, x, hbi);
xue@1: 
xue@1:   //optional data copy from x to lX
xue@1:   if (usex && lx==x) for (int i=0; i<Wid/2+1; i++) lX[i]=x[i];
xue@1: 
xue@1:   //tag this frame as computed and valid
xue@1:   Valid[fr]=1;
xue@1: }//CalculateSpectrum
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   method TQuickSpectrogram::FreeBuffers: discards all computed spectra and free all internal buffers.
xue@1:   This returns the TQuickSpectrogram to its initial state before any frame is accessed. After calling
xue@1:   FreeBuffers() all frames will be recomputed when they are accessed.
xue@1: */
xue@1: void TQuickSpectrogram::FreeBuffers()
xue@1: {
xue@1:   if (fA)
xue@1:   {
xue@1:     for (int i=0; i<BufferCount; i++) free(fA[i*BufferSize]);
xue@1:     FrCount=BufferCount=Capacity=0;
xue@1:     free(Frame); free(Valid);
xue@1:     free(fA);
xue@1:     Frame=Valid=0, fA=0;
xue@1:     if (useph) {free(fPh); fPh=0;}
xue@1:     if (usex) {free(fSpec); fSpec=0;}
xue@1:   }
xue@1: }//FreeBuffers
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   method TQuickSpectrogram::Invalidate: renders all frames that have overlap with interval [From, To],
xue@1:   measured in samples, as invalid. Invalid frames are recomputed when they are accessed again.
xue@1: 
xue@1:   In: [From, To]: an interval spectrogram over which needs to be updated.
xue@1: 
xue@1:   Returns the number of allocated frames affected, no matter if they were valid.
xue@1: */
xue@1: int TQuickSpectrogram::Invalidate(int From, int To)
xue@1: {
xue@1:   int result=0;
xue@1:   if (Frame)
xue@1:   {
xue@1:     int fr1=ceil((From-Wid+1.0)/Offst), fr2=floor(1.0*To/Offst);
xue@1:     if (fr1<0) fr1=0;
xue@1:     if (fr2>=Capacity) fr2=Capacity-1;
xue@1:     for (int fr=fr1; fr<=fr2; fr++) if (Frame[fr]>=0) Valid[fr]=false, result++;
xue@1:   }
xue@1:   return result;
xue@1: }//Invalidate
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   method TQuickSpectrogram::Ph: accesses phase spectrogram by frame
xue@1: 
xue@1:   In: fr: frame index, 0-based
xue@1: 
xue@1:   Returns pointer to phase spectrum of the fr'th frame, NULL if N/A
xue@1: */
xue@1: QSPEC_FORMAT* __fastcall TQuickSpectrogram::Ph(int fr)
xue@1: {
xue@1:   if (Capacity==0) SetFrCapacity((DataLength-Wid)/Offst+2);
xue@1:   if (fr<0 || fr>=Capacity) return NULL;
xue@1:   if (Frame[fr]<0 || !Valid[fr]) CalculateSpectrum(fr);
xue@1:   return fPh[Frame[fr]];
xue@1: }//Ph
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   method TQuickSpectrogram::SetFrCapacity: sets the capacity, i.e. the maximal number of frames handled
xue@1:   by this TQuickSpectrogram.
xue@1: 
xue@1:   In: AnFrCapacity: the new Capacity, in frames
xue@1: 
xue@1:   This method should not be called to set Capacity to a smaller value.
xue@1: */
xue@1: void TQuickSpectrogram::SetFrCapacity(int AnFrCapacity)
xue@1: {
xue@1:   //adjusting the size of index and validity arrays
xue@1:   Frame=(int*)realloc(Frame, sizeof(int)*AnFrCapacity);
xue@1:   Valid=(int*)realloc(Valid, sizeof(int)*AnFrCapacity);
xue@1: 
xue@1:   //
xue@1:   fA=(QSPEC_FORMAT**)realloc(fA, sizeof(QSPEC_FORMAT*)*AnFrCapacity);
xue@1:   if (usex) fSpec=(cmplx<QSPEC_FORMAT>**)realloc(fSpec, sizeof(cmplx<QSPEC_FORMAT>*)*AnFrCapacity);
xue@1:   if (useph) fPh=(QSPEC_FORMAT**)realloc(fPh, sizeof(QSPEC_FORMAT*)*AnFrCapacity);
xue@1:   if (AnFrCapacity>Capacity)
xue@1:   {
xue@1:     memset(&Frame[Capacity], 0xFF, sizeof(int)*(AnFrCapacity-Capacity));
xue@1:     memset(&Valid[Capacity], 0x00, sizeof(int)*(AnFrCapacity-Capacity));
xue@1: 
xue@1:     if (Capacity<BufferCount*BufferSize)
xue@1:     {
xue@1:       for (int fr=Capacity; fr<AnFrCapacity; fr++)
xue@1:       {
xue@1:         int bufferno=fr/BufferSize;
xue@1:         if (bufferno<BufferCount)
xue@1:         {
xue@1:           QSPEC_FORMAT* thisbuffer=fA[BufferSize*bufferno];
xue@1:           int lfr=fr%BufferSize, base=1, Dim=Wid/2+1;
xue@1:           fA[fr]=&thisbuffer[lfr*Dim];
xue@1:           if (usex) fSpec[fr]=(cmplx<QSPEC_FORMAT>*)(&thisbuffer[(BufferSize+lfr*2)*Dim]), base+=2;
xue@1:           if (useph) fPh[fr]=&thisbuffer[(BufferSize*base+lfr)*Dim];
xue@1:         }
xue@1:         else break;
xue@1:       }
xue@1:     }
xue@1:   }
xue@1:   Capacity=AnFrCapacity;
xue@1: }//SetFrCapacity
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   method TQuickSpectrogram::Ph: accesses complex spectrogram by frame
xue@1: 
xue@1:   In: fr: frame index, 0-based
xue@1: 
xue@1:   Returns pointer to complex spectrum of the fr'th frame, NULL if N/A
xue@1: */
xue@1: cmplx<QSPEC_FORMAT>* __fastcall TQuickSpectrogram::Spec(int fr)
xue@1: {
xue@1:   if (Capacity==0) SetFrCapacity((DataLength-Wid)/Offst+2);
xue@1:   if (fr<0 || fr>=Capacity) return NULL;
xue@1:   if (Frame[fr]<0 || !Valid[fr]) CalculateSpectrum(fr);
xue@1:   return fSpec[Frame[fr]];
xue@1: }//Spec
xue@1: 
xue@1: