qm-dsp
1.8
|
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