Essentia  2.1-beta6-dev
essentiamath.h
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2006-2020 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) {
106  T sum = 0.0;
107  for (size_t i = 0; i < array.size(); ++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 // returns the sum of frames
289 template <typename T>
290 std::vector<T> sumFrames(const std::vector<std::vector<T> >& frames) {
291  if (frames.empty()) {
292  throw EssentiaException("sumFrames: trying to calculate sum of empty input frames");
293  }
294  size_t nframes = frames.size();
295  size_t vsize = frames[0].size();
296  std::vector<T> result(vsize, (T)0);
297  for (size_t j=0; j<vsize; j++) {
298  for (size_t i=0; i<nframes; i++) {
299  result[j] += frames[i][j];
300  }
301  }
302  return result;
303 }
304 
305 
306 template <typename T>
307 std::vector<T> skewnessFrames(const std::vector<std::vector<T> >& frames) {
308  if (frames.empty()) {
309  throw EssentiaException("trying to calculate skewness of empty array of frames");
310  }
311 
312  uint nframes = frames.size();
313  uint vsize = frames[0].size();
314 
315  std::vector<T> m = meanFrames(frames);
316 
317  std::vector<T> result(vsize, (T)0.0);
318  std::vector<T> m3(vsize, (T)0.0);
319  std::vector<T> m2(vsize, (T)0.0);
320  T diff;
321  for (uint i=0; i<nframes; i++) {
322  for (uint j=0; j<vsize; j++) {
323  diff = frames[i][j] - m[j];
324  m2[j] += diff*diff;
325  m3[j] += diff*diff*diff;
326  }
327  }
328  for (uint j=0; j<vsize; j++) {
329  m2[j] /= nframes;
330  m3[j] /= nframes;
331  if (m2[j] == (T)0.) result[j] = (T)0.;
332  else result[j] = m3[j] / pow(m2[j], (T)1.5);
333  }
334 
335  return result;
336 }
337 
338 template <typename T>
339 std::vector<T> kurtosisFrames(const std::vector<std::vector<T> >& frames) {
340  if (frames.empty()) {
341  throw EssentiaException("trying to calculate kurtosis of empty array of frames");
342  }
343 
344  uint nframes = frames.size();
345  uint vsize = frames[0].size();
346 
347  std::vector<T> m = meanFrames(frames);
348 
349  std::vector<T> result(vsize, (T)0.0);
350  std::vector<T> m2(vsize, (T)0.0);
351  std::vector<T> m4(vsize, (T)0.0);
352  T diff;
353  for (uint i=0; i<nframes; i++) {
354  for (uint j=0; j<vsize; j++) {
355  diff = frames[i][j] - m[j];
356  m2[j] += diff*diff;
357  m4[j] += diff*diff*diff*diff;
358  }
359  }
360  for (uint j=0; j<vsize; j++) {
361  m2[j] /= nframes;
362  m4[j] /= nframes;
363  if (m2[j] == (T)0.) result[j] = (T)(-3.);
364  else result[j] = m4[j] / (m2[j]*m2[j]) - 3;
365  }
366 
367  return result;
368 }
369 
370 
371 // returns the median of an array
372 template <typename T> T median(const std::vector<T>& array) {
373  if (array.empty())
374  throw EssentiaException("trying to calculate median of empty array");
375 
376  // median has sense only on sorted array
377  std::vector<T> sorted_array = array;
378  std::sort(sorted_array.begin(), sorted_array.end());
379 
380  uint size = sorted_array.size();
381 
382  // array size is an odd number
383  if (size % 2 == 0.0) {
384  return (sorted_array[uint(size/2 - 1)] + sorted_array[uint(size/2)]) / 2;
385  }
386  // array size is an even number
387  else {
388  return sorted_array[uint(size/2)];
389  }
390 }
391 
392 // returns the absolute value of each element of the array
393 template <typename T>
394 void rectify(std::vector<T>& array) {
395  for (int i=0; i<(int)array.size(); i++) {
396  array[i] = fabs(array[i]);
397  }
398 }
399 
400 // returns the sum of the squared array = the energy of the array
401 template <typename T> T energy(const std::vector<T>& array) {
402  if (array.empty())
403  throw EssentiaException("trying to calculate energy of empty array");
404 
405  return inner_product(array.begin(), array.end(), array.begin(), (T)0.0);
406 }
407 
408 // returns the instantaneous power of an array
409 template <typename T> T instantPower(const std::vector<T>& array) {
410  return energy(array) / array.size();
411 }
412 
413 // silence_cutoff_dB = -60
414 // silence_cutoff = 10 ^ (silence_cutoff_dB / 10)
415 // silence cutoff has been set to -90 dB. The rationale behind it is that in
416 // principle we won't have absolute silence anywhere except for those few seconds
417 // left between one song and the next one on a CD, all the other frames will
418 // never have complete digital silence as any equipment will produce some kind
419 // of noise and even when music is rendered to file it is normally dithered
420 // first. For this reason we set the silence cutoff as what should be silence
421 // on a 16bit pcm file, that is (16bit - 1bit)*6.02 which yields -90.3 dB thus
422 // aproximately -90
423 #define SILENCE_CUTOFF 1e-10
424 #define DB_SILENCE_CUTOFF -100
425 #define LOG_SILENCE_CUTOFF -23.025850929940457
426 
427 // returns true if the signal average energy is below a cutoff value, here -90dB
428 template <typename T> bool isSilent(const std::vector<T>& array) {
429  return instantPower(array) < SILENCE_CUTOFF;
430 }
431 
432 // returns the variance of an array of TNT::Array2D<T> elements
433 template <typename T>
434  TNT::Array2D<T> varianceMatrix(const std::vector<TNT::Array2D<T> >& array, const TNT::Array2D<T> & mean) {
435  if (array.empty())
436  throw EssentiaException("trying to calculate variance of empty array");
437 
438  TNT::Array2D<T> variance(array[0].dim1(), array[0].dim2());
439  matinit(variance);
440 
441  for (int i=0; i<(int)array.size(); i++) {
442  TNT::Array2D<T> temp = array[i] - mean;
443  variance += temp * temp;
444  }
445 
446  return variance / (T)array.size();
447 }
448 // returns the variance of an array of TNT::Array2D<T>* elements
449 template <typename T>
450  TNT::Array2D<T> varianceMatrix(const std::vector<TNT::Array2D<T>* >& array, const TNT::Array2D<T> & mean) {
451  if (array.empty())
452  throw EssentiaException("trying to calculate variance of empty array");
453 
454  TNT::Array2D<T> variance(array[0]->dim1(), array[0]->dim2());
455  matinit(variance);
456 
457  for (int i=0; i<(int)array.size(); i++) {
458  TNT::Array2D<T> temp = *array[i] - mean;
459  variance += temp * temp;
460  }
461 
462  return variance / (T)array.size();
463 }
464 
465 // returns the variance of an array
466 template <typename T> T variance(const std::vector<T>& array, const T mean) {
467  if (array.empty())
468  throw EssentiaException("trying to calculate variance of empty array");
469 
470  T variance = (T) 0.0;
471 
472  for (uint i=0; i<array.size(); i++) {
473  T temp = array[i] - mean;
474  variance += temp * temp;
475  }
476 
477  return variance / array.size();
478 }
479 
480 // returns the skewness of an array
481 template <typename T> T skewness(const std::vector<T>& array, const T mean) {
482  if (array.empty())
483  throw EssentiaException("trying to calculate skewness of empty array");
484 
485  const int n = (int)array.size();
486  T m2 = (T)0.0, m3 = (T)0.0;
487 
488  for (int i=0; i<n; i++) {
489  T temp = array[i] - mean;
490  m2 += temp * temp;
491  m3 += temp * temp * temp;
492  }
493 
494  m2 /= n; m3 /= n;
495 
496  T result;
497  //if (std::isnan(result) || std::isinf(result)) return 0;
498  if (m2 == (T)0.) result = (T)0.;
499  else result = m3 / pow(m2, (T)1.5);
500 
501  return result;
502 }
503 
504 // returns the kurtosis of an array
505 template <typename T> T kurtosis(const std::vector<T>& array, const T mean) {
506  if (array.empty())
507  throw EssentiaException("trying to calculate kurtosis of empty array");
508 
509  const int n = (int)array.size();
510  T m2 = (T)0.0, m4 = (T)0.0;
511 
512  for (int i=0; i<n; i++) {
513  T temp = array[i] - mean;
514  m2 += temp * temp;
515  m4 += temp * temp * temp * temp;
516  }
517 
518  m2 /= n; m4 /= n;
519 
520  T result;
521  //if (std::isnan(result) || std::isinf(result)) return 0;
522  if (m2 == (T)0.) result = (T)(-3.);
523  else result = m4 / (m2*m2) - 3;
524 
525  return result;
526 }
527 
528 
529 // returns the standard deviation of an array
530 template <typename T> T stddev(const std::vector<T>& array, const T mean) {
531  if (array.empty())
532  throw EssentiaException("trying to calculate stddev of empty array");
533 
534  return (T)sqrt(variance(array, mean));
535 }
536 
537 // round a value to the nearest integer value
538 template <typename T> T round(const T value) {
539  return (T)std::floor(value + (T)0.5);
540 }
541 
542 
543 inline Real lin2db(Real value) {
544  return value < SILENCE_CUTOFF ? DB_SILENCE_CUTOFF : (Real)10.0 * log10(value);
545 }
546 
547 inline Real lin2db(Real value, Real silenceCutoff, Real dbSilenceCutoff) {
548  return value < silenceCutoff ? dbSilenceCutoff : (Real)10.0 * log10(value);
549 }
550 
551 inline Real db2lin(Real value) {
552  return pow((Real)10.0, value/(Real)10.0);
553 }
554 
555 inline Real pow2db(Real power) {
556  return lin2db(power);
557 }
558 
559 inline Real pow2db(Real power, Real silenceCutoff, Real dbSilenceCutoff) {
560  return lin2db(power, silenceCutoff, dbSilenceCutoff);
561 }
562 
563 inline Real db2pow(Real power) {
564  return db2lin(power);
565 }
566 
567 inline Real amp2db(Real amplitude) {
568  return Real(2.0)*lin2db(amplitude);
569 }
570 
571 inline Real amp2db(Real amplitude, Real silenceCutoff, Real dbSilenceCutoff) {
572  return Real(2.0)*lin2db(amplitude, silenceCutoff, dbSilenceCutoff);
573 }
574 
575 inline Real db2amp(Real amplitude) {
576  return db2lin(0.5*amplitude);
577 }
578 
579 inline Real linear(Real input) {
580  return input;
581 }
582 
583 inline Real lin2log(Real value) {
584  return value < SILENCE_CUTOFF ? LOG_SILENCE_CUTOFF : log(value);
585 }
586 
587 inline Real lin2log(Real input, Real silenceCutoff, Real logSilenceCutoff) {
588  return input < silenceCutoff ? logSilenceCutoff : log(input);
589 }
590 
591 
592 #ifdef OS_WIN32
593 // The following function hz2bark needs the function asinh,
594 // which is not included in ANSI math.h and thus does not
595 // "compile in windows". I copied the function asinh from boost
596 // http://www.boost.org/boost/math/special_functions/asinh.hpp
597 // Joachim, 25. June 2007
598 template<typename T>
599 inline T asinh(const T x)
600 {
601  using ::std::abs;
602  using ::std::sqrt;
603  using ::std::log;
604  using ::std::numeric_limits;
605 
606  T const one = static_cast<T>(1);
607  T const two = static_cast<T>(2);
608 
609  static T const taylor_2_bound = sqrt(numeric_limits<T>::epsilon());
610  static T const taylor_n_bound = sqrt(taylor_2_bound);
611  static T const upper_taylor_2_bound = one/taylor_2_bound;
612  static T const upper_taylor_n_bound = one/taylor_n_bound;
613 
614  if (x >= +taylor_n_bound) {
615  if (x > upper_taylor_n_bound) {
616  if (x > upper_taylor_2_bound) {
617  // approximation by laurent series in 1/x at 0+ order from -1 to 0
618  return( log( x * two) );
619  }
620  else {
621  // approximation by laurent series in 1/x at 0+ order from -1 to 1
622  return( log( x*two + (one/(x*two)) ) );
623  }
624  }
625  else {
626  return( log( x + sqrt(x*x+one) ) );
627  }
628  }
629  else if (x <= -taylor_n_bound) {
630  return(-asinh(-x));
631  }
632  else {
633  // approximation by taylor series in x at 0 up to order 2
634  T result = x;
635 
636  if (abs(x) >= taylor_2_bound)
637  {
638  T x3 = x*x*x;
639 
640  // approximation by taylor series in x at 0 up to order 4
641  result -= x3/static_cast<T>(6);
642  }
643  return(result);
644  }
645 }
646 #endif
647 
656 inline Real hz2bark(Real f) {
657  Real b = ((26.81*f)/(1960 + f)) - 0.53;
658 
659  if (b < 2) b += 0.15*(2-b);
660  if (b > 20.1) b += 0.22*(b - 20.1);
661 
662  return b;
663 }
664 
673 inline Real bark2hz(Real z) {
674  // note: these conditions have been deduced by inverting the ones from hz2bark
675  if (z < 2) z = (z - 0.3) / 0.85;
676  if (z > 20.1) z = (z - 4.422) / 1.22;
677 
678  // this part comes from Traunmüller's paper (could have deduced it also by myself... ;-) )
679  return 1960.0 * (z + 0.53) / (26.28 - z);
680 }
681 
683  return 52548.0 / (z*z - 52.56 * z + 690.39);
684 }
685 
686 
687 inline Real mel2hz(Real mel) {
688  return 700.0 * (exp(mel/1127.01048) - 1.0);
689 }
690 
691 inline Real mel102hz(Real mel) {
692  return 700.0 * (pow(10.0, mel/2595.0) - 1.0);
693 }
694 
695 // Convert Mel to Hz based on Slaney's formula in MATLAB Auditory Toolbox
696 inline Real mel2hzSlaney(Real mel) {
697  const Real minLogHz = 1000.0;
698  const Real linSlope = 3 / 200.;
699  const Real minLogMel = minLogHz * linSlope;
700 
701  if (mel < minLogMel) {
702  // Linear part: 0 - 1000 Hz.
703  return mel / linSlope;
704  }
705  else {
706  // Log-scale part: >= 1000 Hz.
707  const Real logStep = log(6.4) / 27.0;
708  return minLogHz * exp((mel - minLogMel) * logStep);
709  }
710 }
711 
712 inline Real hz2mel(Real hz) {
713  return 1127.01048 * log(hz/700.0 + 1.0);
714 }
715 
716 inline Real hz2mel10(Real hz) {
717  return 2595.0 * log10(hz/700.0 + 1.0);
718 }
719 
720 // Convert Hz to Mel based on Slaney's formula in MATLAB Auditory Toolbox
721 inline Real hz2melSlaney(Real hz) {
722  const Real minLogHz = 1000.0;
723  const Real linSlope = 3 / 200.;
724 
725  if (hz < minLogHz) {
726  // Linear part: 0 - 1000 Hz.
727  return hz * linSlope;
728  }
729  else {
730  // Log-scale part: >= 1000 Hz.
731  const Real minLogMel = minLogHz * linSlope;
732  const Real logStep = log(6.4) / 27.0;
733  return minLogMel + log(hz/minLogHz) / logStep;
734  }
735 }
736 
737 inline Real hz2hz(Real hz){
738  return hz;
739 }
740 
741 inline Real hz2cents(Real hz) {
742  return 12 * std::log(hz/440)/std::log(2.) + 69;
743 }
744 
745 inline int argmin(const std::vector<Real>& input) {
746  if (input.empty())
747  throw EssentiaException("trying to get argmin of empty array");
748  return std::min_element(input.begin(), input.end()) - input.begin();
749 }
750 
751 inline int argmax(const std::vector<Real>& input) {
752  if (input.empty())
753  throw EssentiaException("trying to get argmax of empty array");
754  return std::max_element(input.begin(), input.end()) - input.begin();
755 }
756 
757 // normalize a vector so its largest value gets mapped to 1
758 // if zero, the vector isn't touched
759 template <typename T> void normalize(std::vector<T>& array) {
760  if (array.empty()) return;
761 
762  T maxElement = *std::max_element(array.begin(), array.end());
763 
764  if (maxElement != (T) 0.0) {
765  for (uint i=0; i<array.size(); i++) {
766  array[i] /= maxElement;
767  }
768  }
769 }
770 
771 // normalize to the max(abs(array))
772 template <typename T> void normalizeAbs(std::vector<T>& array) {
773  if (array.empty()) return;
774  std::vector<T> absArray = array;
775  rectify(absArray);
776  T maxElement = *std::max_element(absArray.begin(), absArray.end());
777 
778  if (maxElement != (T) 0.0) {
779  for (uint i=0; i<array.size(); i++) {
780  array[i] /= maxElement;
781  }
782  }
783 }
784 
785 // normalize a vector so it's sum is equal to 1. the vector is not touched if
786 // it contains negative elements or the sum is zero
787 template <typename T> void normalizeSum(std::vector<T>& array) {
788  if (array.empty()) return;
789 
790  //T sumElements = std::accumulate(array.begin(), array.end(), (T) 0.0);
791  T sumElements = (T) 0.;
792  for (size_t i=0; i<array.size(); ++i) {
793  if (array[i] < 0) return;
794  sumElements += array[i];
795  }
796 
797  if (sumElements != (T) 0.0) {
798  for (size_t i=0; i<array.size(); ++i) {
799  array[i] /= sumElements;
800  }
801  }
802 }
803 
804 // returns the difference and approximate derivative vector of a vector
805 // derivative(x), for a vector x, is [x(1)-x(0) x(2)-x(1) ... x(n-1)-x(n-2)]
806 template <typename T>
807 std::vector<T> derivative(const std::vector<T>& array) {
808  if (array.size() < 2) {
809  throw EssentiaException("trying to calculate approximate derivative of empty or single-element array");
810  }
811 
812  std::vector<T> result(array.size()-1, (T)0.0);
813  for (int i=0; i<(int)result.size(); i++) {
814  result[i] = array[i+1] - array[i];
815  }
816  return result;
817 }
818 
819 template<typename T, typename U, typename Comparator=std::greater<T> >
820 class PairCompare : public std::binary_function<T, U, bool> {
821  Comparator _cmp;
822  public:
823  bool operator () (const std::pair<T,U>& p1, const std::pair<T,U>& p2) const {
824  if (_cmp(p1.first, p2.first)) return true;
825  if (_cmp(p2.first, p1.first)) return false;
826  return _cmp(p1.second, p2.second);
827  }
828 };
829 
830 // sorts two vectors by the cmp function. If the first elements of the pairs
831 // are equal, then it sorts by using cmp on the second value of the pair
832 template <typename T, typename U, typename Comparator>
833 void sortpair(std::vector<T>& v1, std::vector<U>& v2) {
834  if (v1.size() != v2.size()) {
835  throw EssentiaException("Cannot sort vectors of different size");
836  }
837  int size = v1.size();
838  std::vector<std::pair<T, U> > tmp(size);
839  for (int i=0; i<size; i++)
840  tmp[i] = std::make_pair(v1[i], v2[i]);
841  std::sort(tmp.begin(), tmp.end(), PairCompare<T, U, Comparator>());
842  for (int i=0; i<size; i++) {
843  v1[i] = tmp[i].first;
844  v2[i] = tmp[i].second;
845  }
846 }
847 
848 
849 // returns whether a number is a denormal number or not
850 inline bool isDenormal(const float& x) {
851  return std::fpclassify(x) == FP_SUBNORMAL;
852 
853  // old version: works only for i386 and float
854  //const int& xbits = reinterpret_cast<const int&>(x);
855  //const int absMantissa = xbits & 0x007FFFFF;
856  //const int biasedExponent = xbits & 0x7F800000;
857  //return (biasedExponent == 0 && absMantissa != 0);
858 }
859 
860 // should always return a positive value, even when a/b is negative
861 template <typename T> T fmod(T a, T b) {
862  T q = floor(a/b);
863  return a - q*b;
864 }
865 
866 // returns the principal phase argument between [-PI,PI]
867 template <typename T> T princarg(T y) {
868  T x = essentia::fmod(y + M_PI, M_2PI);
869  //if (x < 0) x += M_2PI; // should be useless with our implementation of fmod
870  return x - M_PI;
871 }
872 
884 template <typename T>
885 void hist(const T* array, uint n, int* n_array, T* x_array, uint n_bins) {
886  T miny = *std::min_element(array, array+n);
887  T maxy = *std::max_element(array, array+n);
888 
889  // x contains the center of the bins
890  for (uint i=0; i<n_bins; i++) {
891  x_array[i] = (0.5 + i)*(maxy - miny)/n_bins + miny;
892  }
893 
894  // cutoff contains the boundaries between the bins
895  std::vector<T> cutoff(n_bins - 1);
896  for (uint i=0; i<n_bins-1; i++) {
897  cutoff[i] = (x_array[i] + x_array[i+1]) / 2.0;
898  }
899 
900  // algo: either build the cumulated histogram by 2-level
901  // for-loops, and then build the diff, or first
902  // sort the distribution and do it directly.
903  // 1st version: O(n^2) time, O(1) space
904  // 2nd version: O(n·log(n)) time, O(n) space
905 
906  // implementation of 2nd version
907  std::vector<T> dist(array, array+n);
908  std::sort(dist.begin(), dist.end());
909  uint current_cutoff_idx = 0;
910  T current_cutoff = cutoff[0];
911  for (uint i=0; i<n_bins; i++) n_array[i] = 0;
912 
913  for (uint i=0; i<n; i++) {
914  while (dist[i] > current_cutoff) {
915  // last case; skip the rest and fill in the last bin
916  if (current_cutoff_idx == n_bins-2) {
917  n_array[n_bins-1] = n-i; // fill in the last bin with what's left
918  i = n; // to jump out of the 2nd loop (the 'for' one)
919  n_array[n_bins-2]--; // to compensate for the last one that will be added before jumping out of the loop
920  break;
921  }
922  current_cutoff_idx++;
923  current_cutoff = cutoff[current_cutoff_idx];
924  }
925  n_array[current_cutoff_idx]++;
926  }
927 }
928 
929 
933 template <typename T>
934 void bincount(const std::vector<T>& input, std::vector<T>& output) {
935  output.clear();
936  output.resize( (int) ( std::max<Real>( input[argmax(input)], 0.) + 0.5 ) + 1);
937  uint index = 0;
938  for (uint i=0; i< input.size(); i++) {
939  index = int(std::max<Real>(input[i],0) + 0.5);
940  if (index < output.size() ) {
941  output[index] += 1.;
942  }
943  }
944 }
945 
946 
951 template <typename T>
952 std::vector<std::vector<T> > transpose(const std::vector<std::vector<T> >& m) {
953  if (m.empty()) return std::vector<std::vector<T> >();
954 
955  int nrows = m.size();
956  int ncols = m[0].size();
957  for (int i=1; i<nrows; i++) {
958  if ((int)m[i].size() != ncols) {
959  std::ostringstream ss;
960  ss <<"Trying to transpose a non rectangular matrix. Expecting dim2 = " << ncols
961  << " but got " << m[i].size() << ". Cannot transpose!";
962  throw EssentiaException(ss.str());
963  }
964  }
965 
966  std::vector<std::vector<T> > result(ncols, std::vector<Real>(nrows));
967  for (int i=0; i<nrows; i++) {
968  for (int j=0; j<ncols; j++) {
969  result[j][i] = m[i][j];
970  }
971  }
972 
973  return result;
974 }
975 
976 template <typename T>
978  if (m.dim1() == 0) return TNT::Array2D<T>();
979 
980  int nrows = m.dim1();
981  int ncols = m.dim2();
982 
983  TNT::Array2D<T> result(ncols, nrows);
984  for (int i=0; i<nrows; i++) {
985  for (int j=0; j<ncols; j++) {
986  result[j][i] = m[i][j];
987  }
988  }
989 
990  return result;
991 }
992 
993 inline std::string equivalentKey(const std::string key) {
994  if (key == "C")
995  return "C";
996 
997  if (key == "C#")
998  return "Db";
999 
1000  if (key == "Db")
1001  return "C#";
1002 
1003  if (key == "D")
1004  return "D";
1005 
1006  if (key == "D#")
1007  return "Eb";
1008 
1009  if (key == "Eb")
1010  return "D#";
1011 
1012  if (key == "E")
1013  return "E";
1014 
1015  if (key == "F")
1016  return "F";
1017 
1018  if (key == "F#")
1019  return "Gb";
1020 
1021  if (key == "Gb")
1022  return "F#";
1023 
1024  if (key == "G")
1025  return "G";
1026 
1027  if (key == "G#")
1028  return "Ab";
1029 
1030  if (key == "Ab")
1031  return "G#";
1032 
1033  if (key == "A")
1034  return "A";
1035 
1036  if (key == "A#")
1037  return "Bb";
1038 
1039  if (key == "Bb")
1040  return "A#";
1041 
1042  if (key == "B")
1043  return "B";
1044 
1045  return "";
1046 }
1047 
1048 
1053 template <typename T>
1054 void rotateChroma(std::vector<std::vector<T> >& inputMatrix, int oti) {
1055  if (inputMatrix.empty())
1056  throw EssentiaException("rotateChroma: trying to rotate an empty matrix");
1057  for (size_t i=0; i<inputMatrix.size(); i++) {
1058  std::rotate(inputMatrix[i].begin(), inputMatrix[i].end() - oti, inputMatrix[i].end());
1059  }
1060 }
1061 
1062 
1067 template <typename T>
1068 T dotProduct(const std::vector<T>& xArray, const std::vector<T>& yArray) {
1069  if (xArray.empty() || yArray.empty())
1070  throw EssentiaException("dotProduct: trying to calculate the dotProduct of empty arrays!");
1071  return std::inner_product(xArray.begin(), xArray.end(), yArray.begin(), 0.0);
1072 }
1073 
1074 
1079 template <typename T> T percentile(const std::vector<T>& array, Real qpercentile) {
1080  if (array.empty())
1081  throw EssentiaException("percentile: trying to calculate percentile of empty array");
1082 
1083  std::vector<T> sorted_array = array;
1084  // sort the array
1085  std::sort(sorted_array.begin(), sorted_array.end());
1086  qpercentile /= 100.;
1087 
1088  Real k;
1089  int sortArraySize = sorted_array.size();
1090  if (sortArraySize > 1) {
1091  k = (sortArraySize - 1) * qpercentile;
1092  }
1093  else {
1094  // to avoid zero value in arrays with single element
1095  k = sortArraySize * qpercentile;
1096  }
1097  // apply interpolation
1098  Real d0 = sorted_array[int(std::floor(k))] * (std::ceil(k) - k);
1099  Real d1 = sorted_array[int(std::ceil(k))] * (k - std::floor(k));
1100  return d0 + d1;
1101 }
1102 
1103 
1107 template <typename T> T covariance(const std::vector<T>& x, const T xMean, const std::vector<T>& y, const T yMean) {
1108  if (x.empty())
1109  throw EssentiaException("trying to calculate covariance of empty array");
1110  if (y.empty())
1111  throw EssentiaException("trying to calculate covariance of empty array");
1112  if (x.size() != y.size())
1113  throw EssentiaException("x and y should have the same size");
1114 
1115  T cov = (T) 0.0;
1116 
1117  for (uint i=0; i<x.size(); i++) {
1118  cov += (x[i] - xMean) * (y[i] - yMean);
1119  }
1120 
1121  return (T)(cov / (Real)x.size());
1122 }
1123 
1124 
1129 template <typename T> T pearsonCorrelationCoefficient(const std::vector<T>& x, const std::vector<T>& y) {
1130  if (x.empty())
1131  throw EssentiaException("trying to calculate covariance of empty array");
1132  if (y.empty())
1133  throw EssentiaException("trying to calculate covariance of empty array");
1134  if (x.size() != y.size())
1135  throw EssentiaException("x and y should have the same size");
1136 
1137  T xMean = mean(x);
1138  T yMean = mean(y);
1139 
1140  T cov = covariance(x, xMean, y, yMean);
1141 
1142  T xStddev = stddev(x, xMean);
1143  T yStddev = stddev(y, yMean);
1144 
1145  // When dealing with constants corraltion is 0 by convention.
1146  if ((xStddev == (T)0.0) || (xStddev == (T)0.0) || (xStddev == (T)0.0)) return (T) 0.0;
1147 
1148  T corr = cov / (xStddev * yStddev);
1149 
1150  // Numerical error can yield results slightly outside the analytical range [-1, 1].
1151  // Clipping the output is a cheap way to mantain this contrain.
1152  // Seen in https://github.com/numpy/numpy/blob/v1.15.0/numpy/lib/function_base.py#L2403-L2406
1153  return std::max(std::min(corr, (T)1.0), (T)-1.0);
1154 }
1155 
1156 
1163 template <typename T>
1164 void heavisideStepFunction(std::vector<std::vector<T> >& inputArray) {
1165  if (inputArray.empty())
1166  throw EssentiaException("heavisideStepFunction: found empty array as input!");
1167 
1168  for (size_t i=0; i<inputArray.size(); i++) {
1169  for (size_t j=0; j<inputArray[i].size(); j++) {
1170  // initialize all non negative elements as zero otherwise as one
1171  inputArray[i][j] = (inputArray[i][j] < 0) ? 0 : 1;
1172  }
1173  }
1174 }
1175 
1176 
1183 template <typename T>
1184 std::vector<std::vector<T> > pairwiseDistance(const std::vector<std::vector<T> >& m, const std::vector<std::vector<T> >& n) {
1185 
1186  if (m.empty() || n.empty())
1187  throw EssentiaException("pairwiseDistance: found empty array as input!");
1188 
1189  size_t mSize = m.size();
1190  size_t nSize = n.size();
1191  std::vector<std::vector<T> > pdist(mSize, std::vector<T>(nSize));
1192  for (size_t i=0; i<mSize; i++) {
1193  for (size_t j=0; j<nSize; j++) {
1194  Real item = dotProduct(m[i], m[i]) - 2*dotProduct(m[i], n[j]) + dotProduct(n[j], n[j]);
1195  pdist[i][j] = sqrt(item);
1196  }
1197  }
1198  if (pdist.empty())
1199  throw EssentiaException("pairwiseDistance: outputs an empty similarity matrix!");
1200  return pdist;
1201 }
1202 
1203 } // namespace essentia
1204 
1205 #endif // ESSENTIA_MATH_H
T sumSquare(const std::vector< T > array)
Definition: essentiamath.h:105
int dim1() const
Definition: tnt_array2d.h:231
T variance(const std::vector< T > &array, const T mean)
Definition: essentiamath.h:466
std::vector< T > derivative(const std::vector< T > &array)
Definition: essentiamath.h:807
T sum(const std::vector< T > &array, int start, int end)
Definition: essentiamath.h:116
void rectify(std::vector< T > &array)
Definition: essentiamath.h:394
Real lin2db(Real value)
Definition: essentiamath.h:543
std::string equivalentKey(const std::string key)
Definition: essentiamath.h:993
Real hz2mel(Real hz)
Definition: essentiamath.h:712
Real hz2cents(Real hz)
Definition: essentiamath.h:741
T median(const std::vector< T > &array)
Definition: essentiamath.h:372
T percentile(const std::vector< T > &array, Real qpercentile)
Definition: essentiamath.h:1079
Real hz2melSlaney(Real hz)
Definition: essentiamath.h:721
Real lin2log(Real value)
Definition: essentiamath.h:583
Real linear(Real input)
Definition: essentiamath.h:579
Real mel2hz(Real mel)
Definition: essentiamath.h:687
void heavisideStepFunction(std::vector< std::vector< T > > &inputArray)
Definition: essentiamath.h:1164
std::vector< std::vector< T > > pairwiseDistance(const std::vector< std::vector< T > > &m, const std::vector< std::vector< T > > &n)
Definition: essentiamath.h:1184
T stddev(const std::vector< T > &array, const T mean)
Definition: essentiamath.h:530
bool isDenormal(const float &x)
Definition: essentiamath.h:850
void bincount(const std::vector< T > &input, std::vector< T > &output)
Definition: essentiamath.h:934
#define M_2PI
Definition: essentiamath.h:41
#define SILENCE_CUTOFF
Definition: essentiamath.h:423
T round(const T value)
Definition: essentiamath.h:538
std::vector< T > varianceFrames(const std::vector< std::vector< T > > &frames)
Definition: essentiamath.h:264
T dotProduct(const std::vector< T > &xArray, const std::vector< T > &yArray)
Definition: essentiamath.h:1068
Real mel2hzSlaney(Real mel)
Definition: essentiamath.h:696
TNT::Array2D< T > meanMatrix(const std::vector< TNT::Array2D< T > * > &array)
Definition: essentiamath.h:167
Real db2pow(Real power)
Definition: essentiamath.h:563
T instantPower(const std::vector< T > &array)
Definition: essentiamath.h:409
Real amp2db(Real amplitude)
Definition: essentiamath.h:567
int argmin(const std::vector< Real > &input)
Definition: essentiamath.h:745
Definition: tnt_array2d.h:37
T fmod(T a, T b)
Definition: essentiamath.h:861
Real hz2hz(Real hz)
Definition: essentiamath.h:737
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:867
Real bark2hz(Real z)
Definition: essentiamath.h:673
TNT::Array2D< T > & matinit(TNT::Array2D< T > &A)
Definition: tnt2essentiautils.h:45
void normalizeAbs(std::vector< T > &array)
Definition: essentiamath.h:772
int dim2() const
Definition: tnt_array2d.h:234
Definition: essentiamath.h:820
#define DB_SILENCE_CUTOFF
Definition: essentiamath.h:424
float Real
Definition: types.h:69
Real db2lin(Real value)
Definition: essentiamath.h:551
std::vector< std::vector< T > > transpose(const std::vector< std::vector< T > > &m)
Definition: essentiamath.h:952
TNT::Array2D< T > varianceMatrix(const std::vector< TNT::Array2D< T > > &array, const TNT::Array2D< T > &mean)
Definition: essentiamath.h:434
T norm(const std::vector< T > &array)
Definition: essentiamath.h:88
void normalize(std::vector< T > &array)
Definition: essentiamath.h:759
std::vector< T > kurtosisFrames(const std::vector< std::vector< T > > &frames)
Definition: essentiamath.h:339
bool isPowerTwo(T n)
Definition: essentiamath.h:45
#define LOG_SILENCE_CUTOFF
Definition: essentiamath.h:425
Real pow2db(Real power)
Definition: essentiamath.h:555
std::vector< T > medianFrames(const std::vector< std::vector< T > > &frames, int beginIdx=0, int endIdx=-1)
Definition: essentiamath.h:225
std::vector< T > sumFrames(const std::vector< std::vector< T > > &frames)
Definition: essentiamath.h:290
Real hz2mel10(Real hz)
Definition: essentiamath.h:716
T energy(const std::vector< T > &array)
Definition: essentiamath.h:401
Definition: algorithm.h:28
Definition: types.h:77
Real barkCriticalBandwidth(Real z)
Definition: essentiamath.h:682
bool isSilent(const std::vector< T > &array)
Definition: essentiamath.h:428
unsigned int uint
Definition: types.h:49
void rotateChroma(std::vector< std::vector< T > > &inputMatrix, int oti)
Definition: essentiamath.h:1054
T pearsonCorrelationCoefficient(const std::vector< T > &x, const std::vector< T > &y)
Definition: essentiamath.h:1129
T nextPowerTwo(T n)
Definition: essentiamath.h:64
Real db2amp(Real amplitude)
Definition: essentiamath.h:575
T mean(const std::vector< T > &array, int start, int end)
Definition: essentiamath.h:142
Real hz2bark(Real f)
Definition: essentiamath.h:656
T kurtosis(const std::vector< T > &array, const T mean)
Definition: essentiamath.h:505
T skewness(const std::vector< T > &array, const T mean)
Definition: essentiamath.h:481
Comparator _cmp
Definition: essentiamath.h:821
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:823
void normalizeSum(std::vector< T > &array)
Definition: essentiamath.h:787
Real mel102hz(Real mel)
Definition: essentiamath.h:691
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:885
int argmax(const std::vector< Real > &input)
Definition: essentiamath.h:751
std::vector< T > skewnessFrames(const std::vector< std::vector< T > > &frames)
Definition: essentiamath.h:307
void sortpair(std::vector< T > &v1, std::vector< U > &v2)
Definition: essentiamath.h:833
T covariance(const std::vector< T > &x, const T xMean, const std::vector< T > &y, const T yMean)
Definition: essentiamath.h:1107