CLAM-Development
1.1
|
00001 /* 00002 * Copyright (c) 2004 MUSIC TECHNOLOGY GROUP (MTG) 00003 * UNIVERSITAT POMPEU FABRA 00004 * 00005 * 00006 * This program is free software; you can redistribute it and/or modify 00007 * it under the terms of the GNU General Public License as published by 00008 * the Free Software Foundation; either version 2 of the License, or 00009 * (at your option) any later version. 00010 * 00011 * This program is distributed in the hope that it will be useful, 00012 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00013 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00014 * GNU General Public License for more details. 00015 * 00016 * You should have received a copy of the GNU General Public License 00017 * along with this program; if not, write to the Free Software 00018 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00019 * 00020 */ 00021 00022 #ifndef _Stats_ 00023 #define _Stats_ 00024 00025 00026 #include "BasicOps.hxx" 00027 #include "Array.hxx" 00028 #include <algorithm> 00029 00030 namespace CLAM{ 00031 00032 00033 template <unsigned int x,unsigned int y> class GreaterThan 00034 { 00035 public: static StaticBool<(x>y)> mIs; 00036 }; 00037 00038 00039 template <unsigned int x,unsigned int y> StaticBool<(x>y)> GreaterThan<x,y>::mIs; 00040 00041 00050 template <typename T> 00051 class StatMemory 00052 { 00053 public: 00054 StatMemory() : mMemorized(false) {} 00055 const StatMemory & operator = (const T & value) 00056 { 00057 mMemorized=true; 00058 mMemory=value; 00059 00060 return *this; 00061 } 00062 bool HasValue() 00063 { 00064 return mMemorized; 00065 } 00066 void Reset() 00067 { 00068 mMemorized=false; 00069 } 00070 operator T() 00071 { 00072 CLAM_ASSERT(mMemorized,"Using a value that has not been memorized"); 00073 return mMemory; 00074 } 00075 private: 00076 T mMemory; 00077 bool mMemorized; 00078 }; 00079 00080 00089 template <bool abs=false,class T=TData, class U=TData,int initOrder=5> class StatsTmpl 00090 { 00091 00092 public: 00093 00094 /* @internal Only constructor available. We do not want a default constructor because then we could not be sure 00095 * that data is consisten and we would have to be constantly be doing checks.*/ 00096 00097 StatsTmpl(const Array<T>* data):mMoments(initOrder,5),mCentralMoments(initOrder,5),mCenterOfGravities(initOrder,5) 00098 { 00099 CLAM_ASSERT(data!=NULL,"Stats: A constructed array must be passed"); 00100 mData=data; 00101 /*we initialize moments up to initOrder-th order, if higher moments are asked for 00102 arrays are then resized*/ 00103 mMoments.SetSize(initOrder); 00104 mCentralMoments.SetSize(initOrder); 00105 mCenterOfGravities.SetSize(initOrder); 00106 for(unsigned i=0;i<initOrder;i++) 00107 { 00108 mMoments[i]=NULL; 00109 mCentralMoments[i]=NULL; 00110 mCenterOfGravities[i]= NULL; 00111 } 00112 InitMoment((O<initOrder>*)(0)); 00113 } 00114 ~StatsTmpl() 00115 { 00116 for (int i=0;i<mMoments.Size();i++) 00117 { 00118 if(mMoments[i]) delete mMoments[i]; 00119 } 00120 for (int i=0;i<mCentralMoments.Size();i++) 00121 { 00122 if(mCentralMoments[i]) delete mCentralMoments[i]; 00123 } 00124 for (int i=0;i<mCenterOfGravities.Size();i++) 00125 { 00126 if(mCenterOfGravities[i]) delete mCenterOfGravities[i]; 00127 } 00128 } 00129 00131 void SetArray(const Array<T>* data) 00132 { 00133 Reset(); 00134 mData=data; 00135 mCentroid.Reset(); 00136 } 00137 00143 template <int order> U GetMoment(const O<order>*) 00144 { 00145 return GetMoment((const O<order>*)(0),GreaterThan<order,initOrder>::mIs); 00146 } 00147 00149 template<int order> void GetMoments(Array<U>& moments, const O<order>*) 00150 { 00151 if(moments.Size()<order) 00152 { 00153 moments.Resize(order); 00154 moments.SetSize(order); 00155 } 00156 pTmpArray=&moments; 00157 GetChainedMoment((const O<order>*)(0)); 00158 pTmpArray=NULL; 00159 } 00160 00166 template <int order> U GetCentralMoment(const O<order>*) 00167 { 00168 return GetCentralMoment((const O<order>*)(0),GreaterThan<order,initOrder>::mIs); 00169 } 00170 00172 template<int order> void GetCentralMoments(Array<U>& centralMoments, 00173 const O<order>*) 00174 { 00175 if(centralMoments.Size()<order) 00176 { 00177 centralMoments.Resize(order); 00178 centralMoments.SetSize(order); 00179 } 00180 pTmpArray=¢ralMoments; 00181 GetChainedCentralMoment((const O<order>*)(0)); 00182 pTmpArray=NULL; 00183 } 00184 00190 template <int order> U GetCenterOfGravity(const O<order>*) 00191 { 00192 return GetCenterOfGravity((const O<order>*)(0),GreaterThan<order,initOrder>::mIs); 00193 } 00194 00196 template<int order> void GetCenterOfGravities(Array<U>& centerOfGravities, 00197 const O<order>*) 00198 { 00199 if(centerOfGravities.Size()<order) 00200 { 00201 centerOfGravities.Resize(order); 00202 centerOfGravities.SetSize(order); 00203 } 00204 pTmpArray=¢erOfGravities; 00205 GetChainedCenterOfGravity((const O<order>*)(0)); 00206 pTmpArray=NULL; 00207 } 00208 00216 U GetMean() 00217 { 00218 if (mData->Size()<=0) return U(.0); 00219 //FirstOrder* first; 00220 return GetMoment(FirstOrder); 00221 } 00222 00241 U GetCentroid() 00242 { 00243 // return GetCenterOfGravity(FirstOrder); 00244 if (mCentroid.HasValue()) return mCentroid; 00245 unsigned N = mData->Size(); 00246 U mean = GetMean(); 00247 if (mean < 1e-7 ) 00248 { 00249 mCentroid = U(N-1)/2; 00250 return mCentroid; 00251 } 00252 U centroid=0.0; 00253 for (unsigned i = 0; i < N; i++) 00254 { 00255 centroid += (abs?Abs((*mData)[i]):(*mData)[i]) * (i+1); 00256 } 00257 mCentroid=centroid/mean/U(N) - 1; 00258 return mCentroid; 00259 } 00260 00298 U GetSpread() 00299 { 00300 const unsigned N = mData->Size(); 00301 const Array<T> & data = *mData; 00302 const U centroid = GetCentroid(); 00303 00304 // Compute spectrum variance around centroid frequency 00305 TData variance = 0; 00306 TData sumMags = 0; 00307 for (unsigned i=0; i<N; i++) 00308 { 00309 U centroidDistance = i - centroid; 00310 centroidDistance *= centroidDistance; 00311 variance += centroidDistance * data[i]; 00312 sumMags += data[i]; 00313 } 00314 // NaN solving: Silence is like a plain distribution 00315 if (sumMags < 1e-14) return U(N+1) * (N-1) / 12; 00316 00317 return variance / sumMags; 00318 } 00319 00321 U GetStandardDeviation() 00322 { 00323 return mStdDev(*mData,GetCentralMomentFunctor<2>(),true); 00324 } 00325 00355 U GetSkew() 00356 { 00357 return mSkew(*mData,mStdDev,GetCentralMomentFunctor<3>(),true); 00358 } 00359 00383 U GetKurtosis() 00384 { 00385 return mKurtosis(*mData,GetCentralMomentFunctor<2>(),GetCentralMomentFunctor<4>(),true); 00386 } 00387 00398 U GetVariance() 00399 { 00400 return GetCentralMoment(SecondOrder); 00401 } 00402 00411 T GetEnergy() 00412 { 00413 return mEnergy(*mData); 00414 } 00415 00433 U GetGeometricMean() 00434 { 00435 return mGeometricMean(*mData); 00436 } 00437 00445 T GetRMS() 00446 { 00447 return mRMS(*mData); 00448 } 00449 00451 T GetMax() 00452 { 00453 return mMaxElement(*mData); 00454 } 00455 00457 T GetMin() 00458 { 00459 return mMinElement(*mData); 00460 } 00461 00495 U GetSlope() 00496 { 00497 // TODO: Sums where Y is used can be taken from Mean and Centroid 00498 00499 const TSize size = mData->Size(); 00500 00501 // \sum^{i=0}_{N-1}(x_i) 00502 // const TData sumY = GetMean()*size; 00503 // \sum^{i=0}_{N-1}(i x_i) 00504 // const TData sumXY = GetCentroid()*GetMean()*size; 00505 // \sum^{i=0}_{N-1}(i) 00506 // const TData sumX = (size-1)*size/2.0; 00507 // \sum^{i=0}_{N-1}(i^2) 00508 // const TData sumXX = (size-1)*(size)*(size+size-1)/6.0; 00509 00510 //TData num = size*sumXY - sumX*sumY; 00511 // = size Centroid Mean size - (size-1)(size)(size)Mean/2 00512 // = size^2 mean (Centroid - (size-1)/2) 00513 //num = size*size*GetMean()*(GetCentroid()-(size-1)/2.0); 00514 00515 // size*sumXX - sumX*sumX = 00516 // = size (size-1) size (size+size-1)/6 - (size-1)^2(size)^2/4 00517 // = size^2 ( (size-1)(size+size-1)/6 - (size-1)^2/4 ) 00518 // = size^2 (size-1)( (size+size-1)/6 - (size-1)/4 ) 00519 // = size^2 (size-1)( (4*size-2) - (3*size-3) )/12 00520 // = size^2 (size-1) (size+1)/12 00521 00522 //TData denum = (size*sumXX - sumX*sumX)*sumY; 00523 // = size mean size^2 (size-1) (size+1) / 12 00524 // = size^3 mean (size-1) (size+1) / 12 00525 //denum = size*size*size * GetMean() * (size-1) * (size+1) /12.0; 00526 00527 // return num/denum; 00528 // = size^2 mean (Centroid - (size-1)/2) / (size^3 mean (size-1) (size+1) / 12) 00529 // = (Centroid - (size-1)/2) / (size (size-1) (size+1) /12) 00530 // = ( 12*centroid - 6*size + 6 ) / ( size (size-1) (size+1) ) 00531 // = 6 (2*centroid - size + 1)) / ( size (size-1) (size+1) ) 00532 return 6*(2*GetCentroid() - size + 1) / (size * (size-1) * (size+1)); 00533 00534 } 00552 U GetFlatness() 00553 { 00554 U mean = GetMean(); 00555 U geometricMean = GetGeometricMean(); 00556 if (mean<1e-20) mean=TData(1e-20); 00557 if (geometricMean<1e-20 ) geometricMean=TData(1e-20); 00558 return geometricMean/mean; 00559 } 00560 00567 void Reset() 00568 { 00569 //Note: we keep previously allocated data, we just reset computations 00570 for (int i=0;i<mMoments.Size();i++) 00571 if(mMoments[i]!=NULL) mMoments[i]->Reset(); 00572 for (int i=0;i<mCentralMoments.Size();i++) 00573 if(mCentralMoments[i]!=NULL) mCentralMoments[i]->Reset(); 00574 for (int i=0;i<mCenterOfGravities.Size();i++) 00575 if(mCenterOfGravities[i]!=NULL) mCenterOfGravities[i]->Reset(); 00576 00577 mKurtosis.Reset(); 00578 mStdDev.Reset(); 00579 mSkew.Reset(); 00580 mEnergy.Reset(); 00581 mRMS.Reset(); 00582 mGeometricMean.Reset(); 00583 mMaxElement.Reset(); 00584 mMinElement.Reset(); 00585 mCentroid.Reset(); 00586 } 00587 00588 private: 00593 U GetTilt() 00594 { 00595 const Array<T>& Y = *mData; 00596 const TSize size = mData->Size(); 00597 const U m1 = GetMean(); 00598 00599 TData d1=0; 00600 TData d2=0; 00601 for (unsigned i=0;i<size;i++) 00602 { 00603 d1 += i/Y[i]; 00604 d2 += 1/Y[i]; 00605 } 00606 00607 // ti = m1/ai *(n - (d1/d2)) 00608 // SpecTilt = m1²/ti² * SUM[1/ai *(i-d1/d2)] 00609 00610 TData SumTi2 = 0; 00611 TData Tilt = 0; 00612 for (unsigned i=0;i<size;i++) 00613 { 00614 Tilt += (1/Y[i] *(i-d1/d2)); 00615 TData ti = m1/Y[i]*(i - (d1/d2)); 00616 SumTi2 += ti*ti; 00617 } 00618 00619 Tilt*= (m1*m1/SumTi2); 00620 return Tilt; 00621 } 00622 00624 template<int order> void InitMoment(const O<order>*) 00625 { 00626 if(mMoments[order-1]!=NULL) 00627 delete mMoments[order-1]; 00628 mMoments[order-1]=new Moment<order,abs,T,U>; 00629 if(mCentralMoments[order-1]!=NULL) 00630 delete mCentralMoments[order-1]; 00631 mCentralMoments[order-1]=new CentralMoment<order,abs,T,U>; 00632 if(mCenterOfGravities[order-1]!=NULL) 00633 delete mCenterOfGravities[order-1]; 00634 mCenterOfGravities[order-1]= new CenterOfGravity<order,abs,T,U>; 00635 InitMoment((O<order-1>*)(0)); 00636 } 00637 00639 void InitMoment(O<1>*) 00640 { 00641 mMoments[0]=new Moment<1,abs,T,U>; 00642 mCentralMoments[0]=new CentralMoment<1,abs,T,U>; 00643 mCenterOfGravities[0]= new CenterOfGravity<1,abs,T,U>; 00644 } 00645 00647 template<int order> U GetMoment(const O<order>*,StaticFalse&) 00648 { 00649 return (*(dynamic_cast<Moment<order,abs,T,U>*> (mMoments[order-1])))(*mData); 00650 } 00651 00653 template<int order> U GetMoment(const O<order>*,StaticTrue&) 00654 { 00655 if(order>mMoments.Size()) 00656 { 00657 int previousSize=mMoments.Size(); 00658 mMoments.Resize(order); 00659 mMoments.SetSize(order); 00660 for(int i=previousSize;i<order;i++) mMoments[i]=NULL; 00661 } 00662 00663 if(mMoments[order-1]==NULL) 00664 { 00665 mMoments[order-1]=new Moment<order,abs,T,U>; 00666 } 00667 //return GetMoment((const O<order>*)(0),StaticFalse()); 00668 return (*(dynamic_cast<Moment<order,abs,T,U>*> (mMoments[order-1])))(*mData); 00669 } 00670 00672 template<int order> void GetChainedMoment(const O<order>* ) 00673 { 00674 (*pTmpArray)[order-1]=GetMoment((const O<order>*)(0)); 00675 GetChainedMoment((O<order-1>*)(0)); 00676 } 00677 00679 void GetChainedMoment(O<1>* ) 00680 { 00681 (*pTmpArray)[0]=GetMoment((O<1>*)(0)); 00682 } 00683 00685 template<int order> U GetCentralMoment(const O<order>*,StaticFalse&) 00686 { 00687 CentralMoment<order,abs,T,U> & tmpMoment = GetCentralMomentFunctor<order>(); 00688 00689 //first we see if we already have corresponding Raw Moments up to the order demanded 00690 for(int i=0;i<order;i++) 00691 { 00692 //if we don't, we will have to compute them 00693 if(mMoments[i]==NULL) 00694 return tmpMoment(*mData); 00695 } 00696 00697 // if we do, we will use formula that relates Central Moments with Raw Moments 00698 return tmpMoment(*mData,mMoments); 00699 } 00700 00702 template<int order> U GetCentralMoment(const O<order>*,StaticTrue&) 00703 { 00704 if(order>mCentralMoments.Size()) 00705 { 00706 const int previousSize=mCentralMoments.Size(); 00707 mCentralMoments.Resize(order+1); 00708 mCentralMoments.SetSize(order+1); 00709 for(int i=previousSize; i<order; i++) mCentralMoments[i]=NULL; 00710 } 00711 if(mCentralMoments[order-1]==NULL) 00712 { 00713 mCentralMoments[order-1] = new CentralMoment<order,abs,T,U>; 00714 } 00715 00716 return GetCentralMoment((const O<order>*)(0),StaticFalse()); 00717 } 00718 00719 00720 00722 template<int order> void GetChainedCentralMoment(const O<order>* ) 00723 { 00724 (*pTmpArray)[order-1]=GetCentralMoment((const O<order>*)(0)); 00725 GetChainedCentralMoment((O<order-1>*)(0)); 00726 } 00727 00729 void GetChainedCentralMoment(O<1>* ) 00730 { 00731 (*pTmpArray)[0]=GetCentralMoment((O<1>*)(0)); 00732 } 00733 00735 template<int order> U GetCenterOfGravity(const O<order>*,StaticFalse& orderIsGreater) 00736 { 00737 return (*dynamic_cast<CenterOfGravity<order,abs,T,U>*> (mCenterOfGravities[order-1]))(*mData); 00738 } 00739 00741 template<int order> U GetCenterOfGravity(const O<order>*,StaticTrue& orderIsGreater) 00742 { 00743 if(order>mCenterOfGravities.Size()) 00744 { 00745 int previousSize=mCenterOfGravities.Size(); 00746 mCenterOfGravities.Resize(order); 00747 mCenterOfGravities.SetSize(order); 00748 for(int i=previousSize;i<order;i++) mCenterOfGravities[i]=NULL; 00749 } 00750 if(mCenterOfGravities[order-1]=NULL) 00751 { 00752 mCenterOfGravities[order-1]=new CenterOfGravity<order,abs,T,U>; 00753 } 00754 00755 return GetCenterOfGravity((const O<order>*)(0),StaticFalse()); 00756 } 00757 00759 template<int order> void GetChainedCenterOfGravity(const O<order>* ) 00760 { 00761 (*pTmpArray)[order-1]=GetCenterOfGravity((const O<order>*)(0)); 00762 GetChainedCenterOfGravity((O<order-1>*)(0)); 00763 } 00764 00766 void GetChainedCenterOfGravity(O<1>* ) 00767 { 00768 (*pTmpArray)[0]=GetCenterOfGravity((O<1>*)(0)); 00769 } 00770 00771 template <unsigned order> 00772 CentralMoment<order,abs,T,U> & GetCentralMomentFunctor() 00773 { 00774 CLAM_ASSERT( signed(order-1) < mCentralMoments.Size(), 00775 "Calling for a Central Moment order above the configured one"); 00776 00777 typedef CentralMoment<order,abs,T,U> CentralMomentN; 00778 const unsigned int position = order-1; 00779 if (!mCentralMoments[position]) 00780 mCentralMoments[position] = new CentralMomentN; 00781 return *dynamic_cast<CentralMomentN*>(mCentralMoments[position]); 00782 } 00783 00784 00785 Array<BaseMemOp*> mMoments; 00786 Array<BaseMemOp*> mCentralMoments; 00787 Array<BaseMemOp*> mCenterOfGravities; 00788 KurtosisTmpl<abs,T,U> mKurtosis; 00789 SkewTmpl<abs,T,U> mSkew; 00790 StandardDeviationTmpl<abs,T,U> mStdDev; 00791 EnergyTmpl<T> mEnergy; 00792 RMSTmpl<T> mRMS; 00793 GeometricMeanTmpl<T,U> mGeometricMean; 00794 ComplexMaxElement<abs,T> mMaxElement; 00795 ComplexMinElement<abs,T> mMinElement; 00796 StatMemory<U> mCentroid; 00797 00798 const Array<T>* mData; 00799 00801 Array<T>* pTmpArray; 00802 00803 }; 00804 00805 00806 00807 00808 00809 typedef StatsTmpl<> Stats; 00810 00811 00812 00813 };//namespace 00814 00815 00816 00817 #endif 00818