cannam@54: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
cannam@54: 
cannam@54: /*
cannam@54:     QM DSP Library
cannam@54: 
cannam@54:     Centre for Digital Music, Queen Mary, University of London.
cannam@54:     This file copyright 2008-2009 Matthew Davies and QMUL.
cannam@54:     All rights reserved.
cannam@54: */
cannam@54: 
cannam@54: #include "DownBeat.h"
cannam@54: 
cannam@54: #include "maths/MathAliases.h"
cannam@54: #include "maths/MathUtilities.h"
cannam@55: #include "maths/KLDivergence.h"
cannam@54: #include "dsp/transforms/FFT.h"
cannam@54: 
cannam@54: #include <iostream>
cannam@54: #include <cstdlib>
cannam@54: 
cannam@54: DownBeat::DownBeat(float originalSampleRate,
cannam@54:                    size_t decimationFactor,
cannam@54:                    size_t dfIncrement) :
cannam@55:     m_bpb(0),
cannam@54:     m_rate(originalSampleRate),
cannam@54:     m_factor(decimationFactor),
cannam@54:     m_increment(dfIncrement),
cannam@54:     m_decimator1(0),
cannam@54:     m_decimator2(0),
cannam@54:     m_buffer(0),
cannam@54:     m_bufsiz(0),
cannam@54:     m_buffill(0),
cannam@54:     m_beatframesize(0),
cannam@54:     m_beatframe(0)
cannam@54: {
cannam@54:     // beat frame size is next power of two up from 1.3 seconds at the
cannam@54:     // downsampled rate (happens to produce 4096 for 44100 or 48000 at
cannam@54:     // 16x decimation, which is our expected normal situation)
cannam@55:     m_beatframesize = MathUtilities::nextPowerOfTwo
cannam@55:         (int((m_rate / decimationFactor) * 1.3));
cannam@54:     std::cerr << "rate = " << m_rate << ", bfs = " << m_beatframesize << std::endl;
cannam@54:     m_beatframe = new double[m_beatframesize];
cannam@54:     m_fftRealOut = new double[m_beatframesize];
cannam@54:     m_fftImagOut = new double[m_beatframesize];
cannam@54: }
cannam@54: 
cannam@54: DownBeat::~DownBeat()
cannam@54: {
cannam@54:     delete m_decimator1;
cannam@54:     delete m_decimator2;
cannam@54:     if (m_buffer) free(m_buffer);
cannam@54:     delete[] m_decbuf;
cannam@54:     delete[] m_beatframe;
cannam@54:     delete[] m_fftRealOut;
cannam@54:     delete[] m_fftImagOut;
cannam@54: }
cannam@54: 
cannam@54: void
cannam@55: DownBeat::setBeatsPerBar(int bpb)
cannam@55: {
cannam@55:     m_bpb = bpb;
cannam@55: }
cannam@55: 
cannam@55: void
cannam@54: DownBeat::makeDecimators()
cannam@54: {
cannam@54:     if (m_factor < 2) return;
cannam@54:     int highest = Decimator::getHighestSupportedFactor();
cannam@54:     if (m_factor <= highest) {
cannam@54:         m_decimator1 = new Decimator(m_increment, m_factor);
cannam@55:         std::cerr << "DownBeat: decimator 1 factor " << m_factor << ", size " << m_increment << std::endl;
cannam@54:         return;
cannam@54:     }
cannam@54:     m_decimator1 = new Decimator(m_increment, highest);
cannam@55:     std::cerr << "DownBeat: decimator 1 factor " << highest << ", size " << m_increment << std::endl;
cannam@54:     m_decimator2 = new Decimator(m_increment / highest, m_factor / highest);
cannam@55:     std::cerr << "DownBeat: decimator 2 factor " << m_factor / highest << ", size " << m_increment / highest << std::endl;
cannam@55:     m_decbuf = new float[m_increment / highest];
cannam@54: }
cannam@54: 
cannam@54: void
cannam@55: DownBeat::pushAudioBlock(const float *audio)
cannam@54: {
cannam@54:     if (m_buffill + (m_increment / m_factor) > m_bufsiz) {
cannam@54:         if (m_bufsiz == 0) m_bufsiz = m_increment * 16;
cannam@54:         else m_bufsiz = m_bufsiz * 2;
cannam@54:         if (!m_buffer) {
cannam@55:             m_buffer = (float *)malloc(m_bufsiz * sizeof(float));
cannam@54:         } else {
cannam@54:             std::cerr << "DownBeat::pushAudioBlock: realloc m_buffer to " << m_bufsiz << std::endl;
cannam@55:             m_buffer = (float *)realloc(m_buffer, m_bufsiz * sizeof(float));
cannam@54:         }
cannam@54:     }
cannam@54:     if (!m_decimator1) makeDecimators();
cannam@55:     float rmsin = 0, rmsout = 0;
cannam@55:     for (int i = 0; i < m_increment; ++i) {
cannam@55:         rmsin += audio[i] * audio[i];
cannam@55:     }
cannam@54:     if (m_decimator2) {
cannam@54:         m_decimator1->process(audio, m_decbuf);
cannam@54:         m_decimator2->process(m_decbuf, m_buffer + m_buffill);
cannam@54:     } else {
cannam@54:         m_decimator1->process(audio, m_buffer + m_buffill);
cannam@54:     }
cannam@55:     for (int i = 0; i < m_increment / m_factor; ++i) {
cannam@55:         rmsout += m_buffer[m_buffill + i] * m_buffer[m_buffill + i];
cannam@55:     }
cannam@55:     std::cerr << "pushAudioBlock: rms in " << sqrt(rmsin) << ", out " << sqrt(rmsout) << std::endl;
cannam@54:     m_buffill += m_increment / m_factor;
cannam@54: }
cannam@54:     
cannam@55: const float *
cannam@54: DownBeat::getBufferedAudio(size_t &length) const
cannam@54: {
cannam@54:     length = m_buffill;
cannam@54:     return m_buffer;
cannam@54: }
cannam@54: 
cannam@54: void
cannam@55: DownBeat::resetAudioBuffer()
cannam@55: {
cannam@55:     if (m_buffer) free(m_buffer);
cannam@55:     m_buffill = 0;
cannam@55:     m_bufsiz = 0;
cannam@55: }
cannam@55: 
cannam@55: void
cannam@55: DownBeat::findDownBeats(const float *audio,
cannam@54:                         size_t audioLength,
cannam@54:                         const d_vec_t &beats,
cannam@54:                         i_vec_t &downbeats)
cannam@54: {
cannam@54:     // FIND DOWNBEATS BY PARTITIONING THE INPUT AUDIO FILE INTO BEAT SEGMENTS
cannam@54:     // WHERE THE AUDIO FRAMES ARE DOWNSAMPLED  BY A FACTOR OF 16 (fs ~= 2700Hz)
cannam@54:     // THEN TAKING THE JENSEN-SHANNON DIVERGENCE BETWEEN BEAT SYNCHRONOUS SPECTRAL FRAMES
cannam@54: 
cannam@54:     // IMPLEMENTATION (MOSTLY) FOLLOWS:
cannam@54:     //  DAVIES AND PLUMBLEY "A SPECTRAL DIFFERENCE APPROACH TO EXTRACTING DOWNBEATS IN MUSICAL AUDIO"
cannam@54:     //  EUSIPCO 2006, FLORENCE, ITALY
cannam@54: 
cannam@54:     d_vec_t newspec(m_beatframesize / 2); // magnitude spectrum of current beat
cannam@54:     d_vec_t oldspec(m_beatframesize / 2); // magnitude spectrum of previous beat
cannam@56: 
cannam@56:     m_beatsd.clear();
cannam@54: 
cannam@54:     if (audioLength == 0) return;
cannam@54: 
cannam@54:     for (size_t i = 0; i + 1 < beats.size(); ++i) {
cannam@54: 
cannam@54:         // Copy the extents of the current beat from downsampled array
cannam@54:         // into beat frame buffer
cannam@54: 
cannam@54:         size_t beatstart = (beats[i] * m_increment) / m_factor;
cannam@55:         size_t beatend = (beats[i+1] * m_increment) / m_factor;
cannam@54:         if (beatend >= audioLength) beatend = audioLength - 1;
cannam@54:         if (beatend < beatstart) beatend = beatstart;
cannam@54:         size_t beatlen = beatend - beatstart;
cannam@54: 
cannam@54:         // Also apply a Hanning window to the beat frame buffer, sized
cannam@54:         // to the beat extents rather than the frame size.  (Because
cannam@54:         // the size varies, it's easier to do this by hand than use
cannam@54:         // our Window abstraction.)
cannam@54: 
cannam@55:         float rms = 0;
cannam@54:         for (size_t j = 0; j < beatlen; ++j) {
cannam@54:             double mul = 0.5 * (1.0 - cos(TWO_PI * (double(j) / double(beatlen))));
cannam@54:             m_beatframe[j] = audio[beatstart + j] * mul;
cannam@55:             rms += m_beatframe[j] * m_beatframe[j];
cannam@54:         }
cannam@55:         rms = sqrt(rms);
cannam@55:         std::cerr << "beat " << i << ": audio rms " << rms << std::endl;
cannam@54: 
cannam@54:         for (size_t j = beatlen; j < m_beatframesize; ++j) {
cannam@54:             m_beatframe[j] = 0.0;
cannam@54:         }
cannam@54: 
cannam@54:         // Now FFT beat frame
cannam@54:         
cannam@54:         FFT::process(m_beatframesize, false,
cannam@54:                      m_beatframe, 0, m_fftRealOut, m_fftImagOut);
cannam@54:         
cannam@54:         // Calculate magnitudes
cannam@54: 
cannam@54:         for (size_t j = 0; j < m_beatframesize/2; ++j) {
cannam@54:             newspec[j] = sqrt(m_fftRealOut[j] * m_fftRealOut[j] +
cannam@54:                               m_fftImagOut[j] * m_fftImagOut[j]);
cannam@54:         }
cannam@54: 
cannam@54:         // Preserve peaks by applying adaptive threshold
cannam@54: 
cannam@54:         MathUtilities::adaptiveThreshold(newspec);
cannam@54: 
cannam@54:         // Calculate JS divergence between new and old spectral frames
cannam@54: 
cannam@56:         if (i > 0) { // otherwise we have no previous frame
cannam@56:             m_beatsd.push_back(measureSpecDiff(oldspec, newspec));
cannam@56:             std::cerr << "specdiff: " << m_beatsd[m_beatsd.size()-1] << std::endl;
cannam@56:         }
cannam@54: 
cannam@54:         // Copy newspec across to old
cannam@54: 
cannam@54:         for (size_t j = 0; j < m_beatframesize/2; ++j) {
cannam@54:             oldspec[j] = newspec[j];
cannam@54:         }
cannam@54:     }
cannam@54: 
cannam@54:     // We now have all spectral difference measures in specdiff
cannam@54: 
cannam@55:     uint timesig = m_bpb;
cannam@55:     if (timesig == 0) timesig = 4;
cannam@55: 
cannam@54:     d_vec_t dbcand(timesig); // downbeat candidates
cannam@54: 
cannam@55:     for (int beat = 0; beat < timesig; ++beat) {
cannam@55:         dbcand[beat] = 0;
cannam@55:     }
cannam@55: 
cannam@54:     // look for beat transition which leads to greatest spectral change
cannam@54:     for (int beat = 0; beat < timesig; ++beat) {
cannam@56:         int count = 0;
cannam@56:         for (int example = beat - 1; example < m_beatsd.size(); example += timesig) {
cannam@56:             if (example < 0) continue;
cannam@56:             dbcand[beat] += (m_beatsd[example]) / timesig;
cannam@56:             ++count;
cannam@54:         }
cannam@56:         if (count > 0) m_beatsd[beat] /= count;
cannam@55:         std::cerr << "dbcand[" << beat << "] = " << dbcand[beat] << std::endl;
cannam@54:     }
cannam@54: 
cannam@55: 
cannam@54:     // first downbeat is beat at index of maximum value of dbcand
cannam@54:     int dbind = MathUtilities::getMax(dbcand);
cannam@54: 
cannam@54:     // remaining downbeats are at timesig intervals from the first
cannam@54:     for (int i = dbind; i < beats.size(); i += timesig) {
cannam@54:         downbeats.push_back(i);
cannam@54:     }
cannam@54: }
cannam@54: 
cannam@54: double
cannam@54: DownBeat::measureSpecDiff(d_vec_t oldspec, d_vec_t newspec)
cannam@54: {
cannam@54:     // JENSEN-SHANNON DIVERGENCE BETWEEN SPECTRAL FRAMES
cannam@54: 
cannam@54:     uint SPECSIZE = 512;   // ONLY LOOK AT FIRST 512 SAMPLES OF SPECTRUM. 
cannam@54:     if (SPECSIZE > oldspec.size()/4) {
cannam@54:         SPECSIZE = oldspec.size()/4;
cannam@54:     }
cannam@54:     double SD = 0.;
cannam@54:     double sd1 = 0.;
cannam@54: 
cannam@54:     double sumnew = 0.;
cannam@54:     double sumold = 0.;
cannam@54:   
cannam@54:     for (uint i = 0;i < SPECSIZE;i++)
cannam@54:     {
cannam@54:         newspec[i] +=EPS;
cannam@54:         oldspec[i] +=EPS;
cannam@54:         
cannam@54:         sumnew+=newspec[i];
cannam@54:         sumold+=oldspec[i];
cannam@54:     } 
cannam@54:     
cannam@54:     for (uint i = 0;i < SPECSIZE;i++)
cannam@54:     {
cannam@54:         newspec[i] /= (sumnew);
cannam@54:         oldspec[i] /= (sumold);
cannam@54:         
cannam@54:         // IF ANY SPECTRAL VALUES ARE 0 (SHOULDN'T BE ANY!) SET THEM TO 1
cannam@54:         if (newspec[i] == 0)
cannam@54:         {
cannam@54:             newspec[i] = 1.;
cannam@54:         }
cannam@54:         
cannam@54:         if (oldspec[i] == 0)
cannam@54:         {
cannam@54:             oldspec[i] = 1.;
cannam@54:         }
cannam@54:         
cannam@54:         // JENSEN-SHANNON CALCULATION
cannam@54:         sd1 = 0.5*oldspec[i] + 0.5*newspec[i];	
cannam@54:         SD = SD + (-sd1*log(sd1)) + (0.5*(oldspec[i]*log(oldspec[i]))) + (0.5*(newspec[i]*log(newspec[i])));
cannam@54:     }
cannam@54:     
cannam@54:     return SD;
cannam@54: }
cannam@54: 
cannam@56: void
cannam@56: DownBeat::getBeatSD(vector<double> &beatsd) const
cannam@56: {
cannam@56:     for (int i = 0; i < m_beatsd.size(); ++i) beatsd.push_back(m_beatsd[i]);
cannam@56: }
cannam@56: