diff options
Diffstat (limited to 'sms/peakDetection.c')
| -rw-r--r-- | sms/peakDetection.c | 211 |
1 files changed, 211 insertions, 0 deletions
diff --git a/sms/peakDetection.c b/sms/peakDetection.c new file mode 100644 index 0000000..ad94332 --- /dev/null +++ b/sms/peakDetection.c @@ -0,0 +1,211 @@ +/* + * Copyright (c) 2008 MUSIC TECHNOLOGY GROUP (MTG) + * UNIVERSITAT POMPEU FABRA + * + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + * + */ +/*! \file peakDetection.c + * \brief peak detection algorithm and functions + */ + +#include "sms.h" + +/*! \brief function used for the parabolic interpolation of the spectral peaks + * + * it performs the interpolation in a log scale and + * stores the location in pFDiffFromMax and + * + * \param fMaxVal value of max bin + * \param fLeftBinVal value for left bin + * \param fRightBinVal value for right bin + * \param pFDiffFromMax location of the tip as the difference from the top bin + * \return the peak height + */ +static sfloat PeakInterpolation (sfloat fMaxVal, sfloat fLeftBinVal, + sfloat fRightBinVal, sfloat *pFDiffFromMax) +{ + /* get the location of the tip of the parabola */ + *pFDiffFromMax = (.5 * (fLeftBinVal - fRightBinVal) / + (fLeftBinVal - (2*fMaxVal) + fRightBinVal)); + /* return the value at the tip */ + return(fMaxVal - (.25 * (fLeftBinVal - fRightBinVal) * + *pFDiffFromMax)); +} + +/*! \brief detect the next local maximum in the spectrum + * + * stores the value in pFMaxVal + * + * \todo export this to sms.h and wrap in pysms + * + * \param pFMagSpectrum magnitude spectrum + * \param iHighBinBound highest bin to search + * \param pICurrentLoc current bin location + * \param pFMaxVal value of the maximum found + * \param fMinPeakMag minimum magnitude to accept a peak + * \return the bin location of the maximum + */ +static int FindNextMax ( sfloat *pFMagSpectrum, int iHighBinBound, + int *pICurrentLoc, sfloat *pFMaxVal, sfloat fMinPeakMag) +{ + int iCurrentBin = *pICurrentLoc; + sfloat fPrevVal = pFMagSpectrum[iCurrentBin - 1]; + sfloat fCurrentVal = pFMagSpectrum[iCurrentBin]; + sfloat fNextVal = (iCurrentBin >= iHighBinBound) + ? 0 : pFMagSpectrum[iCurrentBin + 1]; + + /* try to find a local maximum */ + while (iCurrentBin <= iHighBinBound) + { + if (fCurrentVal > fMinPeakMag && + fCurrentVal >= fPrevVal && + fCurrentVal >= fNextVal) + break; + iCurrentBin++; + fPrevVal = fCurrentVal; + fCurrentVal = fNextVal; + fNextVal = pFMagSpectrum[1+iCurrentBin]; + } + /* save the current location, value of maximum and return */ + /* location of max */ + *pICurrentLoc = iCurrentBin + 1; + *pFMaxVal = fCurrentVal; + return(iCurrentBin); +} + +/*! \brief function to detect the next spectral peak + * + * \param pFMagSpectrum magnitude spectrum + * \param iHighestBin highest bin to search + * \param pICurrentLoc current bin location + * \param pFPeakMag magnitude value of peak + * \param pFPeakLoc location of peak + * \param fMinPeakMag minimum magnitude to accept a peak + * \return 1 if found, 0 if not + */ +static int FindNextPeak (sfloat *pFMagSpectrum, int iHighestBin, + int *pICurrentLoc, sfloat *pFPeakMag, sfloat *pFPeakLoc, + sfloat fMinPeakMag) +{ + int iPeakBin = 0; /* location of the local peak */ + sfloat fPeakMag = 0; /* value of local peak */ + + /* keep trying to find a good peak while inside the freq range */ + while ((iPeakBin = FindNextMax(pFMagSpectrum, iHighestBin, + pICurrentLoc, &fPeakMag, fMinPeakMag)) + <= iHighestBin) + { + /* get the neighboring samples */ + sfloat fDiffFromMax = 0; + sfloat fLeftBinVal = pFMagSpectrum[iPeakBin - 1]; + sfloat fRightBinVal = pFMagSpectrum[iPeakBin + 1]; + if (fLeftBinVal <= 0 || fRightBinVal <= 0) //ahah! there you are! + continue; + /* interpolate the spectral samples to obtain + a more accurate magnitude and freq */ + *pFPeakMag = PeakInterpolation (fPeakMag, fLeftBinVal, + fRightBinVal, &fDiffFromMax); + *pFPeakLoc = iPeakBin + fDiffFromMax; + return (1); + } + /* if it does not find a peak return 0 */ + return (0); +} + +/*! \brief get the corresponding phase value for a given peak + * + * performs linear interpolation for a more accurate phase + + * \param pPhaseSpectrum phase spectrum + * \param fPeakLoc location of peak + * \return the phase value + */ +static sfloat GetPhaseVal (sfloat *pPhaseSpectrum, sfloat fPeakLoc) +{ + int bin = (int) fPeakLoc; + sfloat fFraction = fPeakLoc - bin, + fLeftPha = pPhaseSpectrum[bin], + fRightPha = pPhaseSpectrum[bin+1]; + + /* check for phase wrapping */ + if (fLeftPha - fRightPha > 1.5 * PI) + fRightPha += TWO_PI; + else if (fRightPha - fLeftPha > 1.5 * PI) + fLeftPha += TWO_PI; + + /* return interpolated phase */ + return (fLeftPha + fFraction * (fRightPha - fLeftPha)); +} + +/*! \brief find the prominent spectral peaks + * + * uses a dB spectrum + * + * \param sizeSpec size of magnitude spectrum + * \param pMag pointer to power spectrum + * \param pPhase pointer to phase spectrum + * \param pSpectralPeaks pointer to array of peaks + * \param pPeakParams peak detection parameters + * \return the number of peaks found + */ +int sms_detectPeaks (int sizeSpec, sfloat *pMag, sfloat *pPhase, + SMS_Peak *pSpectralPeaks, SMS_PeakParams *pPeakParams) +{ + static int iFirstBin, iHighestBin, sizeFft; + static sfloat fInvSizeFft; + static int sizeSpecStatic = 0; + + /* allocate memory if sizeSpec has changed or this is the first run */ + if(sizeSpecStatic != sizeSpec) + { + sizeSpecStatic = sizeSpec; + //printf("sizeSpecStatic: %d \n", sizeSpecStatic); + sizeFft = sizeSpec << 1; + fInvSizeFft = 1.0 / sizeFft; + /* make sure to start on the 2nd bin so interpolation is possible QUESTION: why not allow a peak in bin 1 or 0? */ + /* rte: changed the first argument of MAX from 2 to 1 */ + iFirstBin = MAX (1, sizeFft * pPeakParams->fLowestFreq / + pPeakParams->iSamplingRate); + iHighestBin = MIN (sizeSpec-1, + sizeFft * pPeakParams->fHighestFreq / + pPeakParams->iSamplingRate); + } + + /* clear peak structure */ + memset (pSpectralPeaks, 0, pPeakParams->iMaxPeaks * sizeof(SMS_Peak)); + + /* set starting search values */ + int iCurrentLoc = iFirstBin; + int iPeak = 0; /* index for spectral search */ + sfloat fPeakMag = 0.0; /* magnitude of peak */ + sfloat fPeakLoc = 0.0; /* location of peak */ + + /* find peaks */ + while ((iPeak < pPeakParams->iMaxPeaks) && + (FindNextPeak(pMag, iHighestBin, + &iCurrentLoc, &fPeakMag, &fPeakLoc, pPeakParams->fMinPeakMag) == 1)) + { + /* store peak values */ + pSpectralPeaks[iPeak].fFreq = pPeakParams->iSamplingRate * fPeakLoc * fInvSizeFft; + pSpectralPeaks[iPeak].fMag = fPeakMag; + pSpectralPeaks[iPeak].fPhase = GetPhaseVal(pPhase, fPeakLoc); + iPeak++; + } + + /* return the number of peaks found */ + return (iPeak); +} |