Essentia  2.1-beta6-dev
bpmutil.h
Go to the documentation of this file.
1 
2 /*
3  * Copyright (C) 2006-2021 Music Technology Group - Universitat Pompeu Fabra
4  *
5  * This file is part of Essentia
6  *
7  * Essentia is free software: you can redistribute it and/or modify it under
8  * the terms of the GNU Affero General Public License as published by the Free
9  * Software Foundation (FSF), either version 3 of the License, or (at your
10  * option) any later version.
11  *
12  * This program is distributed in the hope that it will be useful, but WITHOUT
13  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
15  * details.
16  *
17  * You should have received a copy of the Affero GNU General Public License
18  * version 3 along with this program. If not, see http://www.gnu.org/licenses/
19  */
20 
21 #ifndef ESSENTIA_BPMUTILS_H
22 #define ESSENTIA_BPMUTILS_H
23 
24 #include "../essentiamath.h"
25 #include <cassert>
26 
27 namespace essentia {
28 
29 inline
30 Real lagToBpm(Real lag, Real sampleRate, Real hopSize) {
31  return 60.0*sampleRate/lag/hopSize;
32 }
33 
34 inline
35 Real bpmToLag(Real bpm, Real sampleRate, Real hopSize) {
36  return lagToBpm(bpm, sampleRate, hopSize);
37 }
38 
39 inline
40 int longestChain(const std::vector<Real>& dticks, int startpos, Real period, Real tolerance) {
41  int pos = startpos;
42  Real ubound = period*(1+tolerance);
43  Real lbound = period*(1-tolerance);
44  while ((pos < (int)dticks.size()) &&
45  (lbound < dticks[pos] && dticks[pos] < ubound)) {
46  pos++;
47  }
48 
49  return pos - startpos;
50 }
51 
52 inline
53 void bpmDistance(Real x, Real y, Real& error, Real& ratio) {
54  ratio = x/y;
55  error = -1;
56  if (ratio < 1) {
57  ratio = round(1./ratio);
58  error=(x*ratio-y)/std::min(y,Real(x*ratio))*100;
59  }
60  else {
61  ratio = round(ratio);
62  error = (x-y*ratio)/std::min(x,Real(y*ratio))*100;
63  }
64 }
65 
66 inline
67 bool areEqual(Real a, Real b, Real tolerance) {
68  //return fabs(a-b) <= epsilon;
69  Real epsilon = std::max(tolerance, std::numeric_limits<Real>::epsilon());
70  Real error = 0;
71  Real ratio = 0;
72  bpmDistance(a, b, error, ratio);
73  return (std::fabs(error) <= epsilon && int(ratio)==1);
74 }
75 
76 inline
77 bool areHarmonics(Real x, Real y, Real epsilon, bool bPower2) {
78  // FIXME epsilon must be in %. a strict choice could be 3
79  epsilon = std::max(epsilon, std::numeric_limits<Real>::epsilon());
80  Real ratio = 0;
81  Real error = 0;
82  bpmDistance(x, y, error, ratio);
83  error = std::fabs(error);
84  if (error <= epsilon) {
85  if (bPower2) return isPowerTwo(int(fabs(ratio)));
86  return true;
87  }
88  return false;
89 }
90 
91 inline
93  // FIXME epsilon must be in %. a strict choice could be 3
94  epsilon = std::max(epsilon, std::numeric_limits<Real>::epsilon());
95  if (x<y) return greatestCommonDivisor(y, x, epsilon);
96  if (x==0) return 0;
97  Real error = std::numeric_limits<int>::max();
98  Real ratio = std::numeric_limits<int>::max();
99  bpmDistance(x, y, error, ratio);
100  if (fabs(error) <= epsilon) return y;
101  int a = int(x+0.5);
102  int b = int(y+0.5);
103  while (fabs(error) > epsilon) {
104  bpmDistance(a, b, error, ratio);
105  int remainder = a%b;
106  a = b;
107  b = remainder;
108  }
109  return a;
110 }
111 
112 
113 inline
114 std::vector<Real> roundBpms(const std::vector<Real>& bpms) {
115  Real epsilon = 3; // 3%
116  Real mainBpm=bpms[0];
117  std::vector<Real> harmonicBpms;
118  harmonicBpms.reserve(bpms.size());
119  for (int i=0; i<int(bpms.size()); i++) {
120  Real ratio=bpms[0]/mainBpm;
121  if (ratio < Real(1.0)) ratio = 1.0/ratio;
122  ratio = round(ratio*10.)/10.; // rounding to 1 float pos
123  int iRatio = int(ratio);
124  if (ratio-iRatio <= 0.100001) { // allow 2.9, 3.1 be considered as 3
125  harmonicBpms.push_back(bpms[i]);
126  }
127  else {
128  if ((ratio-iRatio) == 0.5) { // only interested in pure halfs
129  harmonicBpms.push_back(bpms[i]);
130  harmonicBpms.push_back(greatestCommonDivisor(bpms[i], mainBpm, epsilon));
131  }
132  }
133  }
134  return harmonicBpms;
135 }
136 
137 
138 // original postprocessticks from essentia 1.0
139 inline
140 std::vector<Real> postProcessTicks(const std::vector<Real>& origticks) {
141  if (origticks.size() < 3) return origticks;
142 
143  // find the most likely beat period
144  const int nticks = origticks.size();
145  std::vector<Real> dticks(nticks-1);
146 
147  for (int i=0; i<nticks-1; i++) dticks[i] = origticks[i+1] - origticks[i];
148 
149  // we might have had 6 secs frames during which we didn't find any beat, in which
150  // case we'll have one huge dtick value, which we actually want to prune
151  for (int i=0; i<(int)dticks.size(); i++) {
152  if (dticks[i] > 2.) {
153  dticks.erase(dticks.begin() + i);
154  i--;
155  }
156  }
157 
158  const int nbins = 100;
159  std::vector<int> dist(nbins);
160  std::vector<Real> distx(nbins);
161 
162  hist(&dticks[0], nticks-1, &dist[0], &distx[0], nbins);
163 
164  int maxidx = max_element(dist.begin(), dist.end()) - dist.begin();
165  Real maxbinCenter = distx[maxidx];
166 
167  // find the longest sequence of beats which has a fixed period of the previously
168  // found value; use a tolerance of about 10%
169  // Note: this favors high BPMs, because they will have more beats in the same amount of time
170  int maxl = 0;
171  int idx = 0;
172  Real period = maxbinCenter;
173 
174  for (int startpos = 0; startpos < nticks-1; startpos++) {
175  int l = longestChain(dticks, startpos, period, 0.1);
176  if (l > maxl) {
177  maxl = l;
178  idx = startpos;
179  }
180  }
181 
182  if (idx == 0 && maxl == 0) {
183  std::cerr << "Internal error while processing the beats, returning the original ones" << std::endl;
184  return origticks;
185  }
186 
187  // let's assume those beats are correct, and try to replace all the other ones
188  // with respect to the fixed period we have and the old positions of the beats
189  std::deque<Real> ticks(origticks.begin() + idx,
190  origticks.begin() + idx + maxl + 1);
191 
192  // take this value as the period for the whole track
193  Real targetPeriod = mean(dticks, idx, idx+maxl);
194  // 0.15, because 0.1 might be too strict, while 0.2 will allow false positives more easily
195  Real tolerance = 0.15 * targetPeriod;
196 
197 
198  // do the beats after the current beat base
199  Real cpos = ticks.back() + targetPeriod;
200  std::deque<Real> remaining(origticks.begin() + idx + maxl + 1,
201  origticks.end());
202 
203  while (!remaining.empty()) {
204  Real nbeat = remaining.front();
205 
206  if (nbeat < cpos - tolerance) {
207  // too far before, drop next beat
208  remaining.pop_front();
209  }
210  else {
211  // right in our expected zone, adjust the estimated beat to the one
212  // we actually found (NB: if we're too far away in front, just keep the
213  // beat as is)
214  if (nbeat < cpos + tolerance) {
215  cpos = nbeat;
216  remaining.pop_front();
217  }
218 
219  // in any case, mark the beat and jump on the next one
220  ticks.push_back(cpos);
221  cpos += targetPeriod;
222  }
223  }
224 
225  // do the beats before the current beat base
226  cpos = ticks.front() - targetPeriod;
227  remaining = std::deque<Real>(origticks.begin(),
228  origticks.begin() + idx);
229 
230  while (!remaining.empty()) {
231  Real nbeat = remaining.back();
232 
233  if (nbeat > cpos + tolerance) {
234  // too far after, drop beat
235  remaining.pop_back();
236  }
237  else {
238  // right in our expected zone, adjust the estimated beat to the one
239  // we actually found
240  if (nbeat > cpos - tolerance) {
241  cpos = nbeat;
242  remaining.pop_back();
243  }
244 
245  // in any case, mark the beat and jump on the next one
246  ticks.push_front(cpos);
247  cpos -= targetPeriod;
248  }
249  }
250 
251 
252  return std::vector<Real>(ticks.begin(), ticks.end());
253 }
254 
255  // modified version of the postprocessticks from tempotapticks, so it does not
256 // tend to favour fast bpms
257 inline
258 std::vector<Real> postProcessTicks(const std::vector<Real>& origticks,
259  const std::vector<Real>& ticksAmplitudes,
260  const Real& preferredPeriod) {
261  if (origticks.size() < 3) return origticks;
262 
263  // find the most likely beat period
264  const int nticks = origticks.size();
265  std::vector<Real> dticks(nticks-1);
266 
267  for (int i=0; i<nticks-1; i++) dticks[i] = origticks[i+1] - origticks[i];
268 
269  // we might have had 6 secs frames during which we didn't find any beat, in which
270  // case we'll have one huge dtick value, which we actually want to prune
271  for (int i=0; i<(int)dticks.size(); i++) {
272  if (dticks[i] > 2.) {
273  dticks.erase(dticks.begin() + i);
274  i--;
275  }
276  }
277 
278  const int nbins = 100;
279  std::vector<int> dist(nbins);
280  std::vector<Real> distx(nbins);
281 
282  hist(&dticks[0], dticks.size(), &dist[0], &distx[0], nbins);
283 
284  int maxidx = max_element(dist.begin(), dist.end()) - dist.begin();
285  Real maxbinCenter = distx[maxidx];
286 
287  // find the longest sequence of beats which has a fixed period of the previously
288  // found value; use a tolerance of about 10%
289  // Note: this favors high BPMs, because they will have more beats in the same amount of time
290  int maxl = 0;
291  int idx = 0;
292  Real period = maxbinCenter;
293 
294  //std::cout << "period: " << period << std::endl;
295  for (int startpos = 0; startpos < nticks-1; startpos++) {
296  int l = longestChain(dticks, startpos, period, 0.05);
297  if (l > maxl) {
298  maxl = l;
299  idx = startpos;
300  }
301  }
302 
303  Real targetPeriod = preferredPeriod;
304  if (idx ==0 && maxl==0) {
305  idx = argmax(ticksAmplitudes);
306  }
307  // take this value as the period for the whole track
308  else targetPeriod = mean(dticks, idx, idx+maxl);
309 
310  targetPeriod = (targetPeriod+preferredPeriod)/2.0;
311 
312  // if the targetPeriod is differs too much from the preferred period we
313  // search for the tick with max amplitude and take that as the reference tick
314  if (targetPeriod < 0.95*preferredPeriod || targetPeriod > 1.05*preferredPeriod) {
315  idx = idx + argmax(std::vector<Real>(ticksAmplitudes.begin()+idx, ticksAmplitudes.begin()+idx+maxl+1));
316  maxl = 0;
317  targetPeriod = preferredPeriod;
318  //std::cout << "Targets differ too much!. New target period will be the preferred one " << targetPeriod << std::endl;
319  }
320 
321  Real originalTargetPeriod = targetPeriod;
322 
323  // let's assume those beats are correct, and try to replace all the other ones
324  // with respect to the fixed period we have and the old positions of the beats
325  std::deque<Real> ticks(origticks.begin() + idx,
326  origticks.begin() + idx + maxl + 1);
327 
328  //fix tolerance at no more than 30ms:
329  Real tolerance =0.03;
330 
331  if (targetPeriod < 0.05) {
332  std::cout << "PostProcessTicks: target Period too short. Returning the original ticks" << std::endl;
333  return origticks;
334  }
335  Real cummulatedPeriod = targetPeriod;
336  int nAccumulations = 1;
337 
338 
339  // do the beats after the current beat base
340  Real cpos = ticks.back() + targetPeriod;
341  std::deque<Real> remaining(origticks.begin() + idx + maxl + 1,
342  origticks.end());
343 
344  while (!remaining.empty()) {
345  Real nbeat = remaining.front();
346 
347  if (nbeat < cpos - tolerance) {
348  // too far before, drop next beat
349  remaining.pop_front();
350  }
351  else {
352  // right in our expected zone, adjust the estimated beat to the one
353  // we actually found (NB: if we're too far away in front, just keep the
354  // beat as is)
355  if (nbeat < cpos + tolerance) {
356  cummulatedPeriod += (nbeat - (cpos - targetPeriod));
357  nAccumulations++;
358  targetPeriod = cummulatedPeriod/nAccumulations;
359  //std::cout << "new target Period: " << targetPeriod << " bpm: " << 60./targetPeriod << std::endl;
360  //std::cout << " \tbeat at: " << nbeat << " belongs to [" << cpos-tolerance << ", " << cpos+tolerance <<"], cpos: " << cpos << std::endl;
361  cpos = nbeat;
362  remaining.pop_front();
363  }
364 
365  // in any case, mark the beat and jump on the next one
366  ticks.push_back(cpos);
367  cpos += targetPeriod;
368  }
369  }
370 
371  // do the beats before the current beat base
372  cpos = ticks.front() - targetPeriod;
373  remaining = std::deque<Real>(origticks.begin(),
374  origticks.begin() + idx);
375 
376  targetPeriod = originalTargetPeriod;
377  cummulatedPeriod = targetPeriod;
378  nAccumulations = 1;
379 
380  while (!remaining.empty()) {
381  Real nbeat = remaining.back();
382 
383  if (nbeat > cpos + tolerance) {
384  // too far after, drop beat
385  remaining.pop_back();
386  }
387  else {
388  // right in our expected zone, adjust the estimated beat to the one
389  // we actually found
390  if (nbeat > cpos - tolerance) {
391  cummulatedPeriod += ((cpos + targetPeriod)-nbeat);
392  nAccumulations++;
393  targetPeriod = cummulatedPeriod/nAccumulations;
394  //std::cout << "new target Period: " << targetPeriod << " bpm: " << 60./targetPeriod << std::endl;
395  //std::cout << " \tbeat at: " << nbeat << " belongs to [" << cpos-tolerance << ", " << cpos+tolerance <<"], cpos: " << cpos << std::endl;
396  cpos = nbeat;
397  remaining.pop_back();
398  }
399 
400  // in any case, mark the beat and jump on the next one
401  ticks.push_front(cpos);
402  cpos -= targetPeriod;
403  }
404  }
405 
406  return std::vector<Real>(ticks.begin(), ticks.end());
407 }
408 } // namespace essentia
409 
410 #endif
Definition: algorithm.h:28
void hist(const T *array, uint n, int *n_array, T *x_array, uint n_bins)
Definition: essentiamath.h:899
T round(const T value)
Definition: essentiamath.h:538
T mean(const std::vector< T > &array, int start, int end)
Definition: essentiamath.h:142
void bpmDistance(Real x, Real y, Real &error, Real &ratio)
Definition: bpmutil.h:53
std::vector< Real > roundBpms(const std::vector< Real > &bpms)
Definition: bpmutil.h:114
std::vector< Real > postProcessTicks(const std::vector< Real > &origticks)
Definition: bpmutil.h:140
Real bpmToLag(Real bpm, Real sampleRate, Real hopSize)
Definition: bpmutil.h:35
Real lagToBpm(Real lag, Real sampleRate, Real hopSize)
Definition: bpmutil.h:30
Real greatestCommonDivisor(Real x, Real y, Real epsilon)
Definition: bpmutil.h:92
bool areEqual(Real a, Real b, Real tolerance)
Definition: bpmutil.h:67
float Real
Definition: types.h:69
bool isPowerTwo(T n)
Definition: essentiamath.h:45
bool areHarmonics(Real x, Real y, Real epsilon, bool bPower2)
Definition: bpmutil.h:77
int argmax(const std::vector< Real > &input)
Definition: essentiamath.h:751
int longestChain(const std::vector< Real > &dticks, int startpos, Real period, Real tolerance)
Definition: bpmutil.h:40