Essentia  2.1-beta5-dev
essentiamath.h
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2006-2016 Music Technology Group - Universitat Pompeu Fabra
3  *
4  * This file is part of Essentia
5  *
6  * Essentia is free software: you can redistribute it and/or modify it under
7  * the terms of the GNU Affero General Public License as published by the Free
8  * Software Foundation (FSF), either version 3 of the License, or (at your
9  * option) any later version.
10  *
11  * This program is distributed in the hope that it will be useful, but WITHOUT
12  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
13  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
14  * details.
15  *
16  * You should have received a copy of the Affero GNU General Public License
17  * version 3 along with this program. If not, see http://www.gnu.org/licenses/
18  */
19 
20 #ifndef ESSENTIA_MATH_H
21 #define ESSENTIA_MATH_H
22 
23 #ifndef _USE_MATH_DEFINES
24 #define _USE_MATH_DEFINES
25 #endif
26 
27 #include <math.h>
28 #include <cmath>
29 #include <vector>
30 #include <numeric>
31 #include <limits>
32 #include <functional>
33 #include <utility> // for pair
34 #include <sstream>
35 #include <algorithm> // for std::sort
36 #include <deque>
37 #include "types.h"
38 #include "utils/tnt/tnt.h"
40 
41 #define M_2PI (2 * M_PI)
42 
43 namespace essentia {
44 
45 template <typename T> bool isPowerTwo(T n) {
46  return (n & (n-1)) == 0;
47 }
48 
49 template <typename T> T log2(T x) {
50  return log(x) / M_LN2;
51 }
52 
53 template <typename T>
54 int ilog10(T n) {
55  if (n < 0) return ilog10(-n);
56  if (n < 10) return 0; // should return -infinite for n == 0
57  return 1 + ilog10(n/10);
58 }
59 
64 template <typename T> T nextPowerTwo(T n) {
65  n--;
66  n |= (n >> 1);
67  n |= (n >> 2);
68  n |= (n >> 4);
69  n |= (n >> 8);
70  n |= (n >> 16);
71  return ++n;
72 }
73 
74 template <> inline long long int nextPowerTwo(long long int n) {
75  n--;
76  n |= (n >> 1);
77  n |= (n >> 2);
78  n |= (n >> 4);
79  n |= (n >> 8);
80  n |= (n >> 16);
81  n |= (n >> 32);
82  return ++n;
83 }
84 
88 template <typename T> T norm(const std::vector<T>& array) {
89  if (array.empty()) {
90  throw EssentiaException("trying to calculate norm of empty array");
91  }
92 
93  T sum = (T) 0.0;
94 
95  for (uint i=0; i<array.size(); i++) {
96  sum += array[i] * array[i];
97  }
98 
99  return sqrt(sum);
100 }
101 
105 template <typename T> T sumSquare(const std::vector<T> array, const size_t start, const size_t end) {
106  T sum = 0.0;
107  for (size_t i = array; i < end; ++i) {
108  sum += array[i] * array[i];
109  }
110  return sum;
111 }
112 
116 template <typename T> T sum(const std::vector<T>& array, int start, int end) {
117  T sum = 0.0;
118  int i = start;
119 
120  for (; i<end-8; i+=8) {
121  sum += array[i];
122  sum += array[i+1];
123  sum += array[i+2];
124  sum += array[i+3];
125  sum += array[i+4];
126  sum += array[i+5];
127  sum += array[i+6];
128  sum += array[i+7];
129  }
130 
131  // do the rest of the loop
132  for (; i<end; i++) {
133  sum += array[i];
134  }
135 
136  return sum;
137 }
138 
142 template <typename T> T mean(const std::vector<T>& array, int start, int end) {
143  return sum(array, start, end) / (end - start);
144 }
145 
149 template <typename T> T sum(const std::vector<T>& array) {
150  if (array.empty()) return 0;
151  return sum(array, 0, array.size());
152 }
153 
157 template <typename T> T mean(const std::vector<T>& array) {
158  if (array.empty())
159  throw EssentiaException("trying to calculate mean of empty array");
160  return mean(array, 0, array.size());
161 }
162 
166 template <typename T>
167  TNT::Array2D<T> meanMatrix(const std::vector<TNT::Array2D<T>* >& array) {
168  if (array.empty())
169  throw EssentiaException("trying to calculate mean of empty array");
170  //return mean(array, 0, array.size());
171  TNT::Array2D<T> mean(array[0]->dim1(), array[0]->dim2());
172  matinit(mean);
173  for (int i = 0; i < (int)array.size(); i++) {
174  mean += *array[i];
175  }
176  mean /= (Real)array.size();
177  return mean;
178 }
179 
183 template <typename T>
184  TNT::Array2D<T> meanMatrix(const std::vector<TNT::Array2D<T> >& array) {
185  if (array.empty())
186  throw EssentiaException("trying to calculate mean of empty array");
187  //return mean(array, 0, array.size());
188  TNT::Array2D<T> mean(array[0].dim1(), array[0].dim2());
189  matinit(mean);
190  for (int i = 0; i < (int)array.size(); i++) {
191  mean += array[i];
192  }
193  mean /= (Real)array.size();
194  return mean;
195 }
196 
197 // returns the mean of frames
198 template <typename T>
199 std::vector<T> meanFrames(const std::vector<std::vector<T> >& frames, int beginIdx=0, int endIdx=-1) {
200  if (frames.empty()) {
201  throw EssentiaException("trying to calculate mean of empty array of frames");
202  }
203 
204  if (endIdx == -1) endIdx = (int)frames.size();
205  uint vsize = frames[0].size();
206 
207  std::vector<T> result(vsize, (T)0.0);
208  typename std::vector<std::vector<T> >::const_iterator it = frames.begin() + beginIdx;
209  typename std::vector<std::vector<T> >::const_iterator end = frames.begin() + endIdx;
210  for (; it!=end; ++it) {
211  typename std::vector<T>::const_iterator itFrame = it->begin();
212  typename std::vector<T>::const_iterator endFrame = it->end();
213  typename std::vector<T>::iterator itResult = result.begin();
214  for (; itFrame != endFrame; ++itFrame, ++itResult) {
215  *itResult += *itFrame;
216  }
217  }
218  for (uint j=0; j<vsize; j++) result[j] /= (endIdx - beginIdx);
219 
220  return result;
221 }
222 
223 // returns the median of frames
224 template <typename T>
225 std::vector<T> medianFrames(const std::vector<std::vector<T> >& frames, int beginIdx=0, int endIdx=-1) {
226  if (frames.empty()) {
227  throw EssentiaException("trying to calculate mean of empty array of frames");
228  }
229 
230  if (endIdx == -1) endIdx = (int)frames.size();
231 
232  uint vsize = frames[0].size();
233  uint fsize = endIdx - beginIdx;
234 
235  std::vector<T> result(vsize, (T)0.0);
236  std::vector<T> temp;
237  temp.reserve(fsize);
238 
239  for (uint i=0; i<vsize; ++i) {
240  typename std::vector<std::vector<T> >::const_iterator it = frames.begin() + beginIdx;
241  typename std::vector<std::vector<T> >::const_iterator end = frames.begin() + endIdx;
242 
243  temp.clear();
244  for (; it!=end; ++it) {
245  temp.push_back((*it)[i]);
246  }
247  std::sort(temp.begin(), temp.end());
248 
249  // array size is an odd number
250  if (fsize % 2 == 0.0) {
251  result[i] = (temp[uint(fsize/2 - 1)] + temp[uint(fsize/2)]) / 2;
252  }
253  // array size is an even number
254  else {
255  result[i] = temp[uint(fsize/2)];
256  }
257  }
258  return result;
259 }
260 
261 
262 // returns the variance of frames
263 template <typename T>
264 std::vector<T> varianceFrames(const std::vector<std::vector<T> >& frames) {
265  if (frames.empty()) {
266  throw EssentiaException("trying to calculate variance of empty array of frames");
267  }
268 
269  uint nframes = frames.size();
270  uint vsize = frames[0].size();
271 
272  std::vector<T> m = meanFrames(frames);
273 
274  std::vector<T> result(vsize, (T)0.0);
275  T diff;
276  for (uint i=0; i<nframes; i++) {
277  for (uint j=0; j<vsize; j++) {
278  diff = frames[i][j] - m[j];
279  result[j] += diff*diff;
280  }
281  }
282  for (uint j=0; j<vsize; j++) result[j] /= nframes;
283 
284  return result;
285 }
286 
287 
288 template <typename T>
289 std::vector<T> skewnessFrames(const std::vector<std::vector<T> >& frames) {
290  if (frames.empty()) {
291  throw EssentiaException("trying to calculate skewness of empty array of frames");
292  }
293 
294  uint nframes = frames.size();
295  uint vsize = frames[0].size();
296 
297  std::vector<T> m = meanFrames(frames);
298 
299  std::vector<T> result(vsize, (T)0.0);
300  std::vector<T> m3(vsize, (T)0.0);
301  std::vector<T> m2(vsize, (T)0.0);
302  T diff;
303  for (uint i=0; i<nframes; i++) {
304  for (uint j=0; j<vsize; j++) {
305  diff = frames[i][j] - m[j];
306  m2[j] += diff*diff;
307  m3[j] += diff*diff*diff;
308  }
309  }
310  for (uint j=0; j<vsize; j++) {
311  m2[j] /= nframes;
312  m3[j] /= nframes;
313  if (m2[j] == (T)0.) result[j] = (T)0.;
314  else result[j] = m3[j] / pow(m2[j], (T)1.5);
315  }
316 
317  return result;
318 }
319 
320 template <typename T>
321 std::vector<T> kurtosisFrames(const std::vector<std::vector<T> >& frames) {
322  if (frames.empty()) {
323  throw EssentiaException("trying to calculate kurtosis of empty array of frames");
324  }
325 
326  uint nframes = frames.size();
327  uint vsize = frames[0].size();
328 
329  std::vector<T> m = meanFrames(frames);
330 
331  std::vector<T> result(vsize, (T)0.0);
332  std::vector<T> m2(vsize, (T)0.0);
333  std::vector<T> m4(vsize, (T)0.0);
334  T diff;
335  for (uint i=0; i<nframes; i++) {
336  for (uint j=0; j<vsize; j++) {
337  diff = frames[i][j] - m[j];
338  m2[j] += diff*diff;
339  m4[j] += diff*diff*diff*diff;
340  }
341  }
342  for (uint j=0; j<vsize; j++) {
343  m2[j] /= nframes;
344  m4[j] /= nframes;
345  if (m2[j] == (T)0.) result[j] = (T)(-3.);
346  else result[j] = m4[j] / (m2[j]*m2[j]) - 3;
347  }
348 
349  return result;
350 }
351 
352 
353 // returns the median of an array
354 template <typename T> T median(const std::vector<T>& array) {
355  if (array.empty())
356  throw EssentiaException("trying to calculate median of empty array");
357 
358  // median has sense only on sorted array
359  std::vector<T> sorted_array = array;
360  std::sort(sorted_array.begin(), sorted_array.end());
361 
362  uint size = sorted_array.size();
363 
364  // array size is an odd number
365  if (size % 2 == 0.0) {
366  return (sorted_array[uint(size/2 - 1)] + sorted_array[uint(size/2)]) / 2;
367  }
368  // array size is an even number
369  else {
370  return sorted_array[uint(size/2)];
371  }
372 }
373 
374 // returns the absolute value of each element of the array
375 template <typename T>
376 void rectify(std::vector<T>& array) {
377  for (int i=0; i<(int)array.size(); i++) {
378  array[i] = fabs(array[i]);
379  }
380 }
381 
382 // returns the sum of the squared array = the energy of the array
383 template <typename T> T energy(const std::vector<T>& array) {
384  if (array.empty())
385  throw EssentiaException("trying to calculate energy of empty array");
386 
387  return inner_product(array.begin(), array.end(), array.begin(), (T)0.0);
388 }
389 
390 // returns the instantaneous power of an array
391 template <typename T> T instantPower(const std::vector<T>& array) {
392  return energy(array) / array.size();
393 }
394 
395 // silence_cutoff_dB = -60
396 // silence_cutoff = 10 ^ (silence_cutoff_dB / 10)
397 // silence cutoff has been set to -90 dB. The rationale behind it is that in
398 // principle we won't have absolute silence anywhere except for those few seconds
399 // left between one song and the next one on a CD, all the other frames will
400 // never have complete digital silence as any equipment will produce some kind
401 // of noise and even when music is rendered to file it is normally dithered
402 // first. For this reason we set the silence cutoff as what should be silence
403 // on a 16bit pcm file, that is (16bit - 1bit)*6.02 which yields -90.3 dB thus
404 // aproximately -90
405 #define silenceCutoff 1e-9
406 #define dbSilenceCutoff -90
407 
408 // returns true if the signal average energy is below a cutoff value, here -90dB
409 template <typename T> bool isSilent(const std::vector<T>& array) {
410  return instantPower(array) < silenceCutoff;
411 }
412 
413 // returns the variance of an array of TNT::Array2D<T> elements
414 template <typename T>
415  TNT::Array2D<T> varianceMatrix(const std::vector<TNT::Array2D<T> >& array, const TNT::Array2D<T> & mean) {
416  if (array.empty())
417  throw EssentiaException("trying to calculate variance of empty array");
418 
419  TNT::Array2D<T> variance(array[0].dim1(), array[0].dim2());
420  matinit(variance);
421 
422  for (int i=0; i<(int)array.size(); i++) {
423  TNT::Array2D<T> temp = array[i] - mean;
424  variance += temp * temp;
425  }
426 
427  return variance / (T)array.size();
428 }
429 // returns the variance of an array of TNT::Array2D<T>* elements
430 template <typename T>
431  TNT::Array2D<T> varianceMatrix(const std::vector<TNT::Array2D<T>* >& array, const TNT::Array2D<T> & mean) {
432  if (array.empty())
433  throw EssentiaException("trying to calculate variance of empty array");
434 
435  TNT::Array2D<T> variance(array[0]->dim1(), array[0]->dim2());
436  matinit(variance);
437 
438  for (int i=0; i<(int)array.size(); i++) {
439  TNT::Array2D<T> temp = *array[i] - mean;
440  variance += temp * temp;
441  }
442 
443  return variance / (T)array.size();
444 }
445 
446 // returns the variance of an array
447 template <typename T> T variance(const std::vector<T>& array, const T mean) {
448  if (array.empty())
449  throw EssentiaException("trying to calculate variance of empty array");
450 
451  T variance = (T) 0.0;
452 
453  for (uint i=0; i<array.size(); i++) {
454  T temp = array[i] - mean;
455  variance += temp * temp;
456  }
457 
458  return variance / array.size();
459 }
460 
461 // returns the skewness of an array
462 template <typename T> T skewness(const std::vector<T>& array, const T mean) {
463  if (array.empty())
464  throw EssentiaException("trying to calculate skewness of empty array");
465 
466  const int n = (int)array.size();
467  T m2 = (T)0.0, m3 = (T)0.0;
468 
469  for (int i=0; i<n; i++) {
470  T temp = array[i] - mean;
471  m2 += temp * temp;
472  m3 += temp * temp * temp;
473  }
474 
475  m2 /= n; m3 /= n;
476 
477  T result;
478  //if (std::isnan(result) || std::isinf(result)) return 0;
479  if (m2 == (T)0.) result = (T)0.;
480  else result = m3 / pow(m2, (T)1.5);
481 
482  return result;
483 }
484 
485 // returns the kurtosis of an array
486 template <typename T> T kurtosis(const std::vector<T>& array, const T mean) {
487  if (array.empty())
488  throw EssentiaException("trying to calculate kurtosis of empty array");
489 
490  const int n = (int)array.size();
491  T m2 = (T)0.0, m4 = (T)0.0;
492 
493  for (int i=0; i<n; i++) {
494  T temp = array[i] - mean;
495  m2 += temp * temp;
496  m4 += temp * temp * temp * temp;
497  }
498 
499  m2 /= n; m4 /= n;
500 
501  T result;
502  //if (std::isnan(result) || std::isinf(result)) return 0;
503  if (m2 == (T)0.) result = (T)(-3.);
504  else result = m4 / (m2*m2) - 3;
505 
506  return result;
507 }
508 
509 
510 // returns the standard deviation of an array
511 template <typename T> T stddev(const std::vector<T>& array, const T mean) {
512  if (array.empty())
513  throw EssentiaException("trying to calculate stddev of empty array");
514 
515  return (T)sqrt(variance(array, mean));
516 }
517 
518 // round a value to the nearest integer value
519 template <typename T> T round(const T value) {
520  return (T)std::floor(value + (T)0.5);
521 }
522 
523 
524 inline Real lin2db(Real value) {
525  return value < silenceCutoff ? dbSilenceCutoff : (Real)10.0 * log10(value);
526 }
527 
528 
529 inline Real db2lin(Real value) {
530  return pow((Real)10.0, value/(Real)10.0);
531 }
532 
533 inline Real pow2db(Real power) {
534  return lin2db(power);
535 }
536 
537 inline Real db2pow(Real power) {
538  return db2lin(power);
539 }
540 
541 inline Real amp2db(Real amplitude) {
542  return Real(2.0)*lin2db(amplitude);
543 }
544 
545 inline Real db2amp(Real amplitude) {
546  return db2lin(0.5*amplitude);
547 }
548 
549 inline Real linear(Real input) {
550  return input;
551 }
552 
553 #ifdef OS_WIN32
554 // The following function hz2bark needs the function asinh,
555 // which is not included in ANSI math.h and thus does not
556 // "compile in windows". I copied the function asinh from boost
557 // http://www.boost.org/boost/math/special_functions/asinh.hpp
558 // Joachim, 25. June 2007
559 template<typename T>
560 inline T asinh(const T x)
561 {
562  using ::std::abs;
563  using ::std::sqrt;
564  using ::std::log;
565  using ::std::numeric_limits;
566 
567  T const one = static_cast<T>(1);
568  T const two = static_cast<T>(2);
569 
570  static T const taylor_2_bound = sqrt(numeric_limits<T>::epsilon());
571  static T const taylor_n_bound = sqrt(taylor_2_bound);
572  static T const upper_taylor_2_bound = one/taylor_2_bound;
573  static T const upper_taylor_n_bound = one/taylor_n_bound;
574 
575  if (x >= +taylor_n_bound) {
576  if (x > upper_taylor_n_bound) {
577  if (x > upper_taylor_2_bound) {
578  // approximation by laurent series in 1/x at 0+ order from -1 to 0
579  return( log( x * two) );
580  }
581  else {
582  // approximation by laurent series in 1/x at 0+ order from -1 to 1
583  return( log( x*two + (one/(x*two)) ) );
584  }
585  }
586  else {
587  return( log( x + sqrt(x*x+one) ) );
588  }
589  }
590  else if (x <= -taylor_n_bound) {
591  return(-asinh(-x));
592  }
593  else {
594  // approximation by taylor series in x at 0 up to order 2
595  T result = x;
596 
597  if (abs(x) >= taylor_2_bound)
598  {
599  T x3 = x*x*x;
600 
601  // approximation by taylor series in x at 0 up to order 4
602  result -= x3/static_cast<T>(6);
603  }
604  return(result);
605  }
606 }
607 #endif
608 
617 inline Real hz2bark(Real f) {
618  Real b = ((26.81*f)/(1960 + f)) - 0.53;
619 
620  if (b < 2) b += 0.15*(2-b);
621  if (b > 20.1) b += 0.22*(b - 20.1);
622 
623  return b;
624 }
625 
634 inline Real bark2hz(Real z) {
635  // note: these conditions have been deduced by inverting the ones from hz2bark
636  if (z < 2) z = (z - 0.3) / 0.85;
637  if (z > 20.1) z = (z - 4.422) / 1.22;
638 
639  // this part comes from Traunmüller's paper (could have deduced it also by myself... ;-) )
640  return 1960.0 * (z + 0.53) / (26.28 - z);
641 }
642 
644  return 52548.0 / (z*z - 52.56 * z + 690.39);
645 }
646 
647 
648 inline Real mel2hz(Real mel) {
649  return 700.0 * (exp(mel/1127.01048) - 1.0);
650 }
651 
652 inline Real mel102hz(Real mel) {
653  return 700.0 * (pow(10.0, mel/2595.0) - 1.0);
654 }
655 
656 inline Real hz2mel(Real hz) {
657  return 1127.01048 * log(hz/700.0 + 1.0);
658 }
659 
660 inline Real hz2mel10(Real hz) {
661  return 2595.0 * log10(hz/700.0 + 1.0);
662 }
663 
664 inline Real hz2hz(Real hz){
665  return hz;
666 }
667 
668 inline Real hz2cents(Real hz) {
669  return 12 * std::log(hz/440)/std::log(2.) + 69;
670 }
671 
672 inline int argmin(const std::vector<Real>& input) {
673  return std::min_element(input.begin(), input.end()) - input.begin();
674 }
675 
676 inline int argmax(const std::vector<Real>& input) {
677  return std::max_element(input.begin(), input.end()) - input.begin();
678 }
679 
680 // normalize a vector so its largest value gets mapped to 1
681 // if zero, the vector isn't touched
682 template <typename T> void normalize(std::vector<T>& array) {
683  if (array.empty()) return;
684 
685  T maxElement = *std::max_element(array.begin(), array.end());
686 
687  if (maxElement != (T) 0.0) {
688  for (uint i=0; i<array.size(); i++) {
689  array[i] /= maxElement;
690  }
691  }
692 }
693 
694 // normalize a vector so it's sum is equal to 1. the vector is not touched if
695 // it contains negative elements or the sum is zero
696 template <typename T> void normalizeSum(std::vector<T>& array) {
697  if (array.empty()) return;
698 
699  //T sumElements = std::accumulate(array.begin(), array.end(), (T) 0.0);
700  T sumElements = (T) 0.;
701  for (size_t i=0; i<array.size(); ++i) {
702  if (array[i] < 0) return;
703  sumElements += array[i];
704  }
705 
706  if (sumElements != (T) 0.0) {
707  for (size_t i=0; i<array.size(); ++i) {
708  array[i] /= sumElements;
709  }
710  }
711 }
712 
713 // returns the difference and approximate derivative vector of a vector
714 // derivative(x), for a vector x, is [x(1)-x(0) x(2)-x(1) ... x(n-1)-x(n-2)]
715 template <typename T>
716 std::vector<T> derivative(const std::vector<T>& array) {
717  if (array.size() < 2) {
718  throw EssentiaException("trying to calculate approximate derivative of empty or single-element array");
719  }
720 
721  std::vector<T> result(array.size()-1, (T)0.0);
722  for (int i=0; i<(int)result.size(); i++) {
723  result[i] = array[i+1] - array[i];
724  }
725  return result;
726 }
727 
728 template<typename T, typename U, typename Comparator=std::greater<T> >
729 class PairCompare : public std::binary_function<T, U, bool> {
730  Comparator _cmp;
731  public:
732  bool operator () (const std::pair<T,U>& p1, const std::pair<T,U>& p2) const {
733  if (_cmp(p1.first, p2.first)) return true;
734  if (_cmp(p2.first, p1.first)) return false;
735  return _cmp(p1.second, p2.second);
736  }
737 };
738 
739 // sorts two vectors by the cmp function. If the first elements of the pairs
740 // are equal, then it sorts by using cmp on the second value of the pair
741 template <typename T, typename U, typename Comparator>
742 void sortpair(std::vector<T>& v1, std::vector<U>& v2) {
743  if (v1.size() != v2.size()) {
744  throw EssentiaException("Cannot sort vectors of different size");
745  }
746  int size = v1.size();
747  std::vector<std::pair<T, U> > tmp(size);
748  for (int i=0; i<size; i++)
749  tmp[i] = std::make_pair(v1[i], v2[i]);
750  std::sort(tmp.begin(), tmp.end(), PairCompare<T, U, Comparator>());
751  for (int i=0; i<size; i++) {
752  v1[i] = tmp[i].first;
753  v2[i] = tmp[i].second;
754  }
755 }
756 
757 
758 // returns whether a number is a denormal number or not
759 inline bool isDenormal(const float& x) {
760  return std::fpclassify(x) == FP_SUBNORMAL;
761 
762  // old version: works only for i386 and float
763  //const int& xbits = reinterpret_cast<const int&>(x);
764  //const int absMantissa = xbits & 0x007FFFFF;
765  //const int biasedExponent = xbits & 0x7F800000;
766  //return (biasedExponent == 0 && absMantissa != 0);
767 }
768 
769 // should always return a positive value, even when a/b is negative
770 template <typename T> T fmod(T a, T b) {
771  T q = floor(a/b);
772  return a - q*b;
773 }
774 
775 // returns the principal phase argument between [-PI,PI]
776 template <typename T> T princarg(T y) {
777  T x = essentia::fmod(y + M_PI, M_2PI);
778  //if (x < 0) x += M_2PI; // should be useless with our implementation of fmod
779  return x - M_PI;
780 }
781 
793 template <typename T>
794 void hist(const T* array, uint n, int* n_array, T* x_array, uint n_bins) {
795  T miny = *std::min_element(array, array+n);
796  T maxy = *std::max_element(array, array+n);
797 
798  // x contains the center of the bins
799  for (uint i=0; i<n_bins; i++) {
800  x_array[i] = (0.5 + i)*(maxy - miny)/n_bins + miny;
801  }
802 
803  // cutoff contains the boundaries between the bins
804  std::vector<T> cutoff(n_bins - 1);
805  for (uint i=0; i<n_bins-1; i++) {
806  cutoff[i] = (x_array[i] + x_array[i+1]) / 2.0;
807  }
808 
809  // algo: either build the cumulated histogram by 2-level
810  // for-loops, and then build the diff, or first
811  // sort the distribution and do it directly.
812  // 1st version: O(n^2) time, O(1) space
813  // 2nd version: O(n·log(n)) time, O(n) space
814 
815  // implementation of 2nd version
816  std::vector<T> dist(array, array+n);
817  std::sort(dist.begin(), dist.end());
818  uint current_cutoff_idx = 0;
819  T current_cutoff = cutoff[0];
820  for (uint i=0; i<n_bins; i++) n_array[i] = 0;
821 
822  for (uint i=0; i<n; i++) {
823  while (dist[i] > current_cutoff) {
824  // last case; skip the rest and fill in the last bin
825  if (current_cutoff_idx == n_bins-2) {
826  n_array[n_bins-1] = n-i; // fill in the last bin with what's left
827  i = n; // to jump out of the 2nd loop (the 'for' one)
828  n_array[n_bins-2]--; // to compensate for the last one that will be added before jumping out of the loop
829  break;
830  }
831  current_cutoff_idx++;
832  current_cutoff = cutoff[current_cutoff_idx];
833  }
834  n_array[current_cutoff_idx]++;
835  }
836 }
837 
838 
842 template <typename T>
843 void bincount(const std::vector<T>& input, std::vector<T>& output) {
844  output.clear();
845  output.resize( (int) ( std::max<Real>( input[argmax(input)], 0.) + 0.5 ) + 1);
846  uint index = 0;
847  for (uint i=0; i< input.size(); i++) {
848  index = int(std::max<Real>(input[i],0) + 0.5);
849  if (index < output.size() ) {
850  output[index] += 1.;
851  }
852  }
853 }
854 
855 
860 template <typename T>
861 std::vector<std::vector<T> > transpose(const std::vector<std::vector<T> >& m) {
862  if (m.empty()) return std::vector<std::vector<T> >();
863 
864  int nrows = m.size();
865  int ncols = m[0].size();
866  for (int i=1; i<nrows; i++) {
867  if ((int)m[i].size() != ncols) {
868  std::ostringstream ss;
869  ss <<"Trying to transpose a non rectangular matrix. Expecting dim2 = " << ncols
870  << " but got " << m[i].size() << ". Cannot transpose!";
871  throw EssentiaException(ss.str());
872  }
873  }
874 
875  std::vector<std::vector<T> > result(ncols, std::vector<Real>(nrows));
876  for (int i=0; i<nrows; i++) {
877  for (int j=0; j<ncols; j++) {
878  result[j][i] = m[i][j];
879  }
880  }
881 
882  return result;
883 }
884 
885 template <typename T>
887  if (m.dim1() == 0) return TNT::Array2D<T>();
888 
889  int nrows = m.dim1();
890  int ncols = m.dim2();
891 
892  TNT::Array2D<T> result(ncols, nrows);
893  for (int i=0; i<nrows; i++) {
894  for (int j=0; j<ncols; j++) {
895  result[j][i] = m[i][j];
896  }
897  }
898 
899  return result;
900 }
901 
902 inline std::string equivalentKey(const std::string key) {
903  if (key == "C")
904  return "C";
905 
906  if (key == "C#")
907  return "Db";
908 
909  if (key == "Db")
910  return "C#";
911 
912  if (key == "D")
913  return "D";
914 
915  if (key == "D#")
916  return "Eb";
917 
918  if (key == "Eb")
919  return "D#";
920 
921  if (key == "E")
922  return "E";
923 
924  if (key == "F")
925  return "F";
926 
927  if (key == "F#")
928  return "Gb";
929 
930  if (key == "Gb")
931  return "F#";
932 
933  if (key == "G")
934  return "G";
935 
936  if (key == "G#")
937  return "Ab";
938 
939  if (key == "Ab")
940  return "G#";
941 
942  if (key == "A")
943  return "A";
944 
945  if (key == "A#")
946  return "Bb";
947 
948  if (key == "Bb")
949  return "A#";
950 
951  if (key == "B")
952  return "B";
953 
954  return "";
955 }
956 
957 } // namespace essentia
958 
959 #endif // ESSENTIA_MATH_H
int dim1() const
Definition: tnt_array2d.h:231
T variance(const std::vector< T > &array, const T mean)
Definition: essentiamath.h:447
std::vector< T > derivative(const std::vector< T > &array)
Definition: essentiamath.h:716
T sum(const std::vector< T > &array, int start, int end)
Definition: essentiamath.h:116
void rectify(std::vector< T > &array)
Definition: essentiamath.h:376
Real lin2db(Real value)
Definition: essentiamath.h:524
std::string equivalentKey(const std::string key)
Definition: essentiamath.h:902
Real hz2mel(Real hz)
Definition: essentiamath.h:656
Real hz2cents(Real hz)
Definition: essentiamath.h:668
T median(const std::vector< T > &array)
Definition: essentiamath.h:354
Real linear(Real input)
Definition: essentiamath.h:549
Real mel2hz(Real mel)
Definition: essentiamath.h:648
T stddev(const std::vector< T > &array, const T mean)
Definition: essentiamath.h:511
bool isDenormal(const float &x)
Definition: essentiamath.h:759
void bincount(const std::vector< T > &input, std::vector< T > &output)
Definition: essentiamath.h:843
#define M_2PI
Definition: essentiamath.h:41
#define dbSilenceCutoff
Definition: essentiamath.h:406
T round(const T value)
Definition: essentiamath.h:519
std::vector< T > varianceFrames(const std::vector< std::vector< T > > &frames)
Definition: essentiamath.h:264
TNT::Array2D< T > meanMatrix(const std::vector< TNT::Array2D< T > * > &array)
Definition: essentiamath.h:167
Real db2pow(Real power)
Definition: essentiamath.h:537
T instantPower(const std::vector< T > &array)
Definition: essentiamath.h:391
Real amp2db(Real amplitude)
Definition: essentiamath.h:541
int argmin(const std::vector< Real > &input)
Definition: essentiamath.h:672
Definition: tnt_array2d.h:37
T fmod(T a, T b)
Definition: essentiamath.h:770
Real hz2hz(Real hz)
Definition: essentiamath.h:664
std::vector< T > meanFrames(const std::vector< std::vector< T > > &frames, int beginIdx=0, int endIdx=-1)
Definition: essentiamath.h:199
T princarg(T y)
Definition: essentiamath.h:776
Real bark2hz(Real z)
Definition: essentiamath.h:634
TNT::Array2D< T > & matinit(TNT::Array2D< T > &A)
Definition: tnt2essentiautils.h:45
T sumSquare(const std::vector< T > array, const size_t start, const size_t end)
Definition: essentiamath.h:105
int dim2() const
Definition: tnt_array2d.h:234
Definition: essentiamath.h:729
float Real
Definition: types.h:68
Real db2lin(Real value)
Definition: essentiamath.h:529
std::vector< std::vector< T > > transpose(const std::vector< std::vector< T > > &m)
Definition: essentiamath.h:861
TNT::Array2D< T > varianceMatrix(const std::vector< TNT::Array2D< T > > &array, const TNT::Array2D< T > &mean)
Definition: essentiamath.h:415
T norm(const std::vector< T > &array)
Definition: essentiamath.h:88
void normalize(std::vector< T > &array)
Definition: essentiamath.h:682
std::vector< T > kurtosisFrames(const std::vector< std::vector< T > > &frames)
Definition: essentiamath.h:321
bool isPowerTwo(T n)
Definition: essentiamath.h:45
Real pow2db(Real power)
Definition: essentiamath.h:533
std::vector< T > medianFrames(const std::vector< std::vector< T > > &frames, int beginIdx=0, int endIdx=-1)
Definition: essentiamath.h:225
Real hz2mel10(Real hz)
Definition: essentiamath.h:660
T energy(const std::vector< T > &array)
Definition: essentiamath.h:383
Definition: algorithm.h:28
Definition: types.h:76
Real barkCriticalBandwidth(Real z)
Definition: essentiamath.h:643
bool isSilent(const std::vector< T > &array)
Definition: essentiamath.h:409
unsigned int uint
Definition: types.h:48
#define silenceCutoff
Definition: essentiamath.h:405
T nextPowerTwo(T n)
Definition: essentiamath.h:64
Real db2amp(Real amplitude)
Definition: essentiamath.h:545
T mean(const std::vector< T > &array, int start, int end)
Definition: essentiamath.h:142
Real hz2bark(Real f)
Definition: essentiamath.h:617
T kurtosis(const std::vector< T > &array, const T mean)
Definition: essentiamath.h:486
T skewness(const std::vector< T > &array, const T mean)
Definition: essentiamath.h:462
Comparator _cmp
Definition: essentiamath.h:730
int ilog10(T n)
Definition: essentiamath.h:54
bool operator()(const std::pair< T, U > &p1, const std::pair< T, U > &p2) const
Definition: essentiamath.h:732
void normalizeSum(std::vector< T > &array)
Definition: essentiamath.h:696
Real mel102hz(Real mel)
Definition: essentiamath.h:652
T log2(T x)
Definition: essentiamath.h:49
void hist(const T *array, uint n, int *n_array, T *x_array, uint n_bins)
Definition: essentiamath.h:794
int argmax(const std::vector< Real > &input)
Definition: essentiamath.h:676
std::vector< T > skewnessFrames(const std::vector< std::vector< T > > &frames)
Definition: essentiamath.h:289
void sortpair(std::vector< T > &v1, std::vector< U > &v2)
Definition: essentiamath.h:742