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 SILENCE_CUTOFF 1e-10
406 #define DB_SILENCE_CUTOFF -100
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) < SILENCE_CUTOFF;
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 < SILENCE_CUTOFF ? DB_SILENCE_CUTOFF : (Real)10.0 * log10(value);
526 }
527 
528 inline Real lin2db(Real value, Real silenceCutoff, Real dbSilenceCutoff) {
529  return value < silenceCutoff ? dbSilenceCutoff : (Real)10.0 * log10(value);
530 }
531 
532 inline Real db2lin(Real value) {
533  return pow((Real)10.0, value/(Real)10.0);
534 }
535 
536 inline Real pow2db(Real power) {
537  return lin2db(power);
538 }
539 
540 inline Real pow2db(Real power, Real silenceCutoff, Real dbSilenceCutoff) {
541  return lin2db(power, silenceCutoff, dbSilenceCutoff);
542 }
543 
544 inline Real db2pow(Real power) {
545  return db2lin(power);
546 }
547 
548 inline Real amp2db(Real amplitude) {
549  return Real(2.0)*lin2db(amplitude);
550 }
551 
552 inline Real amp2db(Real amplitude, Real silenceCutoff, Real dbSilenceCutoff) {
553  return Real(2.0)*lin2db(amplitude, silenceCutoff, dbSilenceCutoff);
554 }
555 
556 inline Real db2amp(Real amplitude) {
557  return db2lin(0.5*amplitude);
558 }
559 
560 inline Real linear(Real input) {
561  return input;
562 }
563 
564 inline Real lin2log(Real input, Real silenceCutoff, Real logSilenceCutoff) {
565  return input < silenceCutoff ? logSilenceCutoff : log(input);
566 }
567 
568 
569 #ifdef OS_WIN32
570 // The following function hz2bark needs the function asinh,
571 // which is not included in ANSI math.h and thus does not
572 // "compile in windows". I copied the function asinh from boost
573 // http://www.boost.org/boost/math/special_functions/asinh.hpp
574 // Joachim, 25. June 2007
575 template<typename T>
576 inline T asinh(const T x)
577 {
578  using ::std::abs;
579  using ::std::sqrt;
580  using ::std::log;
581  using ::std::numeric_limits;
582 
583  T const one = static_cast<T>(1);
584  T const two = static_cast<T>(2);
585 
586  static T const taylor_2_bound = sqrt(numeric_limits<T>::epsilon());
587  static T const taylor_n_bound = sqrt(taylor_2_bound);
588  static T const upper_taylor_2_bound = one/taylor_2_bound;
589  static T const upper_taylor_n_bound = one/taylor_n_bound;
590 
591  if (x >= +taylor_n_bound) {
592  if (x > upper_taylor_n_bound) {
593  if (x > upper_taylor_2_bound) {
594  // approximation by laurent series in 1/x at 0+ order from -1 to 0
595  return( log( x * two) );
596  }
597  else {
598  // approximation by laurent series in 1/x at 0+ order from -1 to 1
599  return( log( x*two + (one/(x*two)) ) );
600  }
601  }
602  else {
603  return( log( x + sqrt(x*x+one) ) );
604  }
605  }
606  else if (x <= -taylor_n_bound) {
607  return(-asinh(-x));
608  }
609  else {
610  // approximation by taylor series in x at 0 up to order 2
611  T result = x;
612 
613  if (abs(x) >= taylor_2_bound)
614  {
615  T x3 = x*x*x;
616 
617  // approximation by taylor series in x at 0 up to order 4
618  result -= x3/static_cast<T>(6);
619  }
620  return(result);
621  }
622 }
623 #endif
624 
633 inline Real hz2bark(Real f) {
634  Real b = ((26.81*f)/(1960 + f)) - 0.53;
635 
636  if (b < 2) b += 0.15*(2-b);
637  if (b > 20.1) b += 0.22*(b - 20.1);
638 
639  return b;
640 }
641 
650 inline Real bark2hz(Real z) {
651  // note: these conditions have been deduced by inverting the ones from hz2bark
652  if (z < 2) z = (z - 0.3) / 0.85;
653  if (z > 20.1) z = (z - 4.422) / 1.22;
654 
655  // this part comes from Traunmüller's paper (could have deduced it also by myself... ;-) )
656  return 1960.0 * (z + 0.53) / (26.28 - z);
657 }
658 
660  return 52548.0 / (z*z - 52.56 * z + 690.39);
661 }
662 
663 
664 inline Real mel2hz(Real mel) {
665  return 700.0 * (exp(mel/1127.01048) - 1.0);
666 }
667 
668 inline Real mel102hz(Real mel) {
669  return 700.0 * (pow(10.0, mel/2595.0) - 1.0);
670 }
671 
672 // Convert Mel to Hz based on Slaney's formula in MATLAB Auditory Toolbox
673 inline Real mel2hzSlaney(Real mel) {
674  const Real minLogHz = 1000.0;
675  const Real linSlope = 3 / 200.;
676  const Real minLogMel = minLogHz * linSlope;
677 
678  if (mel < minLogMel) {
679  // Linear part: 0 - 1000 Hz.
680  return mel / linSlope;
681  }
682  else {
683  // Log-scale part: >= 1000 Hz.
684  const Real logStep = log(6.4) / 27.0;
685  return minLogHz * exp((mel - minLogMel) * logStep);
686  }
687 }
688 
689 inline Real hz2mel(Real hz) {
690  return 1127.01048 * log(hz/700.0 + 1.0);
691 }
692 
693 inline Real hz2mel10(Real hz) {
694  return 2595.0 * log10(hz/700.0 + 1.0);
695 }
696 
697 // Convert Hz to Mel based on Slaney's formula in MATLAB Auditory Toolbox
698 inline Real hz2melSlaney(Real hz) {
699  const Real minLogHz = 1000.0;
700  const Real linSlope = 3 / 200.;
701 
702  if (hz < minLogHz) {
703  // Linear part: 0 - 1000 Hz.
704  return hz * linSlope;
705  }
706  else {
707  // Log-scale part: >= 1000 Hz.
708  const Real minLogMel = minLogHz * linSlope;
709  const Real logStep = log(6.4) / 27.0;
710  return minLogMel + log(hz/minLogHz) / logStep;
711  }
712 }
713 
714 inline Real hz2hz(Real hz){
715  return hz;
716 }
717 
718 inline Real hz2cents(Real hz) {
719  return 12 * std::log(hz/440)/std::log(2.) + 69;
720 }
721 
722 inline int argmin(const std::vector<Real>& input) {
723  return std::min_element(input.begin(), input.end()) - input.begin();
724 }
725 
726 inline int argmax(const std::vector<Real>& input) {
727  return std::max_element(input.begin(), input.end()) - input.begin();
728 }
729 
730 // normalize a vector so its largest value gets mapped to 1
731 // if zero, the vector isn't touched
732 template <typename T> void normalize(std::vector<T>& array) {
733  if (array.empty()) return;
734 
735  T maxElement = *std::max_element(array.begin(), array.end());
736 
737  if (maxElement != (T) 0.0) {
738  for (uint i=0; i<array.size(); i++) {
739  array[i] /= maxElement;
740  }
741  }
742 }
743 
744 // normalize to the max(abs(array))
745 template <typename T> void normalizeAbs(std::vector<T>& array) {
746  if (array.empty()) return;
747  std::vector<T> absArray = array;
748  rectify(absArray);
749  T maxElement = *std::max_element(absArray.begin(), absArray.end());
750 
751  if (maxElement != (T) 0.0) {
752  for (uint i=0; i<array.size(); i++) {
753  array[i] /= maxElement;
754  }
755  }
756 }
757 
758 // normalize a vector so it's sum is equal to 1. the vector is not touched if
759 // it contains negative elements or the sum is zero
760 template <typename T> void normalizeSum(std::vector<T>& array) {
761  if (array.empty()) return;
762 
763  //T sumElements = std::accumulate(array.begin(), array.end(), (T) 0.0);
764  T sumElements = (T) 0.;
765  for (size_t i=0; i<array.size(); ++i) {
766  if (array[i] < 0) return;
767  sumElements += array[i];
768  }
769 
770  if (sumElements != (T) 0.0) {
771  for (size_t i=0; i<array.size(); ++i) {
772  array[i] /= sumElements;
773  }
774  }
775 }
776 
777 // returns the difference and approximate derivative vector of a vector
778 // derivative(x), for a vector x, is [x(1)-x(0) x(2)-x(1) ... x(n-1)-x(n-2)]
779 template <typename T>
780 std::vector<T> derivative(const std::vector<T>& array) {
781  if (array.size() < 2) {
782  throw EssentiaException("trying to calculate approximate derivative of empty or single-element array");
783  }
784 
785  std::vector<T> result(array.size()-1, (T)0.0);
786  for (int i=0; i<(int)result.size(); i++) {
787  result[i] = array[i+1] - array[i];
788  }
789  return result;
790 }
791 
792 template<typename T, typename U, typename Comparator=std::greater<T> >
793 class PairCompare : public std::binary_function<T, U, bool> {
794  Comparator _cmp;
795  public:
796  bool operator () (const std::pair<T,U>& p1, const std::pair<T,U>& p2) const {
797  if (_cmp(p1.first, p2.first)) return true;
798  if (_cmp(p2.first, p1.first)) return false;
799  return _cmp(p1.second, p2.second);
800  }
801 };
802 
803 // sorts two vectors by the cmp function. If the first elements of the pairs
804 // are equal, then it sorts by using cmp on the second value of the pair
805 template <typename T, typename U, typename Comparator>
806 void sortpair(std::vector<T>& v1, std::vector<U>& v2) {
807  if (v1.size() != v2.size()) {
808  throw EssentiaException("Cannot sort vectors of different size");
809  }
810  int size = v1.size();
811  std::vector<std::pair<T, U> > tmp(size);
812  for (int i=0; i<size; i++)
813  tmp[i] = std::make_pair(v1[i], v2[i]);
814  std::sort(tmp.begin(), tmp.end(), PairCompare<T, U, Comparator>());
815  for (int i=0; i<size; i++) {
816  v1[i] = tmp[i].first;
817  v2[i] = tmp[i].second;
818  }
819 }
820 
821 
822 // returns whether a number is a denormal number or not
823 inline bool isDenormal(const float& x) {
824  return std::fpclassify(x) == FP_SUBNORMAL;
825 
826  // old version: works only for i386 and float
827  //const int& xbits = reinterpret_cast<const int&>(x);
828  //const int absMantissa = xbits & 0x007FFFFF;
829  //const int biasedExponent = xbits & 0x7F800000;
830  //return (biasedExponent == 0 && absMantissa != 0);
831 }
832 
833 // should always return a positive value, even when a/b is negative
834 template <typename T> T fmod(T a, T b) {
835  T q = floor(a/b);
836  return a - q*b;
837 }
838 
839 // returns the principal phase argument between [-PI,PI]
840 template <typename T> T princarg(T y) {
841  T x = essentia::fmod(y + M_PI, M_2PI);
842  //if (x < 0) x += M_2PI; // should be useless with our implementation of fmod
843  return x - M_PI;
844 }
845 
857 template <typename T>
858 void hist(const T* array, uint n, int* n_array, T* x_array, uint n_bins) {
859  T miny = *std::min_element(array, array+n);
860  T maxy = *std::max_element(array, array+n);
861 
862  // x contains the center of the bins
863  for (uint i=0; i<n_bins; i++) {
864  x_array[i] = (0.5 + i)*(maxy - miny)/n_bins + miny;
865  }
866 
867  // cutoff contains the boundaries between the bins
868  std::vector<T> cutoff(n_bins - 1);
869  for (uint i=0; i<n_bins-1; i++) {
870  cutoff[i] = (x_array[i] + x_array[i+1]) / 2.0;
871  }
872 
873  // algo: either build the cumulated histogram by 2-level
874  // for-loops, and then build the diff, or first
875  // sort the distribution and do it directly.
876  // 1st version: O(n^2) time, O(1) space
877  // 2nd version: O(n·log(n)) time, O(n) space
878 
879  // implementation of 2nd version
880  std::vector<T> dist(array, array+n);
881  std::sort(dist.begin(), dist.end());
882  uint current_cutoff_idx = 0;
883  T current_cutoff = cutoff[0];
884  for (uint i=0; i<n_bins; i++) n_array[i] = 0;
885 
886  for (uint i=0; i<n; i++) {
887  while (dist[i] > current_cutoff) {
888  // last case; skip the rest and fill in the last bin
889  if (current_cutoff_idx == n_bins-2) {
890  n_array[n_bins-1] = n-i; // fill in the last bin with what's left
891  i = n; // to jump out of the 2nd loop (the 'for' one)
892  n_array[n_bins-2]--; // to compensate for the last one that will be added before jumping out of the loop
893  break;
894  }
895  current_cutoff_idx++;
896  current_cutoff = cutoff[current_cutoff_idx];
897  }
898  n_array[current_cutoff_idx]++;
899  }
900 }
901 
902 
906 template <typename T>
907 void bincount(const std::vector<T>& input, std::vector<T>& output) {
908  output.clear();
909  output.resize( (int) ( std::max<Real>( input[argmax(input)], 0.) + 0.5 ) + 1);
910  uint index = 0;
911  for (uint i=0; i< input.size(); i++) {
912  index = int(std::max<Real>(input[i],0) + 0.5);
913  if (index < output.size() ) {
914  output[index] += 1.;
915  }
916  }
917 }
918 
919 
924 template <typename T>
925 std::vector<std::vector<T> > transpose(const std::vector<std::vector<T> >& m) {
926  if (m.empty()) return std::vector<std::vector<T> >();
927 
928  int nrows = m.size();
929  int ncols = m[0].size();
930  for (int i=1; i<nrows; i++) {
931  if ((int)m[i].size() != ncols) {
932  std::ostringstream ss;
933  ss <<"Trying to transpose a non rectangular matrix. Expecting dim2 = " << ncols
934  << " but got " << m[i].size() << ". Cannot transpose!";
935  throw EssentiaException(ss.str());
936  }
937  }
938 
939  std::vector<std::vector<T> > result(ncols, std::vector<Real>(nrows));
940  for (int i=0; i<nrows; i++) {
941  for (int j=0; j<ncols; j++) {
942  result[j][i] = m[i][j];
943  }
944  }
945 
946  return result;
947 }
948 
949 template <typename T>
951  if (m.dim1() == 0) return TNT::Array2D<T>();
952 
953  int nrows = m.dim1();
954  int ncols = m.dim2();
955 
956  TNT::Array2D<T> result(ncols, nrows);
957  for (int i=0; i<nrows; i++) {
958  for (int j=0; j<ncols; j++) {
959  result[j][i] = m[i][j];
960  }
961  }
962 
963  return result;
964 }
965 
966 inline std::string equivalentKey(const std::string key) {
967  if (key == "C")
968  return "C";
969 
970  if (key == "C#")
971  return "Db";
972 
973  if (key == "Db")
974  return "C#";
975 
976  if (key == "D")
977  return "D";
978 
979  if (key == "D#")
980  return "Eb";
981 
982  if (key == "Eb")
983  return "D#";
984 
985  if (key == "E")
986  return "E";
987 
988  if (key == "F")
989  return "F";
990 
991  if (key == "F#")
992  return "Gb";
993 
994  if (key == "Gb")
995  return "F#";
996 
997  if (key == "G")
998  return "G";
999 
1000  if (key == "G#")
1001  return "Ab";
1002 
1003  if (key == "Ab")
1004  return "G#";
1005 
1006  if (key == "A")
1007  return "A";
1008 
1009  if (key == "A#")
1010  return "Bb";
1011 
1012  if (key == "Bb")
1013  return "A#";
1014 
1015  if (key == "B")
1016  return "B";
1017 
1018  return "";
1019 }
1020 
1021 
1025 template <typename T> T covariance(const std::vector<T>& x, const T xMean, const std::vector<T>& y, const T yMean) {
1026  if (x.empty())
1027  throw EssentiaException("trying to calculate covariance of empty array");
1028  if (y.empty())
1029  throw EssentiaException("trying to calculate covariance of empty array");
1030  if (x.size() != y.size())
1031  throw EssentiaException("x and y should have the same size");
1032 
1033  T cov = (T) 0.0;
1034 
1035  for (uint i=0; i<x.size(); i++) {
1036  cov += (x[i] - xMean) * (y[i] - yMean);
1037  }
1038 
1039  return (T)(cov / (Real)x.size());
1040 }
1041 
1046 template <typename T> T pearsonCorrelationCoefficient(const std::vector<T>& x, const std::vector<T>& y) {
1047  if (x.empty())
1048  throw EssentiaException("trying to calculate covariance of empty array");
1049  if (y.empty())
1050  throw EssentiaException("trying to calculate covariance of empty array");
1051  if (x.size() != y.size())
1052  throw EssentiaException("x and y should have the same size");
1053 
1054  T xMean = mean(x);
1055  T yMean = mean(y);
1056 
1057  T cov = covariance(x, xMean, y, yMean);
1058 
1059  T xStddev = stddev(x, xMean);
1060  T yStddev = stddev(y, yMean);
1061 
1062  // When dealing with constants corraltion is 0 by convention.
1063  if ((xStddev == (T)0.0) || (xStddev == (T)0.0) || (xStddev == (T)0.0)) return (T) 0.0;
1064 
1065  T corr = cov / (xStddev * yStddev);
1066 
1067  // Numerical error can yield results slightly outside the analytical range [-1, 1].
1068  // Clipping the output is a cheap way to mantain this contrain.
1069  // Seen in https://github.com/numpy/numpy/blob/v1.15.0/numpy/lib/function_base.py#L2403-L2406
1070  return std::max(std::min(corr, (T)1.0), (T)-1.0);
1071 }
1072 
1073 } // namespace essentia
1074 
1075 #endif // ESSENTIA_MATH_H
Real lin2log(Real input, Real silenceCutoff, Real logSilenceCutoff)
Definition: essentiamath.h:564
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:780
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:966
Real hz2mel(Real hz)
Definition: essentiamath.h:689
Real hz2cents(Real hz)
Definition: essentiamath.h:718
T median(const std::vector< T > &array)
Definition: essentiamath.h:354
Real hz2melSlaney(Real hz)
Definition: essentiamath.h:698
Real linear(Real input)
Definition: essentiamath.h:560
Real mel2hz(Real mel)
Definition: essentiamath.h:664
T stddev(const std::vector< T > &array, const T mean)
Definition: essentiamath.h:511
bool isDenormal(const float &x)
Definition: essentiamath.h:823
void bincount(const std::vector< T > &input, std::vector< T > &output)
Definition: essentiamath.h:907
#define M_2PI
Definition: essentiamath.h:41
#define SILENCE_CUTOFF
Definition: essentiamath.h:405
T round(const T value)
Definition: essentiamath.h:519
std::vector< T > varianceFrames(const std::vector< std::vector< T > > &frames)
Definition: essentiamath.h:264
Real mel2hzSlaney(Real mel)
Definition: essentiamath.h:673
TNT::Array2D< T > meanMatrix(const std::vector< TNT::Array2D< T > * > &array)
Definition: essentiamath.h:167
Real db2pow(Real power)
Definition: essentiamath.h:544
T instantPower(const std::vector< T > &array)
Definition: essentiamath.h:391
Real amp2db(Real amplitude)
Definition: essentiamath.h:548
int argmin(const std::vector< Real > &input)
Definition: essentiamath.h:722
Definition: tnt_array2d.h:37
T fmod(T a, T b)
Definition: essentiamath.h:834
Real hz2hz(Real hz)
Definition: essentiamath.h:714
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:840
Real bark2hz(Real z)
Definition: essentiamath.h:650
TNT::Array2D< T > & matinit(TNT::Array2D< T > &A)
Definition: tnt2essentiautils.h:45
void normalizeAbs(std::vector< T > &array)
Definition: essentiamath.h:745
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:793
#define DB_SILENCE_CUTOFF
Definition: essentiamath.h:406
float Real
Definition: types.h:68
Real db2lin(Real value)
Definition: essentiamath.h:532
std::vector< std::vector< T > > transpose(const std::vector< std::vector< T > > &m)
Definition: essentiamath.h:925
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:732
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:536
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:693
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:659
bool isSilent(const std::vector< T > &array)
Definition: essentiamath.h:409
unsigned int uint
Definition: types.h:48
T pearsonCorrelationCoefficient(const std::vector< T > &x, const std::vector< T > &y)
Definition: essentiamath.h:1046
T nextPowerTwo(T n)
Definition: essentiamath.h:64
Real db2amp(Real amplitude)
Definition: essentiamath.h:556
T mean(const std::vector< T > &array, int start, int end)
Definition: essentiamath.h:142
Real hz2bark(Real f)
Definition: essentiamath.h:633
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:794
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:796
void normalizeSum(std::vector< T > &array)
Definition: essentiamath.h:760
Real mel102hz(Real mel)
Definition: essentiamath.h:668
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:858
int argmax(const std::vector< Real > &input)
Definition: essentiamath.h:726
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:806
T covariance(const std::vector< T > &x, const T xMean, const std::vector< T > &y, const T yMean)
Definition: essentiamath.h:1025