/external/pysoundtouch14/libsoundtouch/PeakFinder.cpp

  1////////////////////////////////////////////////////////////////////////////////
  2///
  3/// Peak detection routine. 
  4///
  5/// The routine detects highest value on an array of values and calculates the 
  6/// precise peak location as a mass-center of the 'hump' around the peak value.
  7///
  8/// Author        : Copyright (c) Olli Parviainen
  9/// Author e-mail : oparviai 'at' iki.fi
 10/// SoundTouch WWW: http://www.surina.net/soundtouch
 11///
 12////////////////////////////////////////////////////////////////////////////////
 13//
 14// Last changed  : $Date: 2008-12-25 19:54:41 +0200 (Thu, 25 Dec 2008) $
 15// File revision : $Revision: 4 $
 16//
 17// $Id: PeakFinder.cpp 43 2008-12-25 17:54:41Z oparviai $
 18//
 19////////////////////////////////////////////////////////////////////////////////
 20//
 21// License :
 22//
 23//  SoundTouch audio processing library
 24//  Copyright (c) Olli Parviainen
 25//
 26//  This library is free software; you can redistribute it and/or
 27//  modify it under the terms of the GNU Lesser General Public
 28//  License as published by the Free Software Foundation; either
 29//  version 2.1 of the License, or (at your option) any later version.
 30//
 31//  This library is distributed in the hope that it will be useful,
 32//  but WITHOUT ANY WARRANTY; without even the implied warranty of
 33//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 34//  Lesser General Public License for more details.
 35//
 36//  You should have received a copy of the GNU Lesser General Public
 37//  License along with this library; if not, write to the Free Software
 38//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 39//
 40////////////////////////////////////////////////////////////////////////////////
 41
 42#include <math.h>

 43#include <assert.h>

 44
 45#include "PeakFinder.h"

 46
 47using namespace soundtouch;
 48
 49#define max(x, y) (((x) > (y)) ? (x) : (y))

 50
 51
 52PeakFinder::PeakFinder()
 53{
 54}
 55
 56
 57// Finds 'ground level' of a peak hump by starting from 'peakpos' and proceeding
 58// to direction defined by 'direction' until next 'hump' after minimum value will 
 59// begin
 60int PeakFinder::findGround(const float *data, int peakpos, int direction) const
 61{
 62    float refvalue;
 63    int lowpos;
 64    int pos;
 65    int climb_count;
 66    float delta;
 67
 68    climb_count = 0;
 69    refvalue = data[peakpos];
 70    lowpos = peakpos;
 71
 72    pos = peakpos;
 73
 74    while ((pos > minPos) && (pos < maxPos))
 75    {
 76        int prevpos;
 77
 78        prevpos = pos;
 79        pos += direction;
 80
 81        // calculate derivate
 82        delta = data[pos] - data[prevpos];
 83        if (delta <= 0)
 84        {
 85            // going downhill, ok
 86            if (climb_count)
 87            {
 88                climb_count --;  // decrease climb count
 89            }
 90
 91            // check if new minimum found
 92            if (data[pos] < refvalue)
 93            {
 94                // new minimum found
 95                lowpos = pos;
 96                refvalue = data[pos];
 97            }
 98        }
 99        else
100        {
101            // going uphill, increase climbing counter
102            climb_count ++;
103            if (climb_count > 5) break;    // we've been climbing too long => it's next uphill => quit
104        }
105    }
106    return lowpos;
107}
108
109
110// Find offset where the value crosses the given level, when starting from 'peakpos' and
111// proceeds to direction defined in 'direction'
112int PeakFinder::findCrossingLevel(const float *data, float level, int peakpos, int direction) const
113{
114    float peaklevel;
115    int pos;
116
117    peaklevel = data[peakpos];
118    assert(peaklevel >= level);
119    pos = peakpos;
120    while ((pos >= minPos) && (pos < maxPos))
121    {
122        if (data[pos + direction] < level) return pos;   // crossing found
123        pos += direction;
124    }
125    return -1;  // not found
126}
127
128
129// Calculates the center of mass location of 'data' array items between 'firstPos' and 'lastPos'
130double PeakFinder::calcMassCenter(const float *data, int firstPos, int lastPos) const
131{
132    int i;
133    float sum;
134    float wsum;
135
136    sum = 0;
137    wsum = 0;
138    for (i = firstPos; i <= lastPos; i ++)
139    {
140        sum += (float)i * data[i];
141        wsum += data[i];
142    }
143    return sum / wsum;
144}
145
146
147
148/// get exact center of peak near given position by calculating local mass of center
149double PeakFinder::getPeakCenter(const float *data, int peakpos)
150{
151    float peakLevel;            // peak level
152    int crosspos1, crosspos2;   // position where the peak 'hump' crosses cutting level
153    float cutLevel;             // cutting value
154    float groundLevel;          // ground level of the peak
155    int gp1, gp2;               // bottom positions of the peak 'hump'
156
157    // find ground positions.
158    gp1 = findGround(data, peakpos, -1);
159    gp2 = findGround(data, peakpos, 1);
160
161    groundLevel = max(data[gp1], data[gp2]);
162    peakLevel = data[peakpos];
163
164    if (groundLevel < 1e-6) return 0;                // ground level too small => detection failed
165    if ((peakLevel / groundLevel) < 1.3) return 0;   // peak less than 30% of the ground level => no good peak detected
166
167    // calculate 70%-level of the peak
168    cutLevel = 0.70f * peakLevel + 0.30f * groundLevel;
169    // find mid-level crossings
170    crosspos1 = findCrossingLevel(data, cutLevel, peakpos, -1);
171    crosspos2 = findCrossingLevel(data, cutLevel, peakpos, 1);
172
173    if ((crosspos1 < 0) || (crosspos2 < 0)) return 0;   // no crossing, no peak..
174
175    // calculate mass center of the peak surroundings
176    return calcMassCenter(data, crosspos1, crosspos2);
177}
178
179
180
181double PeakFinder::detectPeak(const float *data, int minPos, int maxPos) 
182{
183
184    int i;
185    int peakpos;                // position of peak level
186    double highPeak, peak;
187
188    this->minPos = minPos;
189    this->maxPos = maxPos;
190
191    // find absolute peak
192    peakpos = minPos;
193    peak = data[minPos];
194    for (i = minPos + 1; i < maxPos; i ++)
195    {
196        if (data[i] > peak) 
197        {
198            peak = data[i];
199            peakpos = i;
200        }
201    }
202    
203    // Calculate exact location of the highest peak mass center
204    highPeak = getPeakCenter(data, peakpos);
205    peak = highPeak;
206
207    // Now check if the highest peak were in fact harmonic of the true base beat peak 
208    // - sometimes the highest peak can be Nth harmonic of the true base peak yet 
209    // just a slightly higher than the true base
210    for (i = 2; i < 10; i ++)
211    {
212        double peaktmp, tmp;
213        int i1,i2;
214
215        peakpos = (int)(highPeak / (double)i + 0.5f);
216        if (peakpos < minPos) break;
217
218        // calculate mass-center of possible base peak
219        peaktmp = getPeakCenter(data, peakpos);
220
221        // now compare to highest detected peak
222        i1 = (int)(highPeak + 0.5);
223        i2 = (int)(peaktmp + 0.5);
224        tmp = 2 * (data[i2] - data[i1]) / (data[i2] + data[i1]);
225        if (fabs(tmp) < 0.1)
226        {
227            // The highest peak is harmonic of almost as high base peak,
228            // thus use the base peak instead
229            peak = peaktmp;
230        }
231    }
232
233    return peak;
234}
235
236