PageRenderTime 336ms CodeModel.GetById 161ms app.highlight 20ms RepoModel.GetById 151ms app.codeStats 0ms

/external/pysoundtouch14/libsoundtouch/BPMDetect.cpp

http://echo-nest-remix.googlecode.com/
C++ | 311 lines | 151 code | 60 blank | 100 comment | 16 complexity | b48719eba476d06baea1e61dbb100bb4 MD5 | raw file
  1////////////////////////////////////////////////////////////////////////////////
  2///
  3/// Beats-per-minute (BPM) detection routine.
  4///
  5/// The beat detection algorithm works as follows:
  6/// - Use function 'inputSamples' to input a chunks of samples to the class for
  7///   analysis. It's a good idea to enter a large sound file or stream in smallish
  8///   chunks of around few kilosamples in order not to extinguish too much RAM memory.
  9/// - Inputted sound data is decimated to approx 500 Hz to reduce calculation burden,
 10///   which is basically ok as low (bass) frequencies mostly determine the beat rate.
 11///   Simple averaging is used for anti-alias filtering because the resulting signal
 12///   quality isn't of that high importance.
 13/// - Decimated sound data is enveloped, i.e. the amplitude shape is detected by
 14///   taking absolute value that's smoothed by sliding average. Signal levels that
 15///   are below a couple of times the general RMS amplitude level are cut away to
 16///   leave only notable peaks there.
 17/// - Repeating sound patterns (e.g. beats) are detected by calculating short-term 
 18///   autocorrelation function of the enveloped signal.
 19/// - After whole sound data file has been analyzed as above, the bpm level is 
 20///   detected by function 'getBpm' that finds the highest peak of the autocorrelation 
 21///   function, calculates it's precise location and converts this reading to bpm's.
 22///
 23/// Author        : Copyright (c) Olli Parviainen
 24/// Author e-mail : oparviai 'at' iki.fi
 25/// SoundTouch WWW: http://www.surina.net/soundtouch
 26///
 27////////////////////////////////////////////////////////////////////////////////
 28//
 29// Last changed  : $Date: 2008-12-25 19:54:41 +0200 (Thu, 25 Dec 2008) $
 30// File revision : $Revision: 4 $
 31//
 32// $Id: BPMDetect.cpp 43 2008-12-25 17:54:41Z oparviai $
 33//
 34////////////////////////////////////////////////////////////////////////////////
 35//
 36// License :
 37//
 38//  SoundTouch audio processing library
 39//  Copyright (c) Olli Parviainen
 40//
 41//  This library is free software; you can redistribute it and/or
 42//  modify it under the terms of the GNU Lesser General Public
 43//  License as published by the Free Software Foundation; either
 44//  version 2.1 of the License, or (at your option) any later version.
 45//
 46//  This library is distributed in the hope that it will be useful,
 47//  but WITHOUT ANY WARRANTY; without even the implied warranty of
 48//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 49//  Lesser General Public License for more details.
 50//
 51//  You should have received a copy of the GNU Lesser General Public
 52//  License along with this library; if not, write to the Free Software
 53//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 54//
 55////////////////////////////////////////////////////////////////////////////////
 56
 57#include <math.h>

 58#include <assert.h>

 59#include <string.h>

 60#include "FIFOSampleBuffer.h"

 61#include "PeakFinder.h"

 62#include "BPMDetect.h"

 63
 64using namespace soundtouch;
 65
 66#define INPUT_BLOCK_SAMPLES       2048

 67#define DECIMATED_BLOCK_SAMPLES   256

 68
 69typedef unsigned short ushort;
 70
 71/// decay constant for calculating RMS volume sliding average approximation 
 72/// (time constant is about 10 sec)
 73const float avgdecay = 0.99986f;
 74
 75/// Normalization coefficient for calculating RMS sliding average approximation.
 76const float avgnorm = (1 - avgdecay);
 77
 78
 79
 80BPMDetect::BPMDetect(int numChannels, int sampleRate)
 81{
 82    xcorr = NULL;
 83
 84    buffer = new FIFOSampleBuffer();
 85
 86    decimateSum = 0;
 87    decimateCount = 0;
 88    decimateBy = 0;
 89
 90    this->sampleRate = sampleRate;
 91    this->channels = numChannels;
 92
 93    envelopeAccu = 0;
 94
 95    // Initialize RMS volume accumulator to RMS level of 3000 (out of 32768) that's
 96    // a typical RMS signal level value for song data. This value is then adapted
 97    // to the actual level during processing.
 98#ifdef INTEGER_SAMPLES

 99    // integer samples
100    RMSVolumeAccu = (3000 * 3000) / avgnorm;
101#else

102    // float samples, scaled to range [-1..+1[
103    RMSVolumeAccu = (0.092f * 0.092f) / avgnorm;
104#endif

105
106    init(numChannels, sampleRate);
107}
108
109
110
111BPMDetect::~BPMDetect()
112{
113    delete[] xcorr;
114    delete buffer;
115}
116
117
118/// low-pass filter & decimate to about 500 Hz. return number of outputted samples.
119///
120/// Decimation is used to remove the unnecessary frequencies and thus to reduce 
121/// the amount of data needed to be processed as calculating autocorrelation 
122/// function is a very-very heavy operation.
123///
124/// Anti-alias filtering is done simply by averaging the samples. This is really a 
125/// poor-man's anti-alias filtering, but it's not so critical in this kind of application
126/// (it'd also be difficult to design a high-quality filter with steep cut-off at very 
127/// narrow band)
128int BPMDetect::decimate(SAMPLETYPE *dest, const SAMPLETYPE *src, int numsamples)
129{
130    int count, outcount;
131    LONG_SAMPLETYPE out;
132
133    assert(decimateBy != 0);
134    outcount = 0;
135    for (count = 0; count < numsamples; count ++) 
136    {
137        decimateSum += src[count];
138
139        decimateCount ++;
140        if (decimateCount >= decimateBy) 
141        {
142            // Store every Nth sample only
143            out = (LONG_SAMPLETYPE)(decimateSum / decimateBy);
144            decimateSum = 0;
145            decimateCount = 0;
146#ifdef INTEGER_SAMPLES

147            // check ranges for sure (shouldn't actually be necessary)
148            if (out > 32767) 
149            {
150                out = 32767;
151            } 
152            else if (out < -32768) 
153            {
154                out = -32768;
155            }
156#endif // INTEGER_SAMPLES
157            dest[outcount] = (SAMPLETYPE)out;
158            outcount ++;
159        }
160    }
161    return outcount;
162}
163
164
165
166// Calculates autocorrelation function of the sample history buffer
167void BPMDetect::updateXCorr(int process_samples)
168{
169    int offs;
170    SAMPLETYPE *pBuffer;
171    
172    assert(buffer->numSamples() >= (uint)(process_samples + windowLen));
173
174    pBuffer = buffer->ptrBegin();
175    for (offs = windowStart; offs < windowLen; offs ++) 
176    {
177        LONG_SAMPLETYPE sum;
178        int i;
179
180        sum = 0;
181        for (i = 0; i < process_samples; i ++) 
182        {
183            sum += pBuffer[i] * pBuffer[i + offs];    // scaling the sub-result shouldn't be necessary
184        }
185//        xcorr[offs] *= xcorr_decay;   // decay 'xcorr' here with suitable coefficients 
186                                        // if it's desired that the system adapts automatically to
187                                        // various bpms, e.g. in processing continouos music stream.
188                                        // The 'xcorr_decay' should be a value that's smaller than but 
189                                        // close to one, and should also depend on 'process_samples' value.
190
191        xcorr[offs] += (float)sum;
192    }
193}
194
195
196
197// Calculates envelope of the sample data
198void BPMDetect::calcEnvelope(SAMPLETYPE *samples, int numsamples) 
199{
200    const float decay = 0.7f;               // decay constant for smoothing the envelope
201    const float norm = (1 - decay);
202
203    int i;
204    LONG_SAMPLETYPE out;
205    float val;
206
207    for (i = 0; i < numsamples; i ++) 
208    {
209        // calc average RMS volume
210        RMSVolumeAccu *= avgdecay;
211        val = (float)fabs((float)samples[i]);
212        RMSVolumeAccu += val * val;
213
214        // cut amplitudes that are below 2 times average RMS volume
215        // (we're interested in peak values, not the silent moments)
216        val -= 2 * (float)sqrt(RMSVolumeAccu * avgnorm);
217        val = (val > 0) ? val : 0;
218
219        // smooth amplitude envelope
220        envelopeAccu *= decay;
221        envelopeAccu += val;
222        out = (LONG_SAMPLETYPE)(envelopeAccu * norm);
223
224#ifdef INTEGER_SAMPLES

225        // cut peaks (shouldn't be necessary though)
226        if (out > 32767) out = 32767;
227#endif // INTEGER_SAMPLES
228        samples[i] = (SAMPLETYPE)out;
229    }
230}
231
232
233
234void BPMDetect::inputSamples(SAMPLETYPE *samples, int numSamples)
235{
236    SAMPLETYPE decimated[DECIMATED_BLOCK_SAMPLES];
237
238    // convert from stereo to mono if necessary
239    if (channels == 2)
240    {
241        int i;
242
243        for (i = 0; i < numSamples; i ++)
244        {
245            samples[i] = (samples[i * 2] + samples[i * 2 + 1]) / 2;
246        }
247    }
248    
249    // decimate
250    numSamples = decimate(decimated, samples, numSamples);
251
252    // envelope new samples and add them to buffer
253    calcEnvelope(decimated, numSamples);
254    buffer->putSamples(decimated, numSamples);
255
256    // when the buffer has enought samples for processing...
257    if ((int)buffer->numSamples() > windowLen) 
258    {
259        int processLength;
260
261        // how many samples are processed
262        processLength = buffer->numSamples() - windowLen;
263
264        // ... calculate autocorrelations for oldest samples...
265        updateXCorr(processLength);
266        // ... and remove them from the buffer
267        buffer->receiveSamples(processLength);
268    }
269}
270
271
272void BPMDetect::init(int numChannels, int sampleRate)
273{
274    this->sampleRate = sampleRate;
275
276    // choose decimation factor so that result is approx. 500 Hz
277    decimateBy = sampleRate / 500;
278    assert(decimateBy > 0);
279    assert(INPUT_BLOCK_SAMPLES < decimateBy * DECIMATED_BLOCK_SAMPLES);
280
281    // Calculate window length & starting item according to desired min & max bpms
282    windowLen = (60 * sampleRate) / (decimateBy * MIN_BPM);
283    windowStart = (60 * sampleRate) / (decimateBy * MAX_BPM);
284
285    assert(windowLen > windowStart);
286
287    // allocate new working objects
288    xcorr = new float[windowLen];
289    memset(xcorr, 0, windowLen * sizeof(float));
290
291    // we do processing in mono mode
292    buffer->setChannels(1);
293    buffer->clear();
294}
295
296
297
298float BPMDetect::getBpm()
299{
300    double peakPos;
301    PeakFinder peakFinder;
302
303    // find peak position
304    peakPos = peakFinder.detectPeak(xcorr, windowStart, windowLen);
305
306    assert(decimateBy != 0);
307    if (peakPos < 1e-6) return 0.0; // detection failed.
308
309    // calculate BPM
310    return (float)(60.0 * (((double)sampleRate / (double)decimateBy) / peakPos));
311}