From 5e25edb1b78f28bee09bd7513a80fb500bbd59c4 Mon Sep 17 00:00:00 2001
From: John Glover <glover.john@gmail.com>
Date: Wed, 24 Nov 2010 23:26:43 +0000
Subject: Updated libsms to the latest version (1.15), includes some memory
 management improvements/simplifications but the basic API is the same. Also
 started updating the unit tests, which will now use the nose framework

---
 sms/spectrum.c | 403 +++++++++++++++++++++------------------------------------
 1 file changed, 149 insertions(+), 254 deletions(-)

(limited to 'sms/spectrum.c')

diff --git a/sms/spectrum.c b/sms/spectrum.c
index c7d815a..356811f 100644
--- a/sms/spectrum.c
+++ b/sms/spectrum.c
@@ -25,58 +25,39 @@
 
 /*! \brief compute a complex spectrum from a waveform 
  *              
- * \param sizeWindow	           size of analysis window
- * \param pWaveform	           pointer to input waveform 
- * \param pWindow	                   pointer to input window
- * \param sizeMag	                   size of output magnitude and phase spectrums
+ * \param sizeWindow               size of analysis window
+ * \param pWaveform            pointer to input waveform 
+ * \param pWindow                      pointer to input window
+ * \param sizeMag                      size of output magnitude and phase spectrums
  * \param pMag                          pointer to output magnitude spectrum 
  * \param pPhase                       pointer to output phase spectrum 
  * \return sizeFft, -1 on error \todo remove this return
  */
-int sms_spectrum (int sizeWindow, sfloat *pWaveform, sfloat *pWindow, int sizeMag, 
-                  sfloat *pMag, sfloat *pPhase)
+int sms_spectrum(int sizeWindow, sfloat *pWaveform, sfloat *pWindow, int sizeMag, 
+                 sfloat *pMag, sfloat *pPhase, sfloat *pFftBuffer)
 {
-	int sizeFft = sizeMag << 1;
-        int i, it2;
-        int err = 0;
-        sfloat fReal, fImag;
-  
-        static sfloat *pFftBuffer;
-        static int sizeFftArray = 0;
-        /* if new size fft is larger than old, allocate more memory */
-        if(sizeFftArray < sizeFft)
-        {
-                if(sizeFftArray != 0) free(pFftBuffer);
-                sizeFftArray = sms_power2(sizeFft);
-                if(sizeFftArray != sizeFft)
-                {
-                        sms_error("bad fft size, incremented to power of 2");
-                        err = -1;
-                }
-                if ((pFftBuffer = (sfloat *) malloc(sizeFft * sizeof(sfloat))) == NULL)
-                {
-                        sms_error("could not allocate memory for fft array");
-                        return(-1);
-                }
-        }
-        memset(pFftBuffer, 0, sizeFft * sizeof(sfloat));
-
-	/* apply window to waveform and center window around 0 (zero-phase windowing)*/
-        sms_windowCentered(sizeWindow, pWaveform, pWindow, sizeFft, pFftBuffer);
-
-        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); /* \todo why is fImag negated? */
-	}
-        
-	return (sizeFft);
+    int i, it2;
+    int err = 0;
+    sfloat fReal, fImag;
+    
+    int sizeFft = sizeMag << 1;
+    memset(pFftBuffer, 0, sizeFft * sizeof(sfloat));
+
+    /* apply window to waveform and center window around 0 (zero-phase windowing)*/
+    sms_windowCentered(sizeWindow, pWaveform, pWindow, sizeFft, pFftBuffer);
+    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); /* \todo why is fImag negated? */
+    }
+
+    return sizeFft;
 }
 
 /*! \brief compute the spectrum Magnitude of a waveform
@@ -85,63 +66,42 @@ int sms_spectrum (int sizeWindow, sfloat *pWaveform, sfloat *pWindow, int sizeMa
  * performs a zero-padded FFT (if sizeMag*2 > sizeWindow).
  * The spectra is then converted magnitude (RMS).
  *              
- * \param sizeWindow	       size of analysis window / input wavefrom
- * \param pWaveform	       pointer to input waveform
- * \param pWindow	       pointer to analysis window 
- * \param sizeMag		       size of output magnitude spectrum 
+ * \param sizeWindow           size of analysis window / input wavefrom
+ * \param pWaveform        pointer to input waveform
+ * \param pWindow          pointer to analysis window 
+ * \param sizeMag              size of output magnitude spectrum 
  * \param pMag     pointer to output magnitude spectrum 
  * \return 0 on success, -1 on error
  */
-int sms_spectrumMag (int sizeWindow, sfloat *pWaveform, sfloat *pWindow,  
-                     int sizeMag, sfloat *pMag)
+int sms_spectrumMag(int sizeWindow, sfloat *pWaveform, sfloat *pWindow,  
+                    int sizeMag, sfloat *pMag, sfloat *pFftBuffer)
 {
-        int i,it2;
-        int sizeFft = sizeMag << 1;
-        int err = 0;
-        sfloat fReal, fImag;
-
-        static sfloat *pFftBuffer;
-        static int sizeFftArray = 0;
-
-        if(sizeFftArray != sizeFft)
-        {
-                if(sizeFftArray != 0) free(pFftBuffer);
-                sizeFftArray = sms_power2(sizeFft);
-                if(sizeFftArray != sizeFft)
-                {
-                        sms_error("bad fft size, incremented to power of 2");
-                        err = -1;
-                }
-                if ((pFftBuffer = (sfloat *) malloc(sizeFft * sizeof(sfloat))) == NULL)
-                {
-                        sms_error("could not allocate memory for fft array");
-                        return(-1);
-                }
-        }
-	/* apply window to waveform, zero the rest of the array */
-        //memset(pFftBuffer, 0, sizeFft * sizeof(sfloat));
-	for (i = 0; i < sizeWindow; i++)
-                pFftBuffer[i] =  pWindow[i] * pWaveform[i];
-        for(i = sizeWindow; i < sizeFft; i++)
-                pFftBuffer[i]  = 0.;
-  
-	/* compute real FFT */
-        sms_fft(sizeFft, pFftBuffer); 
-  
-	/* convert from rectangular to polar coordinates */
-	for (i=0; i<sizeMag; i++)
-	{
-		it2 = i << 1;
-		fReal = pFftBuffer[it2];
-		fImag = pFftBuffer[it2+1];
-                pMag[i] = sqrtf(fReal * fReal + fImag * fImag);
-	}
-
-	//return (err);
-	return (sizeFft);
+    int i,it2;
+    int sizeFft = sizeMag << 1;
+    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 = sizeWindow; i < sizeFft; i++)
+        pFftBuffer[i]  = 0.;
+
+    /* compute real FFT */
+    sms_fft(sizeFft, pFftBuffer); 
+
+    /* convert from rectangular to polar coordinates */
+    for (i=0; i<sizeMag; i++)
+    {
+        it2 = i << 1;
+        fReal = pFftBuffer[it2];
+        fImag = pFftBuffer[it2+1];
+        pMag[i] = sqrtf(fReal * fReal + fImag * fImag);
+    }
+
+    return sizeFft;
 }
 
- 
+
 /*! \brief function for a quick inverse spectrum, windowed
  *
  *  Not done yet, but this will be a function that is the inverse of
@@ -150,194 +110,129 @@ int sms_spectrumMag (int sizeWindow, sfloat *pWaveform, sfloat *pWindow,
  * function to perform the inverse FFT, windowing the output
  * sfloat *pFMagSpectrum        input magnitude spectrum
  * sfloat *pFPhaseSpectrum      input phase spectrum
- * int sizeFft		       size of FFT
- * sfloat *pFWaveform	       output waveform
+ * int sizeFft             size of FFT
+ * sfloat *pFWaveform          output waveform
  * int sizeWave                size of output waveform
- * sfloat *pFWindow	       synthesis window
+ * sfloat *pFWindow        synthesis window
  */
-int sms_invSpectrum (int sizeWaveform, sfloat *pWaveform, sfloat *pWindow ,
-                     int sizeMag, sfloat *pMag, sfloat *pPhase)
+int sms_invSpectrum(int sizeWaveform, sfloat *pWaveform, sfloat *pWindow,
+                    int sizeMag, sfloat *pMag, sfloat *pPhase, sfloat *pFftBuffer)
 {
-	int i;
-        int sizeFft = sizeMag << 1;
-        int err = 0;
-        static sfloat *pFftBuffer;
-        static int sizeFftArray = 0;
+    int i;
+    int sizeFft = sizeMag << 1;
 
-        if(sizeFftArray != sizeFft)
-        {
-                if(sizeFftArray != 0) free(pFftBuffer);
-                sizeFftArray = sms_power2(sizeFft);
-                if(sizeFftArray != sizeFft)
-                {
-                        sms_error("bad fft size, incremented to power of 2");
-                        err = -1;
-                }
-                if ((pFftBuffer = (sfloat *) malloc(sizeFft * sizeof(sfloat))) == NULL)
-                {
-                        sms_error("could not allocate memory for fft array");
-                        return(-1);
-                }
-        }
+    sms_PolarToRect(sizeMag, pFftBuffer, pMag, pPhase);
+    sms_ifft(sizeFft, pFftBuffer); 
 
-        sms_PolarToRect(sizeMag, pFftBuffer, pMag, pPhase);
-	/* compute IFFT */
-        sms_ifft(sizeFft, pFftBuffer); 
+    /* assume the output array has been taken care off */
+    /* before, this was multiplied by .5, why? */
+    for (i = 0; i < sizeWaveform; i++)
+        //pWaveform[i] +=  pFftBuffer[i] * pWindow[i];
+        pWaveform[i] = pFftBuffer[i];
 
- 	/* assume the output array has been taken care off */
-        /* before, this was multiplied by .5, why? */
-	for (i = 0; i < sizeWaveform; i++)
-		//pWaveform[i] +=  pFftBuffer[i] * pWindow[i];
-                pWaveform[i] =  pFftBuffer[i];
-  
-	return (sizeFft);
-//	return (err);
+    return sizeFft;
 }
+
 /*! \brief function for a quick inverse spectrum, windowed
  * function to perform the inverse FFT, windowing the output
  * sfloat *pFMagSpectrum        input magnitude spectrum
  * sfloat *pFPhaseSpectrum      input phase spectrum
- * int sizeFft		       size of FFT
- * sfloat *pFWaveform	       output waveform
+ * int sizeFft             size of FFT
+ * sfloat *pFWaveform          output waveform
  * int sizeWave                size of output waveform
- * sfloat *pFWindow	       synthesis window
+ * sfloat *pFWindow        synthesis window
  */
-int sms_invQuickSpectrumW (sfloat *pFMagSpectrum, sfloat *pFPhaseSpectrum, 
-                           int sizeFft, sfloat *pFWaveform, int sizeWave,
-                           sfloat *pFWindow)
+int sms_invQuickSpectrumW(sfloat *pFMagSpectrum, sfloat *pFPhaseSpectrum, 
+                          int sizeFft, sfloat *pFWaveform, int sizeWave,
+                          sfloat *pFWindow, sfloat* pFftBuffer)
 {
-	int sizeMag = sizeFft >> 1, i, it2;
-	sfloat *pFBuffer, fPower;
-  
-	/* allocate buffer */    
-	if ((pFBuffer = (sfloat *) calloc(sizeFft, sizeof(sfloat))) == NULL)
-		return -1;
-
-	/* convert from polar coordinates to rectangular  */
-	for (i = 0; i<sizeMag; i++)
-	{
-		it2 = i << 1;
-		fPower = pFMagSpectrum[i];
-		pFBuffer[it2] =  fPower * cos (pFPhaseSpectrum[i]);
-		pFBuffer[it2+1] = fPower * sin (pFPhaseSpectrum[i]);
-	}    
-	/* compute IFFT */
-        sms_ifft(sizeFft, pFBuffer); 
-
- 	/* assume the output array has been taken care off */
-        /* \todo is a seperate pFBuffer necessary here?
-           it seems like multiplying the window into the waveform
-           would be fine, without pFBuffer */
-	for (i = 0; i < sizeWave; i++)
-		pFWaveform[i] +=  (pFBuffer[i] * pFWindow[i] * .5);
-
-	free (pFBuffer);
-  
-	return (sizeMag);
+    int sizeMag = sizeFft >> 1, i, it2;
+    sfloat fPower;
+
+    /* 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]);
+    }    
+
+    /* 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 */
+    for (i = 0; i < sizeWave; i++)
+        pFWaveform[i] +=  (pFftBuffer[i] * pFWindow[i] * .5);
+
+    return sizeMag;
 }
 
 /*! \brief convert spectrum from Rectangular to Polar form
  *              
- * \param sizeMag	       size of spectrum (pMag and pPhase arrays)
- * \param pRect	       pointer output spectrum in rectangular form (2x sizeSpec)
- * \param pMag	       pointer to sfloat array of magnitude spectrum
- * \param pPhase	       pointer to sfloat array of phase spectrum
+ * \param sizeMag          size of spectrum (pMag and pPhase arrays)
+ * \param pRect        pointer output spectrum in rectangular form (2x sizeSpec)
+ * \param pMag         pointer to sfloat array of magnitude spectrum
+ * \param pPhase           pointer to sfloat array of phase spectrum
  */ 
 void sms_RectToPolar( int sizeMag, sfloat *pRect, sfloat *pMag, sfloat *pPhase)
 {
-        int i, it2;
-        sfloat fReal, fImag;
-
-	for (i=0; i<sizeMag; i++)
-	{
-		it2 = i << 1;
-		fReal = pRect[it2];
-		fImag = pRect[it2+1];
-      
-                pMag[i] = sqrtf(fReal * fReal + fImag * fImag);
-                if (pPhase)
-                        pPhase[i] = atan2f(fImag, fReal);
-	}
-
-
+    int i, it2;
+    sfloat fReal, fImag;
+
+    for (i=0; i<sizeMag; i++)
+    {
+        it2 = i << 1;
+        fReal = pRect[it2];
+        fImag = pRect[it2+1];
+
+        pMag[i] = sqrtf(fReal * fReal + fImag * fImag);
+        if (pPhase)
+            pPhase[i] = atan2f(fImag, fReal);
+    }
 }
 
 /*! \brief convert spectrum from Rectangular to Polar form
  *              
- * \param sizeSpec	       size of spectrum (pMag and pPhase arrays)
- * \param pRect	       pointer output spectrum in rectangular form (2x sizeSpec)
- * \param pMag	       pointer to sfloat array of magnitude spectrum
- * \param pPhase	       pointer to sfloat array of phase spectrum
+ * \param sizeSpec         size of spectrum (pMag and pPhase arrays)
+ * \param pRect        pointer output spectrum in rectangular form (2x sizeSpec)
+ * \param pMag         pointer to sfloat array of magnitude spectrum
+ * \param pPhase           pointer to sfloat array of phase spectrum
  */ 
-void sms_PolarToRect( int sizeSpec, sfloat *pRect, sfloat *pMag, sfloat *pPhase)
+void sms_PolarToRect(int sizeSpec, sfloat *pRect, sfloat *pMag, sfloat *pPhase)
 {
-        int i, it2;
-        sfloat fMag;
-
-	for (i = 0; i<sizeSpec; i++)
-	{
-		it2 = i << 1;
-		fMag = pMag[i];
-		pRect[it2] =  fMag * cos (pPhase[i]);
-		pRect[it2+1] = fMag * sin (pPhase[i]);
-	}    
+    int i, it2;
+    sfloat fMag;
+
+    for(i = 0; i<sizeSpec; i++)
+    {
+        it2 = i << 1;
+        fMag = pMag[i];
+        pRect[it2] =  fMag * cos (pPhase[i]);
+        pRect[it2+1] = fMag * sin (pPhase[i]);
+    }    
 }
 
 /*! \brief compute magnitude spectrum of a DFT in rectangular coordinates
  *              
- * \param sizeMag	       size of output Magnitude (half of input real FFT)
- * \param pInRect	       pointer to input DFT array (real/imag sfloats)
- * \param pOutMag	       pointer to of magnitude spectrum array
+ * \param sizeMag          size of output Magnitude (half of input real FFT)
+ * \param pInRect          pointer to input DFT array (real/imag sfloats)
+ * \param pOutMag          pointer to of magnitude spectrum array
  */
-void sms_spectrumRMS( int sizeMag, sfloat *pInRect, sfloat *pOutMag)
+void sms_spectrumRMS(int sizeMag, sfloat *pInRect, sfloat *pOutMag)
 {
-        int i, it2;
-        sfloat fReal, fImag;
-
-	for (i=0; i<sizeMag; i++)
-	{
-		it2 = i << 1;
-		fReal = pInRect[it2];
-		fImag = pInRect[it2+1];
-                pOutMag[i] = sqrtf(fReal * fReal + fImag * fImag);
-	}
+    int i, it2;
+    sfloat fReal, fImag;
+
+    for(i=0; i<sizeMag; i++)
+    {
+        it2 = i << 1;
+        fReal = pInRect[it2];
+        fImag = pInRect[it2+1];
+        pOutMag[i] = sqrtf(fReal * fReal + fImag * fImag);
+    }
 }
 
-
-/*! \brief convert from Polar spectrum to waveform
- * function to perform the inverse FFT
- * sfloat *pFMagSpectrum        input magnitude spectrum
- * sfloat *pFPhaseSpectrum      input phase spectrum
- * int sizeFft                 size of FFT
- * sfloat *pFWaveform           output waveform
- * int sizeWave                size of output waveform
- */
-/* int sms_invSpectrum (sfloat *pFMagSpectrum, sfloat *pFPhaseSpectrum, */
-/*                           int sizeFft, sfloat *pFWaveform, int sizeWave) */
-/* { */
-/* 	int sizeMag = sizeFft >> 1, i, it2; */
-/* 	sfloat *pFBuffer, fPower; */
-  
-/* 	/\* allocate buffer *\/ */
-/* 	if ((pFBuffer = (sfloat *) calloc(sizeFft+1, sizeof(sfloat))) == NULL) */
-/* 		return -1; */
-   
-/* 	/\* convert from polar coordinates to rectangular  *\/ */
-/* 	for (i = 0; i < sizeMag; i++) */
-/* 	{ */
-/* 		it2 = i << 1; */
-/* 		fPower = pFMagSpectrum[i]; */
-/* 		pFBuffer[it2] =  fPower * cos (pFPhaseSpectrum[i]); */
-/* 		pFBuffer[it2+1] = fPower * sin (pFPhaseSpectrum[i]); */
-/* 	} */
-/* 	/\* compute IFFT *\/ */
-/*         sms_ifft(sizeFft, pFBuffer);  */
- 
-/* 	/\* assume the output array has been taken care off *\/ */
-/* 	for (i = 0; i < sizeWave; i++) */
-/* 		pFWaveform[i] +=  pFBuffer[i]; */
- 
-/* 	free(pFBuffer); */
-  
-/* 	return (sizeMag); */
-/* } */
-- 
cgit v1.2.3