Fixed SMSResidual. Also removed blank frames produced by sms at the beginning of analysis which was causing synthesised audio to be out of sync with the original.

7 files changed, 210 insertions, 131 deletions

diff --git a/sms/analysis.c b/sms/analysis.c
index ad8420e..89d4b4a 100644
--- a/sms/analysis.c
+++ b/sms/analysis.c
@@ -210,7 +210,6 @@ int sms_findPeaks(int sizeWaveform, sfloat *pWaveform, SMS_AnalParams *pAnalPara
                 pSpectralPeaks->pSpectralPeaks[i].fPhase = 0.0;
             }
         }
-        /*printf("\n");*/
         return pSpectralPeaks->nPeaks;
     }
     else
@@ -328,13 +327,16 @@ int sms_findResidual(int sizeSynthesis, sfloat* pSynthesis,
                      int sizeOriginal, sfloat* pOriginal,
                      SMS_ResidualParams *residualParams)
 {
-    if(residualParams->residualSize < sizeOriginal)
+    if(residualParams->hopSize < sizeOriginal)
     {
         sms_error("Residual signal length is smaller than the original signal length");
         return -1;
     }
 
-    sms_residual(residualParams->residualSize, pSynthesis, pOriginal, residualParams);
+    sms_residual(residualParams->hopSize, pSynthesis, pOriginal, residualParams);
+    sms_filterHighPass(residualParams->hopSize,
+                       residualParams->residual,
+                       residualParams->samplingRate);
     return 0;
 }
 
@@ -467,8 +469,9 @@ int sms_analyze(int sizeWaveform, sfloat *pWaveform, SMS_Data *pSmsData, SMS_Ana
             else if(sizeData < pAnalParams->residualParams.residualSize)
             {
                 /* should only happen if we're at the end of a sound, unless hop size changes */
-                sms_getWindow(sizeData, pAnalParams->residualParams.residualWindow, SMS_WIN_HAMMING);
-                sms_scaleWindow(sizeData, pAnalParams->residualParams.residualWindow);
+                /* TODO: should the window type be set to pAnalParams->iWindowType? */
+                sms_getWindow(sizeData, pAnalParams->residualParams.fftWindow, SMS_WIN_HAMMING);
+                sms_scaleWindow(sizeData, pAnalParams->residualParams.fftWindow);
             }
 
             /* obtain residual sound from original and synthesized sounds.  accumulate the residual percentage.*/
@@ -483,13 +486,13 @@ int sms_analyze(int sizeWaveform, sfloat *pWaveform, SMS_Data *pSmsData, SMS_Ana
                 sms_filterHighPass(sizeData, pAnalParams->residualParams.residual, pAnalParams->iSamplingRate);
 
                 /* approximate residual */
-                sms_stocAnalysis(sizeData, pAnalParams->residualParams.residual, pAnalParams->residualParams.residualWindow, 
+                sms_stocAnalysis(sizeData, pAnalParams->residualParams.residual, pAnalParams->residualParams.fftWindow, 
                                  pSmsData, pAnalParams);
             }
             else if(pAnalParams->iStochasticType == SMS_STOC_IFFT)
             {
                 int sizeMag = sms_power2(sizeData >> 1);
-                sms_spectrum(sizeData, pAnalParams->residualParams.residual, pAnalParams->residualParams.residualWindow, 
+                sms_spectrum(sizeData, pAnalParams->residualParams.residual, pAnalParams->residualParams.fftWindow, 
                              sizeMag, pSmsData->pFStocCoeff, pSmsData->pResPhase,
                              pAnalParams->fftBuffer);
             }
diff --git a/sms/sineSynth.c b/sms/sineSynth.c
index c91bf65..9882062 100644
--- a/sms/sineSynth.c
+++ b/sms/sineSynth.c
@@ -45,18 +45,18 @@ static void SinePhaSynth(sfloat fFreq, sfloat fMag, sfloat fPhase,
     sfloat fAlpha, fBeta, fTmp1, fTmp2;
 
     /* if no mag in last frame copy freq from current and make phase */
-    if (pLastFrame->pFSinAmp[iTrack] <= 0)
+    if(pLastFrame->pFSinAmp[iTrack] <= 0)
     {
         pLastFrame->pFSinFreq[iTrack] = fFreq;
         fTmp = fPhase - (fFreq * sizeBuffer);
         pLastFrame->pFSinPha[iTrack] = fTmp - floor(fTmp / TWO_PI) * TWO_PI;
     }
     /* and the other way */
-    else if (fMag <= 0)
+    else if(fMag <= 0)
     {
         fFreq = pLastFrame->pFSinFreq[iTrack];
         fTmp = pLastFrame->pFSinPha[iTrack] + 
-            (pLastFrame->pFSinFreq[iTrack] * sizeBuffer);
+               (pLastFrame->pFSinFreq[iTrack] * sizeBuffer);
         fPhase = fTmp - floor(fTmp / TWO_PI) * TWO_PI;
     }
 
@@ -67,27 +67,27 @@ static void SinePhaSynth(sfloat fFreq, sfloat fMag, sfloat fPhase,
     /* create instantaneous phase from freq. and phase values */
     fTmp1 = fFreq - pLastFrame->pFSinFreq[iTrack];
     fTmp2 = ((pLastFrame->pFSinPha[iTrack] + 
-                pLastFrame->pFSinFreq[iTrack] * sizeBuffer - fPhase) +
-            fTmp1 * sizeBuffer / 2.0) / TWO_PI;
-    iM = (int) (fTmp2 + .5);
+              pLastFrame->pFSinFreq[iTrack] * sizeBuffer - fPhase) +
+             fTmp1 * sizeBuffer / 2.0) / TWO_PI;
+    iM = (int)(fTmp2 + .5);
     fTmp2 = fPhase - pLastFrame->pFSinPha[iTrack] - 
-        pLastFrame->pFSinFreq[iTrack] * sizeBuffer +
-        TWO_PI * iM;
+            pLastFrame->pFSinFreq[iTrack] * sizeBuffer + TWO_PI * iM;
     fAlpha = (3.0 / (sfloat)(sizeBuffer * sizeBuffer)) * 
-        fTmp2 - fTmp1 / sizeBuffer;
+             fTmp2 - fTmp1 / sizeBuffer;
     fBeta = (-2.0 / ((sfloat) (sizeBuffer * sizeBuffer * sizeBuffer))) * 
-        fTmp2 + fTmp1 / ((sfloat) (sizeBuffer * sizeBuffer));
+            fTmp2 + fTmp1 / ((sfloat) (sizeBuffer * sizeBuffer));
 
     for(i=0; i<sizeBuffer; i++)
     {
         fInstMag += fMagIncr;
         fInstPhase = pLastFrame->pFSinPha[iTrack] + 
-            pLastFrame->pFSinFreq[iTrack] * i + 
-            fAlpha * i * i + fBeta * i * i * i;
+                     pLastFrame->pFSinFreq[iTrack] * i + 
+                     fAlpha * i * i + fBeta * i * i * i;
 
-        /*     pFWaveform[i] += sms_dBToMag(fInstMag) * sms_sine(fInstPhase + PI_2); */
+        /*pFWaveform[i] += sms_dBToMag(fInstMag) * sms_sine(fInstPhase + PI_2);*/
         pFWaveform[i] += sms_dBToMag(fInstMag) * sinf(fInstPhase + PI_2);
     }
+
     /* save current values into buffer */
     pLastFrame->pFSinFreq[iTrack] = fFreq;
     pLastFrame->pFSinAmp[iTrack] = fMag;
@@ -110,14 +110,13 @@ static void SineSynth(sfloat fFreq, sfloat fMag, SMS_Data *pLastFrame,
     int i;
 
     /* if no mag in last frame copy freq from current */
-    if (pLastFrame->pFSinAmp[iTrack] <= 0)
+    if(pLastFrame->pFSinAmp[iTrack] <= 0)
     {
         pLastFrame->pFSinFreq[iTrack] = fFreq;
-        pLastFrame->pFSinPha[iTrack] = 
-            TWO_PI * sms_random();
+        pLastFrame->pFSinPha[iTrack] = TWO_PI * sms_random();
     }
     /* and the other way */
-    else if (fMag <= 0)
+    else if(fMag <= 0)
         fFreq = pLastFrame->pFSinFreq[iTrack];
 
     /* calculate the instantaneous amplitude */
@@ -129,20 +128,18 @@ static void SineSynth(sfloat fFreq, sfloat fMag, SMS_Data *pLastFrame,
     fInstPhase = pLastFrame->pFSinPha[iTrack];
 
     /* generate all the samples */    
-    for (i = 0; i < sizeBuffer; i++)
+    for(i = 0; i < sizeBuffer; i++)
     {
         fInstMag += fMagIncr;
         fInstFreq += fFreqIncr;
         fInstPhase += fInstFreq;
-
         pFBuffer[i] += sms_dBToMag(fInstMag) * sms_sine(fInstPhase);
     }
 
     /* save current values into last values */
     pLastFrame->pFSinFreq[iTrack] = fFreq;
     pLastFrame->pFSinAmp[iTrack] = fMag;
-    pLastFrame->pFSinPha[iTrack] = fInstPhase - 
-        floor(fInstPhase / TWO_PI) * TWO_PI;
+    pLastFrame->pFSinPha[iTrack] = fInstPhase - floor(fInstPhase / TWO_PI) * TWO_PI;
 }
 
 /*! \brief generate all the sinusoids for a given frame
@@ -191,4 +188,3 @@ void sms_sineSynthFrame(SMS_Data *pSmsData, sfloat *pFBuffer,
         }
     }
 }     
-
diff --git a/sms/sms.c b/sms/sms.c
index 27b6724..8501d4d 100644
--- a/sms/sms.c
+++ b/sms/sms.c
@@ -319,7 +319,7 @@ int sms_initAnalysis(SMS_AnalParams *pAnalParams)
     }
 
     /* memory for residual */
-    pAnalParams->residualParams.residualSize = pAnalParams->sizeHop * 2;
+    pAnalParams->residualParams.hopSize = pAnalParams->sizeHop;
     sms_initResidual(&pAnalParams->residualParams);
 
     /* memory for guide states */
@@ -467,15 +467,20 @@ int sms_initSynth(SMS_SynthParams *pSynthParams)
 void sms_initResidualParams(SMS_ResidualParams *residualParams)
 {
     residualParams->samplingRate = 44100;
+    residualParams->hopSize = 256;
     residualParams->residualSize = 0;
     residualParams->residual = NULL;
-    residualParams->residualWindow = NULL;
+    residualParams->fftWindow = NULL;
+    residualParams->ifftWindow = NULL;
+    residualParams->windowScale = 0.0;
     residualParams->residualMag = 0.0;
     residualParams->originalMag = 0.0;
     residualParams->nCoeffs = 128;
     residualParams->stocCoeffs = NULL;
     residualParams->sizeStocMagSpectrum = 0;
     residualParams->stocMagSpectrum = NULL;
+    residualParams->stocPhaseSpectrum = NULL;
+    residualParams->approx = NULL;
     residualParams->approxEnvelope = NULL;
     int i;
     for(i = 0; i < SMS_MAX_SPEC; i++)
@@ -492,13 +497,14 @@ void sms_initResidualParams(SMS_ResidualParams *residualParams)
  */
 int sms_initResidual(SMS_ResidualParams *residualParams)
 {
-    if(residualParams->residualSize <= 0)
+    if(residualParams->hopSize <= 0)
     {
-        sms_error("Residual size must be a positive integer");
+        sms_error("Residual hop size must be a positive integer");
         return -1;
     }
 
     /* residual signal */
+    residualParams->residualSize = residualParams->hopSize * 2;
     residualParams->residual = (sfloat *)calloc(residualParams->residualSize, sizeof(sfloat));
     if(residualParams->residual == NULL)
     {
@@ -506,15 +512,34 @@ int sms_initResidual(SMS_ResidualParams *residualParams)
         return -1;
     }
 
-    /* residual window */
-    residualParams->residualWindow = (sfloat *)calloc(residualParams->residualSize, sizeof(sfloat));
-    if(residualParams->residualWindow == NULL)
+    /* residual fft/ifft windows */
+    residualParams->fftWindow = (sfloat *)calloc(residualParams->residualSize, sizeof(sfloat));
+    if(residualParams->fftWindow == NULL)
     {
-        sms_error("Could not allocate memory for residualWindow");
+        sms_error("Could not allocate memory for residual FFT window");
         return -1;
     }
-    sms_getWindow(residualParams->residualSize, residualParams->residualWindow, SMS_WIN_HAMMING);
-    sms_scaleWindow(residualParams->residualSize, residualParams->residualWindow);
+    sms_getWindow(residualParams->residualSize, residualParams->fftWindow, SMS_WIN_BH_70);
+    sms_scaleWindow(residualParams->residualSize, residualParams->fftWindow);
+
+    residualParams->ifftWindow = (sfloat *)calloc(residualParams->residualSize, sizeof(sfloat));
+    if(residualParams->ifftWindow == NULL)
+    {
+        sms_error("Could not allocate memory for residual IFFT window");
+        return -1;
+    }
+    sms_getWindow(residualParams->residualSize, residualParams->ifftWindow, SMS_WIN_HANNING);
+    /* compute IFFT window scaling:
+     * windows per hop = hop size / window size = 0.5
+     * overlap = 50% => 1 window total in each hop/frame
+     * => windowScale = window size / sum(window samples) = 1.85
+     *    for a 1024 sized hamming window
+     */
+    int i;
+    sfloat sum = 0.0;
+    for(i = 0; i < residualParams->residualSize; i++)
+        sum += residualParams->ifftWindow[i];
+    residualParams->windowScale = (sfloat)residualParams->residualSize / sum;
 
     /* stochastic analysis */
     residualParams->stocCoeffs = (sfloat *)calloc(residualParams->nCoeffs, sizeof(sfloat));
@@ -531,7 +556,19 @@ int sms_initResidual(SMS_ResidualParams *residualParams)
         sms_error("Could not allocate memory for stochastic magnitude spectrum");
         return -1;
     }
+    residualParams->stocPhaseSpectrum = (sfloat *)calloc(residualParams->sizeStocMagSpectrum, sizeof(sfloat));
+    if(residualParams->stocPhaseSpectrum == NULL)
+    {
+        sms_error("Could not allocate memory for stochastic magnitude spectrum");
+        return -1;
+    }
 
+    residualParams->approx = (sfloat *)calloc(residualParams->residualSize, sizeof(sfloat));
+    if(residualParams->approx == NULL)
+    {
+        sms_error("Could not allocate memory for spectral approximation");
+        return -1;
+    }
     residualParams->approxEnvelope = (sfloat *)calloc(residualParams->nCoeffs, sizeof(sfloat));
     if(residualParams->approxEnvelope == NULL)
     {
@@ -553,19 +590,28 @@ void sms_freeResidual(SMS_ResidualParams *residualParams)
 {
     if(residualParams->residual)
         free(residualParams->residual);
-    if(residualParams->residualWindow)
-        free(residualParams->residualWindow);
+    if(residualParams->fftWindow)
+        free(residualParams->fftWindow);
+    if(residualParams->ifftWindow)
+        free(residualParams->ifftWindow);
     if(residualParams->stocCoeffs)
         free(residualParams->stocCoeffs);
     if(residualParams->stocMagSpectrum)
         free(residualParams->stocMagSpectrum);
+    if(residualParams->stocPhaseSpectrum)
+        free(residualParams->stocPhaseSpectrum);
+    if(residualParams->approx)
+        free(residualParams->approx);
     if(residualParams->approxEnvelope)
         free(residualParams->approxEnvelope);
 
     residualParams->residual = NULL;
-    residualParams->residualWindow = NULL;
+    residualParams->fftWindow = NULL;
+    residualParams->ifftWindow = NULL;
     residualParams->stocCoeffs = NULL;
     residualParams->stocMagSpectrum = NULL;
+    residualParams->stocPhaseSpectrum = NULL;
+    residualParams->approx = NULL;
     residualParams->approxEnvelope = NULL;
 }
 
diff --git a/sms/sms.h b/sms/sms.h
index 7e80ee4..0613ab3 100644
--- a/sms/sms.h
+++ b/sms/sms.h
@@ -196,15 +196,20 @@ typedef struct
 typedef struct
 {
     int samplingRate;
+    int hopSize;
     int residualSize;
     sfloat *residual;
-    sfloat *residualWindow;
+    sfloat *fftWindow;
+    sfloat *ifftWindow;
+    sfloat windowScale;
     sfloat residualMag;
     sfloat originalMag;
     int nCoeffs;
     sfloat *stocCoeffs;
     int sizeStocMagSpectrum;
     sfloat *stocMagSpectrum;
+    sfloat *stocPhaseSpectrum;
+    sfloat *approx;
     sfloat *approxEnvelope;
     sfloat fftBuffer[SMS_MAX_SPEC * 2];
 } SMS_ResidualParams;
@@ -589,13 +594,17 @@ void sms_arrayScalarTempered(int sizeArray, sfloat *pArray);
 /* function declarations */ 
 void sms_setPeaks(SMS_AnalParams *pAnalParams, int numamps, sfloat* amps,
                   int numfreqs, sfloat* freqs, int numphases, sfloat* phases);
-int sms_findPeaks(int sizeWaveform, sfloat *pWaveform, SMS_AnalParams *pAnalParams, SMS_SpectralPeaks *pSpectralPeaks);
+int sms_findPeaks(int sizeWaveform, sfloat *pWaveform, 
+                  SMS_AnalParams *pAnalParams, SMS_SpectralPeaks *pSpectralPeaks);
 int sms_findPartials(SMS_Data *pSmsFrame, SMS_AnalParams *pAnalParams);
 int sms_findResidual(int sizeSynthesis, sfloat* pSynthesis,
                      int sizeOriginal, sfloat* pOriginal,
                      SMS_ResidualParams *residualParams);
-void sms_approxResidual(SMS_ResidualParams *residualParams);
-int sms_analyze(int sizeWaveform, sfloat *pWaveform, SMS_Data *pSmsData, SMS_AnalParams *pAnalParams);
+void sms_approxResidual(int sizeResidual, sfloat* residual,
+                        int sizeApprox, sfloat* approx, 
+                        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);
 
 int sms_init();  
@@ -615,6 +624,8 @@ void sms_getWindow(int sizeWindow, sfloat *pWindow, int iWindowType);
 void sms_scaleWindow(int sizeWindow, sfloat *pWindow);
 int sms_spectrum(int sizeWindow, sfloat *pWaveform, sfloat *pWindow, int sizeMag, 
                  sfloat *pMag, sfloat *pPhase, sfloat *pFftBuffer);
+int sms_spectrumW(int sizeWindow, sfloat *pWaveform, sfloat *pWindow, int sizeMag, 
+                  sfloat *pMag, sfloat *pPhase, sfloat *pFftBuffer);
 int sms_invSpectrum(int sizeWaveform, sfloat *pWaveform, sfloat *pWindow ,
                     int sizeMag, sfloat *pMag, sfloat *pPhase, sfloat *pFftBuffer);
 /* \todo remove this once invSpectrum is completely implemented */
diff --git a/sms/sms.i b/sms/sms.i
index afde07a..0b3ee5a 100644
--- a/sms/sms.i
+++ b/sms/sms.i
@@ -31,7 +31,6 @@
 %apply(int DIM1, double* INPLACE_ARRAY1) {(int sizeAmp, double* pAmp)};
 %apply(int DIM1, double* INPLACE_ARRAY1) {(int sizeMag, double* pMag)};
 %apply(int DIM1, double* INPLACE_ARRAY1) {(int sizePhase, double* pPhase)};
-%apply(int DIM1, double* INPLACE_ARRAY1) {(int sizeRes, double* pRes)};
 %apply(int DIM1, double* INPLACE_ARRAY1) {(int sizeCepstrum, double* pCepstrum)};
 %apply(int DIM1, double* INPLACE_ARRAY1) {(int sizeEnv, double* pEnv)};
 %apply(int DIM1, double* INPLACE_ARRAY1) {(int sizeTrack, double* pTrack)};
@@ -39,6 +38,11 @@
 %apply(int DIM1, double* IN_ARRAY1) {(int sizeInArray, double* pInArray)};
 %apply(int DIM1, double* INPLACE_ARRAY1) {(int sizeOutArray, double* pOutArray)};
 %apply(int DIM1, double* INPLACE_ARRAY1) {(int sizeHop, double* pSynthesis)};
+%apply(int DIM1, double* INPLACE_ARRAY1)
+{
+    (int sizeResidual, double* residual),
+    (int sizeApprox, double* approx)
+}
 %apply(int DIM1, double* IN_ARRAY1)
 {
     (int numamps, double* amps),
@@ -504,26 +508,15 @@
 {
     void getResidual(int sizeArray, sfloat *pArray)
     {
-        if(sizeArray < $self->residualSize)
+        if(sizeArray < $self->hopSize)
         {
             sms_error("numpy array not big enough");
             return;
         }
         int i;
-        for(i = 0; i < $self->residualSize; i++)
+        for(i = 0; i < $self->hopSize; 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
diff --git a/sms/spectrum.c b/sms/spectrum.c
index 666b42d..1b8e053 100644
--- a/sms/spectrum.c
+++ b/sms/spectrum.c
@@ -37,7 +37,6 @@ int sms_spectrum(int sizeWindow, sfloat *pWaveform, sfloat *pWindow, int sizeMag
                  sfloat *pMag, sfloat *pPhase, sfloat *pFftBuffer)
 {
     int i, it2;
-    int err = 0;
     sfloat fReal, fImag;
     
     int sizeFft = sizeMag << 1;
@@ -55,6 +54,38 @@ int sms_spectrum(int sizeWindow, sfloat *pWaveform, sfloat *pWindow, int sizeMag
         fImag = pFftBuffer[it2 + 1]; /*even numbers 2->N+2 */
         pMag[i] = sqrt(fReal * fReal + fImag * fImag);
         pPhase[i] = atan2(-fImag, fReal); /* \todo why is fImag negated? */
+        /*pPhase[i] = atan2(fImag, fReal);*/
+    }
+
+    return sizeFft;
+}
+
+/* sms_spectrum, but without zero-phase windowing, and with phase calculated
+ * according by arctan(imag/real) instead of arctan2(-imag/real)
+ */
+int sms_spectrumW(int sizeWindow, sfloat *pWaveform, sfloat *pWindow, int sizeMag, 
+                  sfloat *pMag, sfloat *pPhase, sfloat *pFftBuffer)
+{
+    int i, it2;
+    sfloat fReal, fImag;
+    
+    int sizeFft = sizeMag << 1;
+    memset(pFftBuffer, 0, sizeFft * sizeof(sfloat));
+
+    /* apply window to waveform */
+    for(i = 0; i < sizeWindow; i++)
+        pFftBuffer[i] = pWaveform[i] * pWindow[i];
+
+    sms_fft(sizeFft, pFftBuffer);
+
+    /* convert from rectangular to polar coordinates */
+    for(i = 0; i < sizeMag; i++)
+    { 
+        it2 = i << 1; //even numbers 0-N
+        fReal = pFftBuffer[it2]; /*odd numbers 1->N+1 */
+        fImag = pFftBuffer[it2 + 1]; /*even numbers 2->N+2 */
+        pMag[i] = sqrt(fReal * fReal + fImag * fImag);
+        pPhase[i] = atan2(fImag, fReal);
     }
 
     return sizeFft;
@@ -81,8 +112,8 @@ int sms_spectrumMag(int sizeWindow, sfloat *pWaveform, sfloat *pWindow,
     sfloat fReal, fImag;
 
     /* apply window to waveform, zero the rest of the array */
-    for (i = 0; i < sizeWindow; i++)
-        pFftBuffer[i] =  pWindow[i] * pWaveform[i];
+    for(i = 0; i < sizeWindow; i++)
+        pFftBuffer[i] =  pWaveform[i] * pWindow[i];
     for(i = sizeWindow; i < sizeFft; i++)
         pFftBuffer[i]  = 0.;
 
@@ -90,7 +121,7 @@ int sms_spectrumMag(int sizeWindow, sfloat *pWaveform, sfloat *pWindow,
     sms_fft(sizeFft, pFftBuffer); 
 
     /* convert from rectangular to polar coordinates */
-    for (i=0; i<sizeMag; i++)
+    for(i = 0; i < sizeMag; i++)
     {
         it2 = i << 1;
         fReal = pFftBuffer[it2];
@@ -123,7 +154,7 @@ int sms_invSpectrum(int sizeWaveform, sfloat *pWaveform, sfloat *pWindow,
     sms_PolarToRect(sizeMag, pFftBuffer, pMag, pPhase);
     sms_ifft(sizeFft, pFftBuffer); 
 
-    /* assume the output array has been taken care off */
+    /* assume that the output array does not need to be cleared */
     /* before, this was multiplied by .5, why? */
     for(i = 0; i < sizeWaveform; i++)
         //pWaveform[i] +=  pFftBuffer[i] * pWindow[i];
@@ -145,25 +176,21 @@ int sms_invQuickSpectrumW(sfloat *pFMagSpectrum, sfloat *pFPhaseSpectrum,
                           int sizeFft, sfloat *pFWaveform, int sizeWave,
                           sfloat *pFWindow, sfloat* pFftBuffer)
 {
-    int sizeMag = sizeFft >> 1, i, it2;
-    sfloat fPower;
+    int i, it2;
+    int sizeMag = sizeFft >> 1;
 
-    /* convert from polar coordinates to rectangular  */
-    for(i = 0; i<sizeMag; i++)
+    /* convert from polar coordinates to rectangular */
+    for(i = 0; i < sizeMag; i++)
     {
         it2 = i << 1;
-        fPower = pFMagSpectrum[i];
-        pFftBuffer[it2] =  fPower * cos (pFPhaseSpectrum[i]);
-        pFftBuffer[it2+1] = fPower * sin (pFPhaseSpectrum[i]);
+        pFftBuffer[it2] = pFMagSpectrum[i] * cos(pFPhaseSpectrum[i]);
+        pFftBuffer[it2+1] = pFMagSpectrum[i] * sin(pFPhaseSpectrum[i]);
     }    
 
     /* compute IFFT */
     sms_ifft(sizeFft, pFftBuffer); 
 
-    /* assume the output array has been taken care off */
-    /* \todo is a seperate pFftBuffer necessary here?
-       it seems like multiplying the window into the waveform
-       would be fine, without pFftBuffer */
+    /* assume that the output array does not need to be cleared */
     for(i = 0; i < sizeWave; i++)
         pFWaveform[i] +=  (pFftBuffer[i] * pFWindow[i] * .5);
 
diff --git a/sms/synthesis.c b/sms/synthesis.c
index 9523ee6..bc739d0 100644
--- a/sms/synthesis.c
+++ b/sms/synthesis.c
@@ -120,14 +120,12 @@ static int StocSynthApprox(SMS_Data *pSmsData, SMS_SynthParams *pSynthParams)
     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)
+    if(*(pSmsData->pFStocGain) <= 0)
         return 0;
 
-    sizeSpec1Used = sizeSpec1 * pSynthParams->iSamplingRate /
-                    pSynthParams->iOriginalSRate;
+    sizeSpec1Used = sizeSpec1 * pSynthParams->iSamplingRate / pSynthParams->iOriginalSRate;
 
     /* sizeSpec1Used cannot be more than what is available  \todo check by graph */
     if(sizeSpec1Used  > sizeSpec1) sizeSpec1Used = sizeSpec1;
@@ -150,74 +148,79 @@ static int StocSynthApprox(SMS_Data *pSmsData, SMS_SynthParams *pSynthParams)
  *
  * \param residualParams Parameters and memory for residual synthesis
  */
-void sms_approxResidual(SMS_ResidualParams *residualParams)
+void sms_approxResidual(int sizeResidual, sfloat* residual,
+                        int sizeApprox, sfloat* approx, 
+                        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 < residualParams->sizeStocMagSpectrum; i++)
-        fMag += (residualParams->stocMagSpectrum[i] * pAnalParams->stocMagSpectrum[i]);
 
-    /* if no gain or no coefficients return  */
-    sfloat stocGain = fMag / residualParams->sizeStocMagSpectrum;
-    if(stocGain <= 0)
-       return;
+    /* shift buffers */
+    memcpy(residualParams->residual,
+           residualParams->residual + residualParams->hopSize,
+           sizeof(sfloat) * residualParams->hopSize);
+    memcpy(residualParams->residual + residualParams->hopSize, residual,
+           sizeof(sfloat) * residualParams->hopSize);
 
-    int i, sizeSpec1Used;
-    int sizeSpec1 = residualParams->nCoeffs;
-    /*int sizeSpec2 = pSynthParams->sizeHop;*/
-    int sizeSpec2 = residualParams->residualSize;
-    int sizeFft = sizeSpec2 << 1; /* 50% overlap, so sizeFft is 2x sizeHop */
+    memcpy(residualParams->approx,
+           residualParams->approx + residualParams->hopSize,
+           sizeof(sfloat) * residualParams->hopSize);
+    memset(residualParams->approx + residualParams->hopSize, 0, 
+           sizeof(sfloat) * residualParams->hopSize);
 
-    /*sizeSpec1Used = sizeSpec1 * pSynthParams->iSamplingRate / pSynthParams->iOriginalSRate;*/
+    sms_spectrumMag(residualParams->residualSize,
+                    residualParams->residual,
+                    residualParams->fftWindow,
+                    residualParams->sizeStocMagSpectrum,
+                    residualParams->stocMagSpectrum,
+                    residualParams->fftBuffer);
+
+    if(residualParams->sizeStocMagSpectrum != residualParams->nCoeffs)
+    {
+        sms_spectralApprox(residualParams->stocMagSpectrum,
+                           residualParams->sizeStocMagSpectrum,
+                           residualParams->sizeStocMagSpectrum,
+                           residualParams->stocCoeffs,
+                           residualParams->nCoeffs,
+                           residualParams->nCoeffs,
+                           residualParams->approxEnvelope);
 
-    /*[> sizeSpec1Used cannot be more than what is available  \todo check by graph <]*/
-    /*if(sizeSpec1Used  > sizeSpec1) sizeSpec1Used = sizeSpec1;*/
+        sms_spectralApprox(residualParams->stocCoeffs,
+                           residualParams->nCoeffs,
+                           residualParams->nCoeffs,
+                           residualParams->stocMagSpectrum,
+                           residualParams->sizeStocMagSpectrum,
+                           residualParams->sizeStocMagSpectrum,
+                           residualParams->approxEnvelope);
+    }
 
-    /*sms_spectralApprox(pSmsData->pFStocCoeff, sizeSpec1, sizeSpec1Used,*/
-    /*                   pSynthParams->pMagBuff, sizeSpec2, sizeSpec1Used,*/
-    /*                   pSynthParams->approxEnvelope);*/
+    /* generate random phases */
+    for(i = 0; i < residualParams->sizeStocMagSpectrum; i++)
+        residualParams->stocPhaseSpectrum[i] =  TWO_PI * sms_random();
 
-    /*[> generate random phases <]*/
-    /*for(i = 0; i < sizeSpec2; i++)*/
-    /*    pSynthParams->pPhaseBuff[i] =  TWO_PI * sms_random();*/
+    /* IFFT with 50% overlap */
+    sms_invQuickSpectrumW(residualParams->stocMagSpectrum,
+                          residualParams->stocPhaseSpectrum,
+                          residualParams->sizeStocMagSpectrum*2,
+                          residualParams->approx,
+                          residualParams->residualSize,
+                          residualParams->ifftWindow,
+                          residualParams->fftBuffer);
 
-    /*sms_invQuickSpectrumW(pSynthParams->pMagBuff, pSynthParams->pPhaseBuff,*/
-    /*                      sizeFft, pSynthParams->pSynthBuff, sizeFft,*/
-    /*                      pSynthParams->pFStocWindow, pSynthParams->pSpectra);*/
+    /* output */
+    for(i = 0; i < sizeApprox; i++)
+        approx[i] = residualParams->approx[i] * residualParams->windowScale;
 }
 
 /*! \brief  synthesizes one frame of SMS data
  *
- * \param pSmsData      input SMS data
- * \param pFSynthesis      output sound buffer  
- * \param pSynthParams   synthesis parameters
+ * \param pSmsData     input SMS data
+ * \param pFSynthesis  output sound buffer  
+ * \param pSynthParams synthesis parameters
  */
 void sms_synthesize(SMS_Data *pSmsData, sfloat *pFSynthesis,  SMS_SynthParams *pSynthParams)
 {
-    int i, k;
+    int i;
     int sizeHop = pSynthParams->sizeHop;
-    int sizeFft = sizeHop << 1;
 
     memcpy(pSynthParams->pSynthBuff, (sfloat *)(pSynthParams->pSynthBuff+sizeHop), 
            sizeof(sfloat) * sizeHop);