qm-dsp  1.8
DownBeat.cpp
Go to the documentation of this file.
1 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */
2 
3 /*
4  QM DSP Library
5 
6  Centre for Digital Music, Queen Mary, University of London.
7  This file copyright 2008-2009 Matthew Davies and QMUL.
8 
9  This program is free software; you can redistribute it and/or
10  modify it under the terms of the GNU General Public License as
11  published by the Free Software Foundation; either version 2 of the
12  License, or (at your option) any later version. See the file
13  COPYING included with this distribution for more information.
14 */
15 
16 #include "DownBeat.h"
17 
18 #include "maths/MathAliases.h"
19 #include "maths/MathUtilities.h"
20 #include "maths/KLDivergence.h"
21 #include "dsp/transforms/FFT.h"
22 
23 #include <iostream>
24 #include <cstdlib>
25 
26 DownBeat::DownBeat(float originalSampleRate,
27  size_t decimationFactor,
28  size_t dfIncrement) :
29  m_bpb(0),
30  m_rate(originalSampleRate),
31  m_factor(decimationFactor),
32  m_increment(dfIncrement),
33  m_decimator1(0),
34  m_decimator2(0),
35  m_buffer(0),
36  m_decbuf(0),
37  m_bufsiz(0),
38  m_buffill(0),
39  m_beatframesize(0),
40  m_beatframe(0)
41 {
42  // beat frame size is next power of two up from 1.3 seconds at the
43  // downsampled rate (happens to produce 4096 for 44100 or 48000 at
44  // 16x decimation, which is our expected normal situation)
46  (int((m_rate / decimationFactor) * 1.3));
47  if (m_beatframesize < 2) {
48  m_beatframesize = 2;
49  }
50 // std::cerr << "rate = " << m_rate << ", dec = " << decimationFactor << ", bfs = " << m_beatframesize << std::endl;
51  m_beatframe = new double[m_beatframesize];
52  m_fftRealOut = new double[m_beatframesize];
53  m_fftImagOut = new double[m_beatframesize];
55 }
56 
58 {
59  delete m_decimator1;
60  delete m_decimator2;
61  if (m_buffer) free(m_buffer);
62  delete[] m_decbuf;
63  delete[] m_beatframe;
64  delete[] m_fftRealOut;
65  delete[] m_fftImagOut;
66  delete m_fft;
67 }
68 
69 void
71 {
72  m_bpb = bpb;
73 }
74 
75 void
77 {
78 // std::cerr << "m_factor = " << m_factor << std::endl;
79  if (m_factor < 2) return;
80  size_t highest = Decimator::getHighestSupportedFactor();
81  if (m_factor <= highest) {
83 // std::cerr << "DownBeat: decimator 1 factor " << m_factor << ", size " << m_increment << std::endl;
84  return;
85  }
86  m_decimator1 = new Decimator(m_increment, highest);
87 // std::cerr << "DownBeat: decimator 1 factor " << highest << ", size " << m_increment << std::endl;
88  m_decimator2 = new Decimator(m_increment / highest, m_factor / highest);
89 // std::cerr << "DownBeat: decimator 2 factor " << m_factor / highest << ", size " << m_increment / highest << std::endl;
90  m_decbuf = new float[m_increment / highest];
91 }
92 
93 void
94 DownBeat::pushAudioBlock(const float *audio)
95 {
96  if (m_buffill + (m_increment / m_factor) > m_bufsiz) {
97  if (m_bufsiz == 0) m_bufsiz = m_increment * 16;
98  else m_bufsiz = m_bufsiz * 2;
99  if (!m_buffer) {
100  m_buffer = (float *)malloc(m_bufsiz * sizeof(float));
101  } else {
102 // std::cerr << "DownBeat::pushAudioBlock: realloc m_buffer to " << m_bufsiz << std::endl;
103  m_buffer = (float *)realloc(m_buffer, m_bufsiz * sizeof(float));
104  }
105  }
106  if (!m_decimator1 && m_factor > 1) makeDecimators();
107 // float rmsin = 0, rmsout = 0;
108 // for (int i = 0; i < m_increment; ++i) {
109 // rmsin += audio[i] * audio[i];
110 // }
111  if (m_decimator2) {
112  m_decimator1->process(audio, m_decbuf);
114  } else if (m_decimator1) {
116  } else {
117  // just copy across (m_factor is presumably 1)
118  for (size_t i = 0; i < m_increment; ++i) {
119  (m_buffer + m_buffill)[i] = audio[i];
120  }
121  }
122 // for (int i = 0; i < m_increment / m_factor; ++i) {
123 // rmsout += m_buffer[m_buffill + i] * m_buffer[m_buffill + i];
124 // }
125 // std::cerr << "pushAudioBlock: rms in " << sqrt(rmsin) << ", out " << sqrt(rmsout) << std::endl;
127 }
128 
129 const float *
130 DownBeat::getBufferedAudio(size_t &length) const
131 {
132  length = m_buffill;
133  return m_buffer;
134 }
135 
136 void
138 {
139  if (m_buffer) free(m_buffer);
140  m_buffer = 0;
141  m_buffill = 0;
142  m_bufsiz = 0;
143 }
144 
145 void
146 DownBeat::findDownBeats(const float *audio,
147  size_t audioLength,
148  const d_vec_t &beats,
149  i_vec_t &downbeats)
150 {
151  // FIND DOWNBEATS BY PARTITIONING THE INPUT AUDIO FILE INTO BEAT SEGMENTS
152  // WHERE THE AUDIO FRAMES ARE DOWNSAMPLED BY A FACTOR OF 16 (fs ~= 2700Hz)
153  // THEN TAKING THE JENSEN-SHANNON DIVERGENCE BETWEEN BEAT SYNCHRONOUS SPECTRAL FRAMES
154 
155  // IMPLEMENTATION (MOSTLY) FOLLOWS:
156  // DAVIES AND PLUMBLEY "A SPECTRAL DIFFERENCE APPROACH TO EXTRACTING DOWNBEATS IN MUSICAL AUDIO"
157  // EUSIPCO 2006, FLORENCE, ITALY
158 
159  d_vec_t newspec(m_beatframesize / 2); // magnitude spectrum of current beat
160  d_vec_t oldspec(m_beatframesize / 2); // magnitude spectrum of previous beat
161 
162  m_beatsd.clear();
163 
164  if (audioLength == 0) return;
165 
166  for (size_t i = 0; i + 1 < beats.size(); ++i) {
167 
168  // Copy the extents of the current beat from downsampled array
169  // into beat frame buffer
170 
171  size_t beatstart = (beats[i] * m_increment) / m_factor;
172  size_t beatend = (beats[i+1] * m_increment) / m_factor;
173  if (beatend >= audioLength) beatend = audioLength - 1;
174  if (beatend < beatstart) beatend = beatstart;
175  size_t beatlen = beatend - beatstart;
176 
177  // Also apply a Hanning window to the beat frame buffer, sized
178  // to the beat extents rather than the frame size. (Because
179  // the size varies, it's easier to do this by hand than use
180  // our Window abstraction.)
181 
182 // std::cerr << "beatlen = " << beatlen << std::endl;
183 
184 // float rms = 0;
185  for (size_t j = 0; j < beatlen && j < m_beatframesize; ++j) {
186  double mul = 0.5 * (1.0 - cos(TWO_PI * (double(j) / double(beatlen))));
187  m_beatframe[j] = audio[beatstart + j] * mul;
188 // rms += m_beatframe[j] * m_beatframe[j];
189  }
190 // rms = sqrt(rms);
191 // std::cerr << "beat " << i << ": audio rms " << rms << std::endl;
192 
193  for (size_t j = beatlen; j < m_beatframesize; ++j) {
194  m_beatframe[j] = 0.0;
195  }
196 
197  // Now FFT beat frame
198 
200 
201  // Calculate magnitudes
202 
203  for (size_t j = 0; j < m_beatframesize/2; ++j) {
204  newspec[j] = sqrt(m_fftRealOut[j] * m_fftRealOut[j] +
205  m_fftImagOut[j] * m_fftImagOut[j]);
206  }
207 
208  // Preserve peaks by applying adaptive threshold
209 
211 
212  // Calculate JS divergence between new and old spectral frames
213 
214  if (i > 0) { // otherwise we have no previous frame
215  m_beatsd.push_back(measureSpecDiff(oldspec, newspec));
216 // std::cerr << "specdiff: " << m_beatsd[m_beatsd.size()-1] << std::endl;
217  }
218 
219  // Copy newspec across to old
220 
221  for (size_t j = 0; j < m_beatframesize/2; ++j) {
222  oldspec[j] = newspec[j];
223  }
224  }
225 
226  // We now have all spectral difference measures in specdiff
227 
228  int timesig = m_bpb;
229  if (timesig == 0) timesig = 4;
230 
231  d_vec_t dbcand(timesig); // downbeat candidates
232 
233  for (int beat = 0; beat < timesig; ++beat) {
234  dbcand[beat] = 0;
235  }
236 
237  // look for beat transition which leads to greatest spectral change
238  for (int beat = 0; beat < timesig; ++beat) {
239  int count = 0;
240  for (int example = beat-1; example < (int)m_beatsd.size(); example += timesig) {
241  if (example < 0) continue;
242  dbcand[beat] += (m_beatsd[example]) / timesig;
243  ++count;
244  }
245  if (count > 0) dbcand[beat] /= count;
246 // std::cerr << "dbcand[" << beat << "] = " << dbcand[beat] << std::endl;
247  }
248 
249  // first downbeat is beat at index of maximum value of dbcand
250  int dbind = MathUtilities::getMax(dbcand);
251 
252  // remaining downbeats are at timesig intervals from the first
253  for (int i = dbind; i < (int)beats.size(); i += timesig) {
254  downbeats.push_back(i);
255  }
256 }
257 
258 double
260 {
261  // JENSEN-SHANNON DIVERGENCE BETWEEN SPECTRAL FRAMES
262 
263  unsigned int SPECSIZE = 512; // ONLY LOOK AT FIRST 512 SAMPLES OF SPECTRUM.
264  if (SPECSIZE > oldspec.size()/4) {
265  SPECSIZE = oldspec.size()/4;
266  }
267  double SD = 0.;
268  double sd1 = 0.;
269 
270  double sumnew = 0.;
271  double sumold = 0.;
272 
273  for (unsigned int i = 0;i < SPECSIZE;i++)
274  {
275  newspec[i] +=EPS;
276  oldspec[i] +=EPS;
277 
278  sumnew+=newspec[i];
279  sumold+=oldspec[i];
280  }
281 
282  for (unsigned int i = 0;i < SPECSIZE;i++)
283  {
284  newspec[i] /= (sumnew);
285  oldspec[i] /= (sumold);
286 
287  // IF ANY SPECTRAL VALUES ARE 0 (SHOULDN'T BE ANY!) SET THEM TO 1
288  if (newspec[i] == 0)
289  {
290  newspec[i] = 1.;
291  }
292 
293  if (oldspec[i] == 0)
294  {
295  oldspec[i] = 1.;
296  }
297 
298  // JENSEN-SHANNON CALCULATION
299  sd1 = 0.5*oldspec[i] + 0.5*newspec[i];
300  SD = SD + (-sd1*log(sd1)) + (0.5*(oldspec[i]*log(oldspec[i]))) + (0.5*(newspec[i]*log(newspec[i])));
301  }
302 
303  return SD;
304 }
305 
306 void
307 DownBeat::getBeatSD(vector<double> &beatsd) const
308 {
309  for (int i = 0; i < (int)m_beatsd.size(); ++i) beatsd.push_back(m_beatsd[i]);
310 }
311 
void process(const double *src, double *dst)
Process inLength samples (as supplied to constructor) from src and write inLength / decFactor samples...
Definition: Decimator.cpp:200
const float * getBufferedAudio(size_t &length) const
Retrieve the accumulated audio produced by pushAudioBlock calls.
Definition: DownBeat.cpp:130
vector< int > i_vec_t
Definition: DownBeat.h:110
static void adaptiveThreshold(std::vector< double > &data)
Threshold the input/output vector data against a moving-mean average filter.
double * m_fftRealOut
Definition: DownBeat.h:131
size_t m_bufsiz
Definition: DownBeat.h:126
size_t m_beatframesize
Definition: DownBeat.h:128
static int getHighestSupportedFactor()
Definition: Decimator.h:57
void pushAudioBlock(const float *audio)
For your downsampling convenience: call this function repeatedly with input audio blocks containing d...
Definition: DownBeat.cpp:94
DownBeat(float originalSampleRate, size_t decimationFactor, size_t dfIncrement)
Construct a downbeat locator that will operate on audio at the downsampled by the given decimation fa...
Definition: DownBeat.cpp:26
size_t m_factor
Definition: DownBeat.h:120
vector< double > d_vec_t
Definition: DownBeat.h:112
void forward(const double *realIn, double *realOut, double *imagOut)
Carry out a forward real-to-complex transform of size nsamples, where nsamples is the value provided ...
Definition: FFT.cpp:185
float * m_decbuf
Definition: DownBeat.h:125
float * m_buffer
Definition: DownBeat.h:124
#define EPS
size_t m_buffill
Definition: DownBeat.h:127
double * m_fftImagOut
Definition: DownBeat.h:132
void makeDecimators()
Definition: DownBeat.cpp:76
void getBeatSD(vector< double > &beatsd) const
Return the beat spectral difference function.
Definition: DownBeat.cpp:307
Decimator * m_decimator2
Definition: DownBeat.h:123
double * m_beatframe
Definition: DownBeat.h:129
void resetAudioBuffer()
Clear any buffered downsampled audio data.
Definition: DownBeat.cpp:137
#define TWO_PI
Definition: MathAliases.h:30
Definition: FFT.h:46
Decimator * m_decimator1
Definition: DownBeat.h:122
static int nextPowerOfTwo(int x)
Return the next higher integer power of two from x, e.g.
size_t m_increment
Definition: DownBeat.h:121
void findDownBeats(const float *audio, size_t audioLength, const vector< double > &beats, vector< int > &downbeats)
Estimate which beats are down-beats.
Definition: DownBeat.cpp:146
int m_bpb
Definition: DownBeat.h:118
Decimator carries out a fast downsample by a power-of-two factor.
Definition: Decimator.h:24
void setBeatsPerBar(int bpb)
Definition: DownBeat.cpp:70
d_vec_t m_beatsd
Definition: DownBeat.h:133
FFTReal * m_fft
Definition: DownBeat.h:130
double measureSpecDiff(d_vec_t oldspec, d_vec_t newspec)
Definition: DownBeat.cpp:259
static int getMax(double *data, unsigned int length, double *max=0)
float m_rate
Definition: DownBeat.h:119
~DownBeat()
Definition: DownBeat.cpp:57