annotate dsp/tempotracking/DownBeat.cpp @ 54:5bec06ecc88a

* First cut at Matthew's downbeat estimator -- untested so far
author cannam
date Tue, 10 Feb 2009 12:52:43 +0000
parents
children 7fe29d8a7eaf
rev   line source
cannam@54 1 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */
cannam@54 2
cannam@54 3 /*
cannam@54 4 QM DSP Library
cannam@54 5
cannam@54 6 Centre for Digital Music, Queen Mary, University of London.
cannam@54 7 This file copyright 2008-2009 Matthew Davies and QMUL.
cannam@54 8 All rights reserved.
cannam@54 9 */
cannam@54 10
cannam@54 11 #include "DownBeat.h"
cannam@54 12
cannam@54 13 #include "maths/MathAliases.h"
cannam@54 14 #include "maths/MathUtilities.h"
cannam@54 15 #include "dsp/transforms/FFT.h"
cannam@54 16
cannam@54 17 #include <iostream>
cannam@54 18 #include <cstdlib>
cannam@54 19
cannam@54 20 DownBeat::DownBeat(float originalSampleRate,
cannam@54 21 size_t decimationFactor,
cannam@54 22 size_t dfIncrement) :
cannam@54 23 m_rate(originalSampleRate),
cannam@54 24 m_factor(decimationFactor),
cannam@54 25 m_increment(dfIncrement),
cannam@54 26 m_decimator1(0),
cannam@54 27 m_decimator2(0),
cannam@54 28 m_buffer(0),
cannam@54 29 m_bufsiz(0),
cannam@54 30 m_buffill(0),
cannam@54 31 m_beatframesize(0),
cannam@54 32 m_beatframe(0)
cannam@54 33 {
cannam@54 34 // beat frame size is next power of two up from 1.3 seconds at the
cannam@54 35 // downsampled rate (happens to produce 4096 for 44100 or 48000 at
cannam@54 36 // 16x decimation, which is our expected normal situation)
cannam@54 37 int bfs = int((m_rate / decimationFactor) * 1.3);
cannam@54 38 m_beatframesize = 1;
cannam@54 39 while (bfs) { bfs >>= 1; m_beatframesize <<= 1; }
cannam@54 40 std::cerr << "rate = " << m_rate << ", bfs = " << m_beatframesize << std::endl;
cannam@54 41 m_beatframe = new double[m_beatframesize];
cannam@54 42 m_fftRealOut = new double[m_beatframesize];
cannam@54 43 m_fftImagOut = new double[m_beatframesize];
cannam@54 44 }
cannam@54 45
cannam@54 46 DownBeat::~DownBeat()
cannam@54 47 {
cannam@54 48 delete m_decimator1;
cannam@54 49 delete m_decimator2;
cannam@54 50 if (m_buffer) free(m_buffer);
cannam@54 51 delete[] m_decbuf;
cannam@54 52 delete[] m_beatframe;
cannam@54 53 delete[] m_fftRealOut;
cannam@54 54 delete[] m_fftImagOut;
cannam@54 55 }
cannam@54 56
cannam@54 57 void
cannam@54 58 DownBeat::makeDecimators()
cannam@54 59 {
cannam@54 60 if (m_factor < 2) return;
cannam@54 61 int highest = Decimator::getHighestSupportedFactor();
cannam@54 62 if (m_factor <= highest) {
cannam@54 63 m_decimator1 = new Decimator(m_increment, m_factor);
cannam@54 64 return;
cannam@54 65 }
cannam@54 66 m_decimator1 = new Decimator(m_increment, highest);
cannam@54 67 m_decimator2 = new Decimator(m_increment / highest, m_factor / highest);
cannam@54 68 m_decbuf = new double[m_factor / highest];
cannam@54 69 }
cannam@54 70
cannam@54 71 void
cannam@54 72 DownBeat::pushAudioBlock(const double *audio)
cannam@54 73 {
cannam@54 74 if (m_buffill + (m_increment / m_factor) > m_bufsiz) {
cannam@54 75 if (m_bufsiz == 0) m_bufsiz = m_increment * 16;
cannam@54 76 else m_bufsiz = m_bufsiz * 2;
cannam@54 77 if (!m_buffer) {
cannam@54 78 m_buffer = (double *)malloc(m_bufsiz * sizeof(double));
cannam@54 79 } else {
cannam@54 80 std::cerr << "DownBeat::pushAudioBlock: realloc m_buffer to " << m_bufsiz << std::endl;
cannam@54 81 m_buffer = (double *)realloc(m_buffer, m_bufsiz * sizeof(double));
cannam@54 82 }
cannam@54 83 }
cannam@54 84 if (!m_decimator1) makeDecimators();
cannam@54 85 if (m_decimator2) {
cannam@54 86 m_decimator1->process(audio, m_decbuf);
cannam@54 87 m_decimator2->process(m_decbuf, m_buffer + m_buffill);
cannam@54 88 } else {
cannam@54 89 m_decimator1->process(audio, m_buffer + m_buffill);
cannam@54 90 }
cannam@54 91 m_buffill += m_increment / m_factor;
cannam@54 92 }
cannam@54 93
cannam@54 94 const double *
cannam@54 95 DownBeat::getBufferedAudio(size_t &length) const
cannam@54 96 {
cannam@54 97 length = m_buffill;
cannam@54 98 return m_buffer;
cannam@54 99 }
cannam@54 100
cannam@54 101 void
cannam@54 102 DownBeat::findDownBeats(const double *audio,
cannam@54 103 size_t audioLength,
cannam@54 104 const d_vec_t &beats,
cannam@54 105 i_vec_t &downbeats)
cannam@54 106 {
cannam@54 107 // FIND DOWNBEATS BY PARTITIONING THE INPUT AUDIO FILE INTO BEAT SEGMENTS
cannam@54 108 // WHERE THE AUDIO FRAMES ARE DOWNSAMPLED BY A FACTOR OF 16 (fs ~= 2700Hz)
cannam@54 109 // THEN TAKING THE JENSEN-SHANNON DIVERGENCE BETWEEN BEAT SYNCHRONOUS SPECTRAL FRAMES
cannam@54 110
cannam@54 111 // IMPLEMENTATION (MOSTLY) FOLLOWS:
cannam@54 112 // DAVIES AND PLUMBLEY "A SPECTRAL DIFFERENCE APPROACH TO EXTRACTING DOWNBEATS IN MUSICAL AUDIO"
cannam@54 113 // EUSIPCO 2006, FLORENCE, ITALY
cannam@54 114
cannam@54 115 d_vec_t newspec(m_beatframesize / 2); // magnitude spectrum of current beat
cannam@54 116 d_vec_t oldspec(m_beatframesize / 2); // magnitude spectrum of previous beat
cannam@54 117 d_vec_t specdiff;
cannam@54 118
cannam@54 119 if (audioLength == 0) return;
cannam@54 120
cannam@54 121 for (size_t i = 0; i + 1 < beats.size(); ++i) {
cannam@54 122
cannam@54 123 // Copy the extents of the current beat from downsampled array
cannam@54 124 // into beat frame buffer
cannam@54 125
cannam@54 126 size_t beatstart = (beats[i] * m_increment) / m_factor;
cannam@54 127 size_t beatend = (beats[i] * m_increment) / m_factor;
cannam@54 128 if (beatend >= audioLength) beatend = audioLength - 1;
cannam@54 129 if (beatend < beatstart) beatend = beatstart;
cannam@54 130 size_t beatlen = beatend - beatstart;
cannam@54 131
cannam@54 132 // Also apply a Hanning window to the beat frame buffer, sized
cannam@54 133 // to the beat extents rather than the frame size. (Because
cannam@54 134 // the size varies, it's easier to do this by hand than use
cannam@54 135 // our Window abstraction.)
cannam@54 136
cannam@54 137 for (size_t j = 0; j < beatlen; ++j) {
cannam@54 138 double mul = 0.5 * (1.0 - cos(TWO_PI * (double(j) / double(beatlen))));
cannam@54 139 m_beatframe[j] = audio[beatstart + j] * mul;
cannam@54 140 }
cannam@54 141
cannam@54 142 for (size_t j = beatlen; j < m_beatframesize; ++j) {
cannam@54 143 m_beatframe[j] = 0.0;
cannam@54 144 }
cannam@54 145
cannam@54 146 // Now FFT beat frame
cannam@54 147
cannam@54 148 FFT::process(m_beatframesize, false,
cannam@54 149 m_beatframe, 0, m_fftRealOut, m_fftImagOut);
cannam@54 150
cannam@54 151 // Calculate magnitudes
cannam@54 152
cannam@54 153 for (size_t j = 0; j < m_beatframesize/2; ++j) {
cannam@54 154 newspec[j] = sqrt(m_fftRealOut[j] * m_fftRealOut[j] +
cannam@54 155 m_fftImagOut[j] * m_fftImagOut[j]);
cannam@54 156 }
cannam@54 157
cannam@54 158 // Preserve peaks by applying adaptive threshold
cannam@54 159
cannam@54 160 MathUtilities::adaptiveThreshold(newspec);
cannam@54 161
cannam@54 162 // Calculate JS divergence between new and old spectral frames
cannam@54 163
cannam@54 164 specdiff.push_back(measureSpecDiff(oldspec, newspec));
cannam@54 165
cannam@54 166 // Copy newspec across to old
cannam@54 167
cannam@54 168 for (size_t j = 0; j < m_beatframesize/2; ++j) {
cannam@54 169 oldspec[j] = newspec[j];
cannam@54 170 }
cannam@54 171 }
cannam@54 172
cannam@54 173 // We now have all spectral difference measures in specdiff
cannam@54 174
cannam@54 175 uint timesig = 4; // SHOULD REPLACE THIS WITH A FIND_METER FUNCTION - OR USER PARAMETER
cannam@54 176 d_vec_t dbcand(timesig); // downbeat candidates
cannam@54 177
cannam@54 178 // look for beat transition which leads to greatest spectral change
cannam@54 179 for (int beat = 0; beat < timesig; ++beat) {
cannam@54 180 for (int example = beat; example < specdiff.size(); ++example) {
cannam@54 181 dbcand[beat] += (specdiff[example]) / timesig;
cannam@54 182 }
cannam@54 183 }
cannam@54 184
cannam@54 185 // first downbeat is beat at index of maximum value of dbcand
cannam@54 186 int dbind = MathUtilities::getMax(dbcand);
cannam@54 187
cannam@54 188 // remaining downbeats are at timesig intervals from the first
cannam@54 189 for (int i = dbind; i < beats.size(); i += timesig) {
cannam@54 190 downbeats.push_back(i);
cannam@54 191 }
cannam@54 192 }
cannam@54 193
cannam@54 194 double
cannam@54 195 DownBeat::measureSpecDiff(d_vec_t oldspec, d_vec_t newspec)
cannam@54 196 {
cannam@54 197 // JENSEN-SHANNON DIVERGENCE BETWEEN SPECTRAL FRAMES
cannam@54 198
cannam@54 199 uint SPECSIZE = 512; // ONLY LOOK AT FIRST 512 SAMPLES OF SPECTRUM.
cannam@54 200 if (SPECSIZE > oldspec.size()/4) {
cannam@54 201 SPECSIZE = oldspec.size()/4;
cannam@54 202 }
cannam@54 203 double SD = 0.;
cannam@54 204 double sd1 = 0.;
cannam@54 205
cannam@54 206 double sumnew = 0.;
cannam@54 207 double sumold = 0.;
cannam@54 208
cannam@54 209 for (uint i = 0;i < SPECSIZE;i++)
cannam@54 210 {
cannam@54 211 newspec[i] +=EPS;
cannam@54 212 oldspec[i] +=EPS;
cannam@54 213
cannam@54 214 sumnew+=newspec[i];
cannam@54 215 sumold+=oldspec[i];
cannam@54 216 }
cannam@54 217
cannam@54 218 for (uint i = 0;i < SPECSIZE;i++)
cannam@54 219 {
cannam@54 220 newspec[i] /= (sumnew);
cannam@54 221 oldspec[i] /= (sumold);
cannam@54 222
cannam@54 223 // IF ANY SPECTRAL VALUES ARE 0 (SHOULDN'T BE ANY!) SET THEM TO 1
cannam@54 224 if (newspec[i] == 0)
cannam@54 225 {
cannam@54 226 newspec[i] = 1.;
cannam@54 227 }
cannam@54 228
cannam@54 229 if (oldspec[i] == 0)
cannam@54 230 {
cannam@54 231 oldspec[i] = 1.;
cannam@54 232 }
cannam@54 233
cannam@54 234 // JENSEN-SHANNON CALCULATION
cannam@54 235 sd1 = 0.5*oldspec[i] + 0.5*newspec[i];
cannam@54 236 SD = SD + (-sd1*log(sd1)) + (0.5*(oldspec[i]*log(oldspec[i]))) + (0.5*(newspec[i]*log(newspec[i])));
cannam@54 237 }
cannam@54 238
cannam@54 239 return SD;
cannam@54 240 }
cannam@54 241