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
DownBeat(float originalSampleRate,size_t decimationFactor,size_t dfIncrement)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)
45 m_beatframesize = MathUtilities::nextPowerOfTwo
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];
54 m_fft = new FFTReal(m_beatframesize);
55 }
56
~DownBeat()57 DownBeat::~DownBeat()
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
setBeatsPerBar(int bpb)70 DownBeat::setBeatsPerBar(int bpb)
71 {
72 m_bpb = bpb;
73 }
74
75 void
makeDecimators()76 DownBeat::makeDecimators()
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) {
82 m_decimator1 = new Decimator(m_increment, m_factor);
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
pushAudioBlock(const float * audio)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);
113 m_decimator2->process(m_decbuf, m_buffer + m_buffill);
114 } else if (m_decimator1) {
115 m_decimator1->process(audio, m_buffer + m_buffill);
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;
126 m_buffill += m_increment / m_factor;
127 }
128
129 const float *
getBufferedAudio(size_t & length) const130 DownBeat::getBufferedAudio(size_t &length) const
131 {
132 length = m_buffill;
133 return m_buffer;
134 }
135
136 void
resetAudioBuffer()137 DownBeat::resetAudioBuffer()
138 {
139 if (m_buffer) free(m_buffer);
140 m_buffer = 0;
141 m_buffill = 0;
142 m_bufsiz = 0;
143 }
144
145 void
findDownBeats(const float * audio,size_t audioLength,const d_vec_t & beats,i_vec_t & downbeats)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
199 m_fft->forward(m_beatframe, m_fftRealOut, m_fftImagOut);
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
210 MathUtilities::adaptiveThreshold(newspec);
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
measureSpecDiff(d_vec_t oldspec,d_vec_t newspec)259 DownBeat::measureSpecDiff(d_vec_t oldspec, d_vec_t newspec)
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
getBeatSD(vector<double> & beatsd) const307 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
312