diff options
Diffstat (limited to 'sms/spectrum.c')
| -rw-r--r-- | sms/spectrum.c | 403 |
1 files changed, 149 insertions, 254 deletions
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); */ -/* } */ |