Change to using distutils.

Currently only builds the simplsndobj module

1 files changed, 202 insertions, 0 deletions

diff --git a/src/sms/peakDetection.c b/src/sms/peakDetection.c
new file mode 100644
index 0000000..4015c27
--- /dev/null
+++ b/src/sms/peakDetection.c
@@ -0,0 +1,202 @@
+/* 
+ * 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 pAnalParams peak detection parameters
+ * \return the number of peaks found
+ */
+int sms_detectPeaks(int sizeSpec, sfloat *pMag, sfloat *pPhase,
+                    SMS_Peak *pSpectralPeaks, SMS_AnalParams *pAnalParams)
+{
+    int sizeFft = sizeSpec << 1;
+    sfloat fInvSizeFft = 1.0 / sizeFft;
+    int iFirstBin = MAX(1, sizeFft * pAnalParams->fLowestFreq / pAnalParams->iSamplingRate);
+    int iHighestBin = MIN(sizeSpec-1, 
+                          sizeFft * pAnalParams->fHighestFreq / pAnalParams->iSamplingRate);
+    int iPeak = 0;
+
+    /* clear peak structure */
+    for(iPeak = 0; iPeak < pAnalParams->maxPeaks; iPeak++)
+    {
+        pSpectralPeaks[iPeak].fFreq = 0.0;
+        pSpectralPeaks[iPeak].fMag = 0.0;
+        pSpectralPeaks[iPeak].fPhase = 0.0;
+    }
+
+    /* set starting search values */
+    int iCurrentLoc = iFirstBin;
+    sfloat fPeakMag = 0.0;  /* magnitude of peak */
+    sfloat fPeakLoc = 0.0;  /* location of peak */
+
+    /* find peaks */
+    iPeak = 0;
+    while((iPeak < pAnalParams->maxPeaks) &&
+          (FindNextPeak(pMag, iHighestBin, &iCurrentLoc, &fPeakMag, 
+                        &fPeakLoc, pAnalParams->fMinPeakMag) == 1))
+    {
+        /* store peak values */
+        pSpectralPeaks[iPeak].fFreq = pAnalParams->iSamplingRate * fPeakLoc * fInvSizeFft;
+        pSpectralPeaks[iPeak].fMag = fPeakMag;
+        pSpectralPeaks[iPeak].fPhase = GetPhaseVal(pPhase, fPeakLoc);
+        iPeak++;
+    }
+
+    /* return the number of peaks found */
+    return iPeak;
+}