CLAM-Development  1.4.0
Stats.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2004 MUSIC TECHNOLOGY GROUP (MTG)
3  * UNIVERSITAT POMPEU FABRA
4  *
5  *
6  * This program is free software; you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation; either version 2 of the License, or
9  * (at your option) any later version.
10  *
11  * This program is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with this program; if not, write to the Free Software
18  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19  *
20  */
21 
22 #ifndef _Stats_
23 #define _Stats_
24 
25 
26 #include "BasicOps.hxx"
27 #include "Array.hxx"
28 #include <algorithm>
29 
30 namespace CLAM{
31 
32 
33 template <unsigned int x,unsigned int y> class GreaterThan
34 {
35  public: static StaticBool<(x>y)> mIs;
36 };
37 
38 
39 template <unsigned int x,unsigned int y> StaticBool<(x>y)> GreaterThan<x,y>::mIs;
40 
41 
50 template <typename T>
52 {
53 public:
54  StatMemory() : mMemorized(false) {}
55  const StatMemory & operator = (const T & value)
56  {
57  mMemorized=true;
58  mMemory=value;
59 
60  return *this;
61  }
62  bool HasValue()
63  {
64  return mMemorized;
65  }
66  void Reset()
67  {
68  mMemorized=false;
69  }
70  operator T()
71  {
72  CLAM_ASSERT(mMemorized,"Using a value that has not been memorized");
73  return mMemory;
74  }
75 private:
76  T mMemory;
77  bool mMemorized;
78 };
79 
80 
89 template <bool abs=false,class T=TData, class U=TData,int initOrder=5> class StatsTmpl
90 {
91 
92 public:
93 
94 /* @internal Only constructor available. We do not want a default constructor because then we could not be sure
95  * that data is consisten and we would have to be constantly be doing checks.*/
96 
97  StatsTmpl(const Array<T>* data):mMoments(initOrder,5),mCentralMoments(initOrder,5),mCenterOfGravities(initOrder,5)
98  {
99  CLAM_ASSERT(data!=NULL,"Stats: A constructed array must be passed");
100  mData=data;
101  /*we initialize moments up to initOrder-th order, if higher moments are asked for
102  arrays are then resized*/
103  mMoments.SetSize(initOrder);
104  mCentralMoments.SetSize(initOrder);
105  mCenterOfGravities.SetSize(initOrder);
106  for(unsigned i=0;i<initOrder;i++)
107  {
108  mMoments[i]=NULL;
109  mCentralMoments[i]=NULL;
110  mCenterOfGravities[i]= NULL;
111  }
112  InitMoment((O<initOrder>*)(0));
113  }
115  {
116  for (int i=0;i<mMoments.Size();i++)
117  {
118  if(mMoments[i]) delete mMoments[i];
119  }
120  for (int i=0;i<mCentralMoments.Size();i++)
121  {
122  if(mCentralMoments[i]) delete mCentralMoments[i];
123  }
124  for (int i=0;i<mCenterOfGravities.Size();i++)
125  {
126  if(mCenterOfGravities[i]) delete mCenterOfGravities[i];
127  }
128  }
129 
131  void SetArray(const Array<T>* data)
132  {
133  Reset();
134  mData=data;
135  mCentroid.Reset();
136  }
137 
143  template <int order> U GetMoment(const O<order>*)
144  {
146  }
147 
149  template<int order> void GetMoments(Array<U>& moments, const O<order>*)
150  {
151  if(moments.Size()<order)
152  {
153  moments.Resize(order);
154  moments.SetSize(order);
155  }
156  pTmpArray=&moments;
157  GetChainedMoment((const O<order>*)(0));
158  pTmpArray=NULL;
159  }
160 
166  template <int order> U GetCentralMoment(const O<order>*)
167  {
169  }
170 
172  template<int order> void GetCentralMoments(Array<U>& centralMoments,
173  const O<order>*)
174  {
175  if(centralMoments.Size()<order)
176  {
177  centralMoments.Resize(order);
178  centralMoments.SetSize(order);
179  }
180  pTmpArray=&centralMoments;
181  GetChainedCentralMoment((const O<order>*)(0));
182  pTmpArray=NULL;
183  }
184 
190  template <int order> U GetCenterOfGravity(const O<order>*)
191  {
193  }
194 
196  template<int order> void GetCenterOfGravities(Array<U>& centerOfGravities,
197  const O<order>*)
198  {
199  if(centerOfGravities.Size()<order)
200  {
201  centerOfGravities.Resize(order);
202  centerOfGravities.SetSize(order);
203  }
204  pTmpArray=&centerOfGravities;
205  GetChainedCenterOfGravity((const O<order>*)(0));
206  pTmpArray=NULL;
207  }
208 
217  {
218  if (mData->Size()<=0) return U(.0);
219  //FirstOrder* first;
220  return GetMoment(FirstOrder);
221  }
222 
242  {
243 // return GetCenterOfGravity(FirstOrder);
244  if (mCentroid.HasValue()) return mCentroid;
245  unsigned N = mData->Size();
246  U mean = GetMean();
247  if (mean < 1e-7 )
248  {
249  mCentroid = U(N-1)/2;
250  return mCentroid;
251  }
252  U centroid=0.0;
253  for (unsigned i = 0; i < N; i++)
254  {
255  centroid += (abs?Abs((*mData)[i]):(*mData)[i]) * (i+1);
256  }
257  mCentroid=centroid/mean/U(N) - 1;
258  return mCentroid;
259  }
260 
299  {
300  const unsigned N = mData->Size();
301  const Array<T> & data = *mData;
302  const U centroid = GetCentroid();
303 
304  // Compute spectrum variance around centroid frequency
305  TData variance = 0;
306  TData sumMags = 0;
307  for (unsigned i=0; i<N; i++)
308  {
309  U centroidDistance = i - centroid;
310  centroidDistance *= centroidDistance;
311  variance += centroidDistance * data[i];
312  sumMags += data[i];
313  }
314  // NaN solving: Silence is like a plain distribution
315  if (sumMags < 1e-14) return U(N+1) * (N-1) / 12;
316 
317  return variance / sumMags;
318  }
319 
322  {
323  return mStdDev(*mData,GetCentralMomentFunctor<2>(),true);
324  }
325 
356  {
357  return mSkew(*mData,mStdDev,GetCentralMomentFunctor<3>(),true);
358  }
359 
384  {
385  return mKurtosis(*mData,GetCentralMomentFunctor<2>(),GetCentralMomentFunctor<4>(),true);
386  }
387 
399  {
401  }
402 
412  {
413  return mEnergy(*mData);
414  }
415 
434  {
435  return mGeometricMean(*mData);
436  }
437 
445  T GetRMS()
446  {
447  return mRMS(*mData);
448  }
449 
451  T GetMax()
452  {
453  return mMaxElement(*mData);
454  }
455 
457  T GetMin()
458  {
459  return mMinElement(*mData);
460  }
461 
496  {
497  // TODO: Sums where Y is used can be taken from Mean and Centroid
498 
499  const TSize size = mData->Size();
500 
501  // \sum^{i=0}_{N-1}(x_i)
502 // const TData sumY = GetMean()*size;
503  // \sum^{i=0}_{N-1}(i x_i)
504 // const TData sumXY = GetCentroid()*GetMean()*size;
505  // \sum^{i=0}_{N-1}(i)
506 // const TData sumX = (size-1)*size/2.0;
507  // \sum^{i=0}_{N-1}(i^2)
508 // const TData sumXX = (size-1)*(size)*(size+size-1)/6.0;
509 
510  //TData num = size*sumXY - sumX*sumY;
511  // = size Centroid Mean size - (size-1)(size)(size)Mean/2
512  // = size^2 mean (Centroid - (size-1)/2)
513  //num = size*size*GetMean()*(GetCentroid()-(size-1)/2.0);
514 
515  // size*sumXX - sumX*sumX =
516  // = size (size-1) size (size+size-1)/6 - (size-1)^2(size)^2/4
517  // = size^2 ( (size-1)(size+size-1)/6 - (size-1)^2/4 )
518  // = size^2 (size-1)( (size+size-1)/6 - (size-1)/4 )
519  // = size^2 (size-1)( (4*size-2) - (3*size-3) )/12
520  // = size^2 (size-1) (size+1)/12
521 
522  //TData denum = (size*sumXX - sumX*sumX)*sumY;
523  // = size mean size^2 (size-1) (size+1) / 12
524  // = size^3 mean (size-1) (size+1) / 12
525  //denum = size*size*size * GetMean() * (size-1) * (size+1) /12.0;
526 
527  // return num/denum;
528  // = size^2 mean (Centroid - (size-1)/2) / (size^3 mean (size-1) (size+1) / 12)
529  // = (Centroid - (size-1)/2) / (size (size-1) (size+1) /12)
530  // = ( 12*centroid - 6*size + 6 ) / ( size (size-1) (size+1) )
531  // = 6 (2*centroid - size + 1)) / ( size (size-1) (size+1) )
532  return 6*(2*GetCentroid() - size + 1) / (size * (size-1) * (size+1));
533 
534  }
553  {
554  U mean = GetMean();
555  U geometricMean = GetGeometricMean();
556  if (mean<1e-20) mean=TData(1e-20);
557  if (geometricMean<1e-20 ) geometricMean=TData(1e-20);
558  return geometricMean/mean;
559  }
560 
567  void Reset()
568  {
569  //Note: we keep previously allocated data, we just reset computations
570  for (int i=0;i<mMoments.Size();i++)
571  if(mMoments[i]!=NULL) mMoments[i]->Reset();
572  for (int i=0;i<mCentralMoments.Size();i++)
573  if(mCentralMoments[i]!=NULL) mCentralMoments[i]->Reset();
574  for (int i=0;i<mCenterOfGravities.Size();i++)
575  if(mCenterOfGravities[i]!=NULL) mCenterOfGravities[i]->Reset();
576 
577  mKurtosis.Reset();
578  mStdDev.Reset();
579  mSkew.Reset();
580  mEnergy.Reset();
581  mRMS.Reset();
582  mGeometricMean.Reset();
583  mMaxElement.Reset();
584  mMinElement.Reset();
585  mCentroid.Reset();
586  }
587 
588 private:
593  U GetTilt()
594  {
595  const Array<T>& Y = *mData;
596  const TSize size = mData->Size();
597  const U m1 = GetMean();
598 
599  TData d1=0;
600  TData d2=0;
601  for (unsigned i=0;i<size;i++)
602  {
603  d1 += i/Y[i];
604  d2 += 1/Y[i];
605  }
606 
607  // ti = m1/ai *(n - (d1/d2))
608  // SpecTilt = m1²/ti² * SUM[1/ai *(i-d1/d2)]
609 
610  TData SumTi2 = 0;
611  TData Tilt = 0;
612  for (unsigned i=0;i<size;i++)
613  {
614  Tilt += (1/Y[i] *(i-d1/d2));
615  TData ti = m1/Y[i]*(i - (d1/d2));
616  SumTi2 += ti*ti;
617  }
618 
619  Tilt*= (m1*m1/SumTi2);
620  return Tilt;
621  }
622 
624  template<int order> void InitMoment(const O<order>*)
625  {
626  if(mMoments[order-1]!=NULL)
627  delete mMoments[order-1];
628  mMoments[order-1]=new Moment<order,abs,T,U>;
629  if(mCentralMoments[order-1]!=NULL)
630  delete mCentralMoments[order-1];
631  mCentralMoments[order-1]=new CentralMoment<order,abs,T,U>;
632  if(mCenterOfGravities[order-1]!=NULL)
633  delete mCenterOfGravities[order-1];
634  mCenterOfGravities[order-1]= new CenterOfGravity<order,abs,T,U>;
635  InitMoment((O<order-1>*)(0));
636  }
637 
639  void InitMoment(O<1>*)
640  {
641  mMoments[0]=new Moment<1,abs,T,U>;
642  mCentralMoments[0]=new CentralMoment<1,abs,T,U>;
643  mCenterOfGravities[0]= new CenterOfGravity<1,abs,T,U>;
644  }
645 
647  template<int order> U GetMoment(const O<order>*,StaticFalse&)
648  {
649  return (*(dynamic_cast<Moment<order,abs,T,U>*> (mMoments[order-1])))(*mData);
650  }
651 
653  template<int order> U GetMoment(const O<order>*,StaticTrue&)
654  {
655  if(order>mMoments.Size())
656  {
657  int previousSize=mMoments.Size();
658  mMoments.Resize(order);
659  mMoments.SetSize(order);
660  for(int i=previousSize;i<order;i++) mMoments[i]=NULL;
661  }
662 
663  if(mMoments[order-1]==NULL)
664  {
665  mMoments[order-1]=new Moment<order,abs,T,U>;
666  }
667  //return GetMoment((const O<order>*)(0),StaticFalse());
668  return (*(dynamic_cast<Moment<order,abs,T,U>*> (mMoments[order-1])))(*mData);
669  }
670 
672  template<int order> void GetChainedMoment(const O<order>* )
673  {
674  (*pTmpArray)[order-1]=GetMoment((const O<order>*)(0));
675  GetChainedMoment((O<order-1>*)(0));
676  }
677 
679  void GetChainedMoment(O<1>* )
680  {
681  (*pTmpArray)[0]=GetMoment((O<1>*)(0));
682  }
683 
685  template<int order> U GetCentralMoment(const O<order>*,StaticFalse&)
686  {
687  CentralMoment<order,abs,T,U> & tmpMoment = GetCentralMomentFunctor<order>();
688 
689  //first we see if we already have corresponding Raw Moments up to the order demanded
690  for(int i=0;i<order;i++)
691  {
692  //if we don't, we will have to compute them
693  if(mMoments[i]==NULL)
694  return tmpMoment(*mData);
695  }
696 
697  // if we do, we will use formula that relates Central Moments with Raw Moments
698  return tmpMoment(*mData,mMoments);
699  }
700 
702  template<int order> U GetCentralMoment(const O<order>*,StaticTrue&)
703  {
704  if(order>mCentralMoments.Size())
705  {
706  const int previousSize=mCentralMoments.Size();
707  mCentralMoments.Resize(order+1);
708  mCentralMoments.SetSize(order+1);
709  for(int i=previousSize; i<order; i++) mCentralMoments[i]=NULL;
710  }
711  if(mCentralMoments[order-1]==NULL)
712  {
713  mCentralMoments[order-1] = new CentralMoment<order,abs,T,U>;
714  }
715 
716  return GetCentralMoment((const O<order>*)(0),StaticFalse());
717  }
718 
719 
720 
722  template<int order> void GetChainedCentralMoment(const O<order>* )
723  {
724  (*pTmpArray)[order-1]=GetCentralMoment((const O<order>*)(0));
725  GetChainedCentralMoment((O<order-1>*)(0));
726  }
727 
729  void GetChainedCentralMoment(O<1>* )
730  {
731  (*pTmpArray)[0]=GetCentralMoment((O<1>*)(0));
732  }
733 
735  template<int order> U GetCenterOfGravity(const O<order>*,StaticFalse& orderIsGreater)
736  {
737  return (*dynamic_cast<CenterOfGravity<order,abs,T,U>*> (mCenterOfGravities[order-1]))(*mData);
738  }
739 
741  template<int order> U GetCenterOfGravity(const O<order>*,StaticTrue& orderIsGreater)
742  {
743  if(order>mCenterOfGravities.Size())
744  {
745  int previousSize=mCenterOfGravities.Size();
746  mCenterOfGravities.Resize(order);
747  mCenterOfGravities.SetSize(order);
748  for(int i=previousSize;i<order;i++) mCenterOfGravities[i]=NULL;
749  }
750  if(mCenterOfGravities[order-1]=NULL)
751  {
752  mCenterOfGravities[order-1]=new CenterOfGravity<order,abs,T,U>;
753  }
754 
755  return GetCenterOfGravity((const O<order>*)(0),StaticFalse());
756  }
757 
759  template<int order> void GetChainedCenterOfGravity(const O<order>* )
760  {
761  (*pTmpArray)[order-1]=GetCenterOfGravity((const O<order>*)(0));
762  GetChainedCenterOfGravity((O<order-1>*)(0));
763  }
764 
766  void GetChainedCenterOfGravity(O<1>* )
767  {
768  (*pTmpArray)[0]=GetCenterOfGravity((O<1>*)(0));
769  }
770 
771  template <unsigned order>
772  CentralMoment<order,abs,T,U> & GetCentralMomentFunctor()
773  {
774  CLAM_ASSERT( signed(order-1) < mCentralMoments.Size(),
775  "Calling for a Central Moment order above the configured one");
776 
777  typedef CentralMoment<order,abs,T,U> CentralMomentN;
778  const unsigned int position = order-1;
779  if (!mCentralMoments[position])
780  mCentralMoments[position] = new CentralMomentN;
781  return *dynamic_cast<CentralMomentN*>(mCentralMoments[position]);
782  }
783 
784 
785  Array<BaseMemOp*> mMoments;
786  Array<BaseMemOp*> mCentralMoments;
787  Array<BaseMemOp*> mCenterOfGravities;
788  KurtosisTmpl<abs,T,U> mKurtosis;
789  SkewTmpl<abs,T,U> mSkew;
790  StandardDeviationTmpl<abs,T,U> mStdDev;
791  EnergyTmpl<T> mEnergy;
792  RMSTmpl<T> mRMS;
793  GeometricMeanTmpl<T,U> mGeometricMean;
794  ComplexMaxElement<abs,T> mMaxElement;
795  ComplexMinElement<abs,T> mMinElement;
796  StatMemory<U> mCentroid;
797 
798  const Array<T>* mData;
799 
801  Array<T>* pTmpArray;
802 
803 };
804 
805 
806 
807 
808 
810 
811 
812 
813 };//namespace
814 
815 
816 
817 #endif
818