Fixing bugs in SMSResidual.

8 files changed, 176 insertions, 113 deletions

diff --git a/sms/analysis.c b/sms/analysis.c
index a6b344d..95ddf60 100644
--- a/sms/analysis.c
+++ b/sms/analysis.c
@@ -331,91 +331,18 @@ int sms_findPartials(SMS_Data *pSmsData, SMS_AnalParams *pAnalParams)
 
 int sms_findResidual(int sizeSynthesis, sfloat* pSynthesis,
                      int sizeOriginal, sfloat* pOriginal,
-                     int sizeResidual, sfloat* pResidual,
-                     SMS_AnalParams *analParams)
+                     SMS_ResidualParams *residualParams)
 {
-    if(sizeResidual < sizeOriginal)
+    if(residualParams->residualSize < sizeOriginal)
     {
         sms_error("Residual signal length is smaller than the original signal length");
         return -1;
     }
 
-    /*sms_residual(sizeResidual, pSynthesis, pOriginal, pResidual);*/
+    sms_residual(residualParams->residualSize, pSynthesis, pOriginal, residualParams);
     return 0;
 }
 
-void sms_approxResidual(int sizeResidual, sfloat* pResidual,
-                        SMS_Data* pSmsData, SMS_SynthParams* pSynthParams)
-{
-    /* perform stochastic analysis after 1 frame of the     */
-    /* deterministic synthesis because it needs two frames  */
-//  if (pAnalParams->ppFrames[0]->iStatus != SMS_FRAME_EMPTY &&
-//      pAnalParams->ppFrames[0]->iStatus != SMS_FRAME_END)
-//  {
-//      int sizeResidual = pAnalParams->sizeHop * 2;
-//      int iSoundLoc = pAnalParams->ppFrames[0]->iFrameSample - pAnalParams->sizeHop;
-//  //          sfloat *pOriginal = &(pAnalParams->soundBuffer.pFBuffer[iSoundLoc - pAnalParams->soundBuffer.iMarker]);
-//      sfloat *pOriginal;
-//      sfloat *pFResidual;
-//
-//      static sfloat *pWindow;
-//      static int sizeWindowArray = 0;
-//
-//  //          int sizeData = MIN(pAnalParams->soundBuffer.sizeBuffer -
-//  //                            (iSoundLoc - pAnalParams->soundBuffer.iMarker),
-//  //                            sizeResidual);
-//      int sizeData = sizeResidual;
-//      if ((pFResidual = (sfloat *) calloc (sizeResidual, sizeof(sfloat))) == NULL)
-//      {
-//          sms_error("sms_analyze: error allocating memory for pFResidual");
-//          return -1;
-//      }
-//      if (sizeWindowArray != sizeData)
-//      {
-//          if(sizeWindowArray != 0) free(pWindow);
-//          if((pWindow = (sfloat *) calloc(sizeData, sizeof(sfloat))) == NULL)
-//          {
-//                  sms_error("sms_analyze: error allocating memory for pWindow");
-//                  return -1;
-//          }
-//          sms_getWindow( sizeData, pWindow, SMS_WIN_HAMMING);
-//          sms_scaleWindow( sizeData, pWindow);
-//          sizeWindowArray = sizeData;
-//      }
-//
-//      /* obtain residual sound from original and synthesized sounds.  accumulate the residual percentage.*/
-//      pAnalParams->fResidualAccumPerc += sms_residual(sizeData,
-//                                                      pAnalParams->synthBuffer.pFBuffer,
-//                                                      pOriginal,
-//                                                      pFResidual,
-//                                                      pWindow);
-//
-//      if (pAnalParams->iStochasticType == SMS_STOC_APPROX)
-//      {
-//          /* filter residual with a high pass filter (it solves some problems) */
-//          sms_filterHighPass (sizeData, pFResidual, pAnalParams->iSamplingRate);
-//
-//          /* approximate residual */
-//          sms_stocAnalysis (sizeData, pFResidual, pWindow, pSmsData);
-//      }
-//      else if  (pAnalParams->iStochasticType == SMS_STOC_IFFT)
-//      {
-//          int sizeMag = sms_power2(sizeData >> 1);
-//          sms_spectrum (sizeData, pFResidual, pWindow, sizeMag, pSmsData->pFStocCoeff,
-//                  pSmsData->pResPhase);
-//      }
-//
-//      /* get sharper transitions in deterministic representation */
-//  //              sms_scaleDet (pAnalParams->synthBuffer.pFBuffer, pOriginal,
-//  //                            pAnalParams->ppFrames[0]->deterministic.pFSinAmp,
-//  //                            pAnalParams, pSmsData->nTracks);
-//
-//      pAnalParams->ppFrames[0]->iStatus = SMS_FRAME_DONE;
-//
-//      free ((char *) pFResidual); /* \todo get rid of this free, manage memory the same as spectrum functions */
-//  }
-}
-
 int sms_analyze(int sizeWaveform, sfloat *pWaveform, SMS_Data *pSmsData, SMS_AnalParams *pAnalParams)
 {
     int iCurrentFrame = pAnalParams->iMaxDelayFrames - 1;  /* frame # of current frame */
diff --git a/sms/residual.c b/sms/residual.c
index 1c217ce..0c720f0 100644
--- a/sms/residual.c
+++ b/sms/residual.c
@@ -36,8 +36,6 @@
 int sms_residual(int sizeWindow, sfloat *pSynthesis, sfloat *pOriginal,
                  SMS_ResidualParams* residualParams)
 {
-    /*static sfloat fResidualMag = 0.;*/
-    /*static sfloat fOriginalMag = 0.;*/
     sfloat fScale = 1.;
     sfloat fCurrentResidualMag = 0.;
     sfloat fCurrentOriginalMag = 0.;
diff --git a/sms/sms.c b/sms/sms.c
index 144f3a6..4fee489 100644
--- a/sms/sms.c
+++ b/sms/sms.c
@@ -466,11 +466,23 @@ int sms_initSynth(SMS_SynthParams *pSynthParams)
  */
 void sms_initResidualParams(SMS_ResidualParams *residualParams)
 {
+    residualParams->samplingRate = 44100;
     residualParams->residualSize = 0;
     residualParams->residual = NULL;
     residualParams->residualWindow = NULL;
     residualParams->residualMag = 0.0;
     residualParams->originalMag = 0.0;
+    residualParams->nCoeffs = 128;
+    residualParams->stocCoeffs = NULL;
+    residualParams->sizeStocMagSpectrum = 0;
+    residualParams->stocMagSpectrum = NULL;
+    residualParams->approxEnvelope = NULL;
+    int i;
+    for(i = 0; i < SMS_MAX_SPEC; i++)
+    {
+        residualParams->fftBuffer[i] = 0.0;
+        residualParams->fftBuffer[i+SMS_MAX_SPEC] = 0.0;
+    }
 }
 
 /*! \brief initialize residual data structure
@@ -504,6 +516,29 @@ int sms_initResidual(SMS_ResidualParams *residualParams)
     sms_getWindow(residualParams->residualSize, residualParams->residualWindow, SMS_WIN_HAMMING);
     sms_scaleWindow(residualParams->residualSize, residualParams->residualWindow);
 
+    /* stochastic analysis */
+    residualParams->stocCoeffs = (sfloat *)calloc(residualParams->nCoeffs, sizeof(sfloat));
+    if(residualParams->stocCoeffs == NULL)
+    {
+        sms_error("Could not allocate memory for stochastic coefficients");
+        return -1;
+    }
+
+    residualParams->sizeStocMagSpectrum = sms_power2(residualParams->residualSize) >> 1;
+    residualParams->stocMagSpectrum = (sfloat *)calloc(residualParams->sizeStocMagSpectrum, sizeof(sfloat));
+    if(residualParams->stocMagSpectrum == NULL)
+    {
+        sms_error("Could not allocate memory for stochastic magnitude spectrum");
+        return -1;
+    }
+
+    residualParams->approxEnvelope = (sfloat *)calloc(residualParams->nCoeffs, sizeof(sfloat));
+    if(residualParams->approxEnvelope == NULL)
+    {
+        sms_error("Could not allocate memory for spectral approximation envelope");
+        return -1;
+    }
+
     return 0;
 }
 
@@ -520,9 +555,18 @@ void sms_freeResidual(SMS_ResidualParams *residualParams)
         free(residualParams->residual);
     if(residualParams->residualWindow)
         free(residualParams->residualWindow);
+    if(residualParams->stocCoeffs)
+        free(residualParams->stocCoeffs);
+    if(residualParams->stocMagSpectrum)
+        free(residualParams->stocMagSpectrum);
+    if(residualParams->approxEnvelope)
+        free(residualParams->approxEnvelope);
 
     residualParams->residual = NULL;
     residualParams->residualWindow = NULL;
+    residualParams->stocCoeffs = NULL;
+    residualParams->stocMagSpectrum = NULL;
+    residualParams->approxEnvelope = NULL;
 }
 
 /*! \brief free analysis data
diff --git a/sms/sms.h b/sms/sms.h
index 4726760..7e80ee4 100644
--- a/sms/sms.h
+++ b/sms/sms.h
@@ -195,11 +195,18 @@ typedef struct
  */
 typedef struct
 {
+    int samplingRate;
     int residualSize;
     sfloat *residual;
     sfloat *residualWindow;
     sfloat residualMag;
     sfloat originalMag;
+    int nCoeffs;
+    sfloat *stocCoeffs;
+    int sizeStocMagSpectrum;
+    sfloat *stocMagSpectrum;
+    sfloat *approxEnvelope;
+    sfloat fftBuffer[SMS_MAX_SPEC * 2];
 } SMS_ResidualParams;
 
 /*! \struct SMS_AnalParams
@@ -561,13 +568,13 @@ sfloat sms_dBToMag(sfloat x);
 void sms_arrayMagToDB(int sizeArray, sfloat *pArray);
 void sms_arrayDBToMag(int sizeArray, sfloat *pArray);
 void sms_setMagThresh(sfloat x);
-sfloat sms_rms ( int sizeArray, sfloat *pArray );
-sfloat sms_sine (sfloat fTheta);
-sfloat sms_sinc (sfloat fTheta);
-sfloat sms_random ( void );
+sfloat sms_rms(int sizeArray, sfloat *pArray);
+sfloat sms_sine(sfloat fTheta);
+sfloat sms_sinc(sfloat fTheta);
+sfloat sms_random(void);
 int sms_power2(int n);
-sfloat sms_scalarTempered( sfloat x);
-void sms_arrayScalarTempered( int sizeArray, sfloat *pArray);
+sfloat sms_scalarTempered(sfloat x);
+void sms_arrayScalarTempered(int sizeArray, sfloat *pArray);
 
 #ifndef MAX
 /*! \brief returns the maximum of a and b */
@@ -586,8 +593,8 @@ int sms_findPeaks(int sizeWaveform, sfloat *pWaveform, SMS_AnalParams *pAnalPara
 int sms_findPartials(SMS_Data *pSmsFrame, SMS_AnalParams *pAnalParams);
 int sms_findResidual(int sizeSynthesis, sfloat* pSynthesis,
                      int sizeOriginal, sfloat* pOriginal,
-                     int sizeResidual, sfloat* pResidual,
-                     SMS_AnalParams *analParams);
+                     SMS_ResidualParams *residualParams);
+void sms_approxResidual(SMS_ResidualParams *residualParams);
 int sms_analyze(int sizeWaveform, sfloat *pWaveform, SMS_Data *pSmsData, SMS_AnalParams *pAnalParams);
 void sms_analyzeFrame(int iCurrentFrame, SMS_AnalParams *pAnalParams, sfloat fRefFundamental);
 
diff --git a/sms/sms.i b/sms/sms.i
index a61d4b8..afde07a 100644
--- a/sms/sms.i
+++ b/sms/sms.i
@@ -500,6 +500,32 @@
     }
 }
 
+%extend SMS_ResidualParams
+{
+    void getResidual(int sizeArray, sfloat *pArray)
+    {
+        if(sizeArray < $self->residualSize)
+        {
+            sms_error("numpy array not big enough");
+            return;
+        }
+        int i;
+        for(i = 0; i < $self->residualSize; i++)
+            pArray[i] = $self->residual[i];
+    }
+    void getApprox(int sizeArray, sfloat *pArray)
+    {
+        if(sizeArray < $self->nCoeffs)
+        {
+            sms_error("numpy array not big enough");
+            return;
+        }
+        int i;
+        for(i = 0; i < $self->nCoeffs; i++)
+            pArray[i] = $self->approxEnvelope[i];
+    }
+}
+
 %extend SMS_ModifyParams
 {
     /* no need to return an error code, if sms_error is called, it will throw an exception in python */
diff --git a/sms/spectralApprox.c b/sms/spectralApprox.c
index c19afbe..82c94e3 100644
--- a/sms/spectralApprox.c
+++ b/sms/spectralApprox.c
@@ -28,14 +28,13 @@
  * interpolation. The output spectrum doesn't have to be the same size as 
  * the input one.
  *
- * \param pFSpec1     magnitude spectrum to approximate
- * \param sizeSpec1       size of input spectrum
- * \param sizeSpec1Used  size of the spectrum to use
- * \param pFSpec2     output envelope
- * \param sizeSpec2      size of output envelope
- * \param nCoefficients  number of coefficients to use in approximation
- * \return error code \see SMS_ERRORS (or -1 if the algorithm just messes up, 
- *         it will print an error of its own.
+ * \param pFSpec1       magnitude spectrum to approximate
+ * \param sizeSpec1     size of input spectrum
+ * \param sizeSpec1Used size of the spectrum to use
+ * \param pFSpec2       output envelope
+ * \param sizeSpec2     size of output envelope
+ * \param nCoefficients number of coefficients to use in approximation
+ * \return error code \see SMS_ERRORS
  */
 int sms_spectralApprox(sfloat *pFSpec1, int sizeSpec1, int sizeSpec1Used,
                        sfloat *pFSpec2, int sizeSpec2, int nCoefficients,
@@ -48,44 +47,49 @@ int sms_spectralApprox(sfloat *pFSpec1, int sizeSpec1, int sizeSpec1Used,
     /* when number of coefficients is smaller than 2 do not approximate */
     if(nCoefficients < 2)
     {
-        for (i = 0; i < sizeSpec2; i++)
+        for(i = 0; i < sizeSpec2; i++)
             pFSpec2[i] = 1;
         return SMS_OK;
     }
 
     /* calculate the hop size */
-    if(sizeSpec1 != sizeSpec1Used)
-        fHopSize = (sfloat) sizeSpec1Used / nCoefficients;
-    else //why is this here, would be the same as sizeSpec1Used / nCoefficients
-        fHopSize = (sfloat) sizeSpec1 / nCoefficients;
-
+    /*if(sizeSpec1 != sizeSpec1Used)*/
+        /*fHopSize = (sfloat) sizeSpec1Used / nCoefficients;*/
+    /*else //why is this here, would be the same as sizeSpec1Used / nCoefficients*/
+        /*fHopSize = (sfloat) sizeSpec1 / nCoefficients;*/
     if(nCoefficients > sizeSpec1)
         nCoefficients = sizeSpec1;
 
-    fHopSize = (sfloat) sizeSpec1Used / nCoefficients;
+    fHopSize = (sfloat)sizeSpec1Used / nCoefficients;
 
     /* approximate by linear interpolation */
-    if (fHopSize > 1)
+    if(fHopSize > 1)
     {
         iFirstGood = 0;
         for(i = 0; i < nCoefficients; i++)
         {
             iLastSample = fLastLocation = fCurrentLoc + fHopSize;
-            iLastSample = MIN (sizeSpec1-1, iLastSample);
+            iLastSample = MIN(sizeSpec1-1, iLastSample);
             if(iLastSample < sizeSpec1-1)
+            {
                 fRight = pFSpec1[iLastSample] +
-                    (pFSpec1[iLastSample+1] - pFSpec1[iLastSample]) * 
-                    (fLastLocation - iLastSample);
+                         (pFSpec1[iLastSample+1] - pFSpec1[iLastSample]) * 
+                         (fLastLocation - iLastSample);
+            }
             else
+            {
                 fRight = pFSpec1[iLastSample];
+            }
+
             fValue = 0;
             for(j = iFirstGood; j <= iLastSample; j++)
                 fValue = MAX (fValue, pFSpec1[j]);
-            fValue = MAX (fValue, MAX (fRight, fLeft));
+            fValue = MAX(fValue, MAX (fRight, fLeft));
             envelope[i] = fValue;
+
             fLeft = fRight;
             fCurrentLoc = fLastLocation;
-            iFirstGood = (int) (1+ fCurrentLoc);
+            iFirstGood = (int)(1+ fCurrentLoc);
         }
     }
     else if(fHopSize == 1)
@@ -95,7 +99,7 @@ int sms_spectralApprox(sfloat *pFSpec1, int sizeSpec1, int sizeSpec1Used,
     }
     else
     {
-        sms_error("SpectralApprox: sizeSpec1 has too many nCoefficients\n"); /* \todo need to increase the frequency? */
+        sms_error("SpectralApprox: sizeSpec1 has too many nCoefficients"); /* \todo need to increase the frequency? */
         return -1;
     }
 
diff --git a/sms/stocAnalysis.c b/sms/stocAnalysis.c
index d073ee7..802e41c 100644
--- a/sms/stocAnalysis.c
+++ b/sms/stocAnalysis.c
@@ -23,8 +23,6 @@
  */
 #include "sms.h"
 
-#define ENV_THRESHOLD     .01 /* \todo was this value for type shorts?? */
-
 /*! \brief main function for the stochastic analysis
  * \param sizeWindow   size of buffer
  * \param pResidual    pointer to residual signal
@@ -38,7 +36,6 @@ int sms_stocAnalysis(int sizeWindow, sfloat *pResidual, sfloat *pWindow,
 {
     int i;
     sfloat fMag = 0.0;
-    sfloat fStocNorm;
 
     sms_spectrumMag(sizeWindow, pResidual, pWindow, pAnalParams->sizeStocMagSpectrum,
                     pAnalParams->stocMagSpectrum, pAnalParams->fftBuffer);
@@ -53,8 +50,6 @@ int sms_stocAnalysis(int sizeWindow, sfloat *pResidual, sfloat *pWindow,
         fMag += (pAnalParams->stocMagSpectrum[i] * pAnalParams->stocMagSpectrum[i]);
 
     *pSmsData->pFStocGain = fMag / pAnalParams->sizeStocMagSpectrum;
-    fStocNorm = 1. / *pSmsData->pFStocGain;
-
     return 0;
 }
 
diff --git a/sms/synthesis.c b/sms/synthesis.c
index 8a679ca..e9b0967 100644
--- a/sms/synthesis.c
+++ b/sms/synthesis.c
@@ -146,6 +146,68 @@ static int StocSynthApprox(SMS_Data *pSmsData, SMS_SynthParams *pSynthParams)
     return 1;
 }
 
+/*! \brief synthesizes one frame of the residual signal
+ *
+ * \param residualParams Parameters and memory for residual synthesis
+ */
+void sms_approxResidual(SMS_ResidualParams *residualParams)
+{
+    /* filter residual with a high pass filter */
+    sms_filterHighPass(residualParams->residualSize, 
+                       residualParams->residual, 
+                       residualParams->samplingRate);
+
+    sms_spectrumMag(residualParams->residualSize, 
+                    residualParams->residual, 
+                    residualParams->residualWindow, 
+                    residualParams->sizeStocMagSpectrum,
+                    residualParams->stocMagSpectrum, 
+                    residualParams->fftBuffer);
+
+    sms_spectralApprox(residualParams->stocMagSpectrum, 
+                       residualParams->sizeStocMagSpectrum,
+                       residualParams->sizeStocMagSpectrum, 
+                       residualParams->stocCoeffs,
+                       residualParams->nCoeffs, 
+                       residualParams->nCoeffs,
+                       residualParams->approxEnvelope);
+
+    /* get energy of spectrum  */
+    int i;
+    sfloat fMag = 0.0;
+    /*for(i = 0; i < pAnalParams->sizeStocMagSpectrum; i++)*/
+    /*    fMag += (pAnalParams->stocMagSpectrum[i] * pAnalParams->stocMagSpectrum[i]);*/
+
+    /**pSmsData->pFStocGain = fMag / pAnalParams->sizeStocMagSpectrum;*/
+
+    /*int i, sizeSpec1Used;*/
+    /*int sizeSpec1 = pSmsData->nCoeff;*/
+    /*int sizeSpec2 = pSynthParams->sizeHop;*/
+    /*int sizeFft = pSynthParams->sizeHop << 1; [> 50% overlap, so sizeFft is 2x sizeHop <]*/
+    /*sfloat fStocGain;*/
+
+    /*[> if no gain or no coefficients return  <]*/
+    /*if (*(pSmsData->pFStocGain) <= 0)*/
+    /*    return 0;*/
+
+    /*sizeSpec1Used = sizeSpec1 * pSynthParams->iSamplingRate / pSynthParams->iOriginalSRate;*/
+
+    /*[> sizeSpec1Used cannot be more than what is available  \todo check by graph <]*/
+    /*if(sizeSpec1Used  > sizeSpec1) sizeSpec1Used = sizeSpec1;*/
+
+    /*sms_spectralApprox(pSmsData->pFStocCoeff, sizeSpec1, sizeSpec1Used,*/
+    /*                   pSynthParams->pMagBuff, sizeSpec2, sizeSpec1Used,*/
+    /*                   pSynthParams->approxEnvelope);*/
+
+    /*[> generate random phases <]*/
+    /*for(i = 0; i < sizeSpec2; i++)*/
+    /*    pSynthParams->pPhaseBuff[i] =  TWO_PI * sms_random();*/
+
+    /*sms_invQuickSpectrumW(pSynthParams->pMagBuff, pSynthParams->pPhaseBuff,*/
+    /*                      sizeFft, pSynthParams->pSynthBuff, sizeFft,*/
+    /*                      pSynthParams->pFStocWindow, pSynthParams->pSpectra);*/
+}
+
 /*! \brief  synthesizes one frame of SMS data
  *
  * \param pSmsData      input SMS data