qm-dsp  1.8
DownBeat.cpp
Go to the documentation of this file.
00001 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
00002 
00003 /*
00004     QM DSP Library
00005 
00006     Centre for Digital Music, Queen Mary, University of London.
00007     This file copyright 2008-2009 Matthew Davies and QMUL.
00008 
00009     This program is free software; you can redistribute it and/or
00010     modify it under the terms of the GNU General Public License as
00011     published by the Free Software Foundation; either version 2 of the
00012     License, or (at your option) any later version.  See the file
00013     COPYING included with this distribution for more information.
00014 */
00015 
00016 #include "DownBeat.h"
00017 
00018 #include "maths/MathAliases.h"
00019 #include "maths/MathUtilities.h"
00020 #include "maths/KLDivergence.h"
00021 #include "dsp/transforms/FFT.h"
00022 
00023 #include <iostream>
00024 #include <cstdlib>
00025 
00026 DownBeat::DownBeat(float originalSampleRate,
00027                    size_t decimationFactor,
00028                    size_t dfIncrement) :
00029     m_bpb(0),
00030     m_rate(originalSampleRate),
00031     m_factor(decimationFactor),
00032     m_increment(dfIncrement),
00033     m_decimator1(0),
00034     m_decimator2(0),
00035     m_buffer(0),
00036     m_decbuf(0),
00037     m_bufsiz(0),
00038     m_buffill(0),
00039     m_beatframesize(0),
00040     m_beatframe(0)
00041 {
00042     // beat frame size is next power of two up from 1.3 seconds at the
00043     // downsampled rate (happens to produce 4096 for 44100 or 48000 at
00044     // 16x decimation, which is our expected normal situation)
00045     m_beatframesize = MathUtilities::nextPowerOfTwo
00046         (int((m_rate / decimationFactor) * 1.3));
00047     if (m_beatframesize < 2) {
00048         m_beatframesize = 2;
00049     }
00050 //    std::cerr << "rate = " << m_rate << ", dec = " << decimationFactor << ", bfs = " << m_beatframesize << std::endl;
00051     m_beatframe = new double[m_beatframesize];
00052     m_fftRealOut = new double[m_beatframesize];
00053     m_fftImagOut = new double[m_beatframesize];
00054     m_fft = new FFTReal(m_beatframesize);
00055 }
00056 
00057 DownBeat::~DownBeat()
00058 {
00059     delete m_decimator1;
00060     delete m_decimator2;
00061     if (m_buffer) free(m_buffer);
00062     delete[] m_decbuf;
00063     delete[] m_beatframe;
00064     delete[] m_fftRealOut;
00065     delete[] m_fftImagOut;
00066     delete m_fft;
00067 }
00068 
00069 void
00070 DownBeat::setBeatsPerBar(int bpb)
00071 {
00072     m_bpb = bpb;
00073 }
00074 
00075 void
00076 DownBeat::makeDecimators()
00077 {
00078 //    std::cerr << "m_factor = " << m_factor << std::endl;
00079     if (m_factor < 2) return;
00080     size_t highest = Decimator::getHighestSupportedFactor();
00081     if (m_factor <= highest) {
00082         m_decimator1 = new Decimator(m_increment, m_factor);
00083 //        std::cerr << "DownBeat: decimator 1 factor " << m_factor << ", size " << m_increment << std::endl;
00084         return;
00085     }
00086     m_decimator1 = new Decimator(m_increment, highest);
00087 //    std::cerr << "DownBeat: decimator 1 factor " << highest << ", size " << m_increment << std::endl;
00088     m_decimator2 = new Decimator(m_increment / highest, m_factor / highest);
00089 //    std::cerr << "DownBeat: decimator 2 factor " << m_factor / highest << ", size " << m_increment / highest << std::endl;
00090     m_decbuf = new float[m_increment / highest];
00091 }
00092 
00093 void
00094 DownBeat::pushAudioBlock(const float *audio)
00095 {
00096     if (m_buffill + (m_increment / m_factor) > m_bufsiz) {
00097         if (m_bufsiz == 0) m_bufsiz = m_increment * 16;
00098         else m_bufsiz = m_bufsiz * 2;
00099         if (!m_buffer) {
00100             m_buffer = (float *)malloc(m_bufsiz * sizeof(float));
00101         } else {
00102 //            std::cerr << "DownBeat::pushAudioBlock: realloc m_buffer to " << m_bufsiz << std::endl;
00103             m_buffer = (float *)realloc(m_buffer, m_bufsiz * sizeof(float));
00104         }
00105     }
00106     if (!m_decimator1 && m_factor > 1) makeDecimators();
00107 //    float rmsin = 0, rmsout = 0;
00108 //    for (int i = 0; i < m_increment; ++i) {
00109 //        rmsin += audio[i] * audio[i];
00110 //    }
00111     if (m_decimator2) {
00112         m_decimator1->process(audio, m_decbuf);
00113         m_decimator2->process(m_decbuf, m_buffer + m_buffill);
00114     } else if (m_decimator1) {
00115         m_decimator1->process(audio, m_buffer + m_buffill);
00116     } else {
00117         // just copy across (m_factor is presumably 1)
00118         for (size_t i = 0; i < m_increment; ++i) {
00119             (m_buffer + m_buffill)[i] = audio[i];
00120         }
00121     }
00122 //    for (int i = 0; i < m_increment / m_factor; ++i) {
00123 //        rmsout += m_buffer[m_buffill + i] * m_buffer[m_buffill + i];
00124 //    }
00125 //    std::cerr << "pushAudioBlock: rms in " << sqrt(rmsin) << ", out " << sqrt(rmsout) << std::endl;
00126     m_buffill += m_increment / m_factor;
00127 }
00128     
00129 const float *
00130 DownBeat::getBufferedAudio(size_t &length) const
00131 {
00132     length = m_buffill;
00133     return m_buffer;
00134 }
00135 
00136 void
00137 DownBeat::resetAudioBuffer()
00138 {
00139     if (m_buffer) free(m_buffer);
00140     m_buffer = 0;
00141     m_buffill = 0;
00142     m_bufsiz = 0;
00143 }
00144 
00145 void
00146 DownBeat::findDownBeats(const float *audio,
00147                         size_t audioLength,
00148                         const d_vec_t &beats,
00149                         i_vec_t &downbeats)
00150 {
00151     // FIND DOWNBEATS BY PARTITIONING THE INPUT AUDIO FILE INTO BEAT SEGMENTS
00152     // WHERE THE AUDIO FRAMES ARE DOWNSAMPLED  BY A FACTOR OF 16 (fs ~= 2700Hz)
00153     // THEN TAKING THE JENSEN-SHANNON DIVERGENCE BETWEEN BEAT SYNCHRONOUS SPECTRAL FRAMES
00154 
00155     // IMPLEMENTATION (MOSTLY) FOLLOWS:
00156     //  DAVIES AND PLUMBLEY "A SPECTRAL DIFFERENCE APPROACH TO EXTRACTING DOWNBEATS IN MUSICAL AUDIO"
00157     //  EUSIPCO 2006, FLORENCE, ITALY
00158 
00159     d_vec_t newspec(m_beatframesize / 2); // magnitude spectrum of current beat
00160     d_vec_t oldspec(m_beatframesize / 2); // magnitude spectrum of previous beat
00161 
00162     m_beatsd.clear();
00163 
00164     if (audioLength == 0) return;
00165 
00166     for (size_t i = 0; i + 1 < beats.size(); ++i) {
00167 
00168         // Copy the extents of the current beat from downsampled array
00169         // into beat frame buffer
00170 
00171         size_t beatstart = (beats[i] * m_increment) / m_factor;
00172         size_t beatend = (beats[i+1] * m_increment) / m_factor;
00173         if (beatend >= audioLength) beatend = audioLength - 1;
00174         if (beatend < beatstart) beatend = beatstart;
00175         size_t beatlen = beatend - beatstart;
00176 
00177         // Also apply a Hanning window to the beat frame buffer, sized
00178         // to the beat extents rather than the frame size.  (Because
00179         // the size varies, it's easier to do this by hand than use
00180         // our Window abstraction.)
00181 
00182 //        std::cerr << "beatlen = " << beatlen << std::endl;
00183 
00184 //        float rms = 0;
00185         for (size_t j = 0; j < beatlen && j < m_beatframesize; ++j) {
00186             double mul = 0.5 * (1.0 - cos(TWO_PI * (double(j) / double(beatlen))));
00187             m_beatframe[j] = audio[beatstart + j] * mul;
00188 //            rms += m_beatframe[j] * m_beatframe[j];
00189         }
00190 //        rms = sqrt(rms);
00191 //        std::cerr << "beat " << i << ": audio rms " << rms << std::endl;
00192 
00193         for (size_t j = beatlen; j < m_beatframesize; ++j) {
00194             m_beatframe[j] = 0.0;
00195         }
00196 
00197         // Now FFT beat frame
00198         
00199         m_fft->forward(m_beatframe, m_fftRealOut, m_fftImagOut);
00200         
00201         // Calculate magnitudes
00202 
00203         for (size_t j = 0; j < m_beatframesize/2; ++j) {
00204             newspec[j] = sqrt(m_fftRealOut[j] * m_fftRealOut[j] +
00205                               m_fftImagOut[j] * m_fftImagOut[j]);
00206         }
00207 
00208         // Preserve peaks by applying adaptive threshold
00209 
00210         MathUtilities::adaptiveThreshold(newspec);
00211 
00212         // Calculate JS divergence between new and old spectral frames
00213 
00214         if (i > 0) { // otherwise we have no previous frame
00215             m_beatsd.push_back(measureSpecDiff(oldspec, newspec));
00216 //            std::cerr << "specdiff: " << m_beatsd[m_beatsd.size()-1] << std::endl;
00217         }
00218 
00219         // Copy newspec across to old
00220 
00221         for (size_t j = 0; j < m_beatframesize/2; ++j) {
00222             oldspec[j] = newspec[j];
00223         }
00224     }
00225 
00226     // We now have all spectral difference measures in specdiff
00227 
00228     int timesig = m_bpb;
00229     if (timesig == 0) timesig = 4;
00230 
00231     d_vec_t dbcand(timesig); // downbeat candidates
00232 
00233     for (int beat = 0; beat < timesig; ++beat) {
00234         dbcand[beat] = 0;
00235     }
00236 
00237    // look for beat transition which leads to greatest spectral change
00238    for (int beat = 0; beat < timesig; ++beat) {
00239        int count = 0;
00240        for (int example = beat-1; example < (int)m_beatsd.size(); example += timesig) {
00241            if (example < 0) continue;
00242            dbcand[beat] += (m_beatsd[example]) / timesig;
00243            ++count;
00244        }
00245        if (count > 0) dbcand[beat] /= count;
00246 //        std::cerr << "dbcand[" << beat << "] = " << dbcand[beat] << std::endl;
00247    }
00248 
00249     // first downbeat is beat at index of maximum value of dbcand
00250     int dbind = MathUtilities::getMax(dbcand);
00251 
00252     // remaining downbeats are at timesig intervals from the first
00253     for (int i = dbind; i < (int)beats.size(); i += timesig) {
00254         downbeats.push_back(i);
00255     }
00256 }
00257 
00258 double
00259 DownBeat::measureSpecDiff(d_vec_t oldspec, d_vec_t newspec)
00260 {
00261     // JENSEN-SHANNON DIVERGENCE BETWEEN SPECTRAL FRAMES
00262 
00263     unsigned int SPECSIZE = 512;   // ONLY LOOK AT FIRST 512 SAMPLES OF SPECTRUM. 
00264     if (SPECSIZE > oldspec.size()/4) {
00265         SPECSIZE = oldspec.size()/4;
00266     }
00267     double SD = 0.;
00268     double sd1 = 0.;
00269 
00270     double sumnew = 0.;
00271     double sumold = 0.;
00272   
00273     for (unsigned int i = 0;i < SPECSIZE;i++)
00274     {
00275         newspec[i] +=EPS;
00276         oldspec[i] +=EPS;
00277         
00278         sumnew+=newspec[i];
00279         sumold+=oldspec[i];
00280     } 
00281     
00282     for (unsigned int i = 0;i < SPECSIZE;i++)
00283     {
00284         newspec[i] /= (sumnew);
00285         oldspec[i] /= (sumold);
00286         
00287         // IF ANY SPECTRAL VALUES ARE 0 (SHOULDN'T BE ANY!) SET THEM TO 1
00288         if (newspec[i] == 0)
00289         {
00290             newspec[i] = 1.;
00291         }
00292         
00293         if (oldspec[i] == 0)
00294         {
00295             oldspec[i] = 1.;
00296         }
00297         
00298         // JENSEN-SHANNON CALCULATION
00299         sd1 = 0.5*oldspec[i] + 0.5*newspec[i];  
00300         SD = SD + (-sd1*log(sd1)) + (0.5*(oldspec[i]*log(oldspec[i]))) + (0.5*(newspec[i]*log(newspec[i])));
00301     }
00302     
00303     return SD;
00304 }
00305 
00306 void
00307 DownBeat::getBeatSD(vector<double> &beatsd) const
00308 {
00309     for (int i = 0; i < (int)m_beatsd.size(); ++i) beatsd.push_back(m_beatsd[i]);
00310 }
00311