Stats.hxx

Go to the documentation of this file.

/*
 * Copyright (c) 2004 MUSIC TECHNOLOGY GROUP (MTG)
 *                         UNIVERSITAT POMPEU FABRA
 *
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */

#ifndef _Stats_
#define _Stats_


#include "BasicOps.hxx"
#include "Array.hxx"
#include <algorithm>

namespace CLAM{


template <unsigned int x,unsigned int y> class GreaterThan
{
        public: static StaticBool<(x>y)> mIs;
};


template <unsigned int x,unsigned int y> StaticBool<(x>y)>  GreaterThan<x,y>::mIs;


template <typename T>
class StatMemory
{
public:
        StatMemory() : mMemorized(false) {}
        const StatMemory & operator = (const T & value)
        {
                mMemorized=true;
                mMemory=value;

                return *this;
        }
        bool HasValue()
        {
                return mMemorized;
        }
        void Reset()
        {
                mMemorized=false;
        }
        operator T()
        {
                CLAM_ASSERT(mMemorized,"Using a value that has not been memorized");
                return mMemory;
        }
private:
        T mMemory;
        bool mMemorized;
};


template <bool abs=false,class T=TData, class U=TData,int initOrder=5> class StatsTmpl
{

public:

/* @internal Only constructor available. We do not want a default constructor because then we could not be sure
 * that data is consisten and we would have to be constantly be doing checks.*/

        StatsTmpl(const Array<T>* data):mMoments(initOrder,5),mCentralMoments(initOrder,5),mCenterOfGravities(initOrder,5)
        {
                CLAM_ASSERT(data!=NULL,"Stats: A constructed array must be passed");
                mData=data;
                /*we initialize moments up to initOrder-th order, if higher moments are asked for
                arrays are then resized*/
                mMoments.SetSize(initOrder);
                mCentralMoments.SetSize(initOrder);
                mCenterOfGravities.SetSize(initOrder);
                for(unsigned i=0;i<initOrder;i++)
                {
                        mMoments[i]=NULL;
                        mCentralMoments[i]=NULL;
                        mCenterOfGravities[i]= NULL;
                }
                InitMoment((O<initOrder>*)(0));
        }
        ~StatsTmpl()
        {
                for (int i=0;i<mMoments.Size();i++)
                {
                        if(mMoments[i]) delete mMoments[i];
                }
                for (int i=0;i<mCentralMoments.Size();i++)
                {
                        if(mCentralMoments[i]) delete mCentralMoments[i];
                }
                for (int i=0;i<mCenterOfGravities.Size();i++)
                {
                        if(mCenterOfGravities[i]) delete mCenterOfGravities[i];
                }
        }

        void SetArray(const Array<T>* data)
        {
                Reset();
                mData=data;
                mCentroid.Reset();
        }

        template <int order> U GetMoment(const O<order>*)
        {
                return GetMoment((const O<order>*)(0),GreaterThan<order,initOrder>::mIs);
        }

        template<int order> void GetMoments(Array<U>& moments, const O<order>*)
        {
                if(moments.Size()<order)
                {
                        moments.Resize(order);
                        moments.SetSize(order);
                }
                pTmpArray=&moments;
                GetChainedMoment((const O<order>*)(0));
                pTmpArray=NULL;
        }

        template <int order> U GetCentralMoment(const O<order>*)
        {
                return GetCentralMoment((const O<order>*)(0),GreaterThan<order,initOrder>::mIs);
        }

        template<int order> void GetCentralMoments(Array<U>& centralMoments,
                const O<order>*)
        {
                if(centralMoments.Size()<order)
                {
                        centralMoments.Resize(order);
                        centralMoments.SetSize(order);
                }
                pTmpArray=&centralMoments;
                GetChainedCentralMoment((const O<order>*)(0));
                pTmpArray=NULL;
        }

        template <int order> U GetCenterOfGravity(const O<order>*)
        {
                return GetCenterOfGravity((const O<order>*)(0),GreaterThan<order,initOrder>::mIs);
        }

        template<int order> void GetCenterOfGravities(Array<U>& centerOfGravities,
                const O<order>*)
        {
                if(centerOfGravities.Size()<order)
                {
                        centerOfGravities.Resize(order);
                        centerOfGravities.SetSize(order);
                }
                pTmpArray=&centerOfGravities;
                GetChainedCenterOfGravity((const O<order>*)(0));
                pTmpArray=NULL;
        }

        U GetMean()
        {
                if (mData->Size()<=0) return U(.0);
                //FirstOrder* first;
                return GetMoment(FirstOrder);
        }

        U GetCentroid()
        {
//              return GetCenterOfGravity(FirstOrder);
                if (mCentroid.HasValue()) return mCentroid;
                unsigned N = mData->Size();
                U mean = GetMean();
                if (mean < 1e-7 ) 
                {
                        mCentroid = U(N-1)/2;
                        return mCentroid;
                }
                U centroid=0.0;
                for (unsigned i = 0; i < N; i++)
                {
                        centroid += (abs?Abs((*mData)[i]):(*mData)[i]) * (i+1);
                }
                mCentroid=centroid/mean/U(N) - 1;
                return mCentroid;
        }

        U GetSpread()
        {
                const unsigned N = mData->Size();
                const Array<T> & data = *mData;
                const U centroid = GetCentroid();

                // Compute spectrum variance around centroid frequency
                TData variance = 0;
                TData sumMags  = 0;
                for (unsigned i=0; i<N; i++)
                {
                        U centroidDistance = i - centroid;
                        centroidDistance *= centroidDistance;
                        variance += centroidDistance * data[i];
                        sumMags  += data[i];
                }
                // NaN solving: Silence is like a plain distribution
                if (sumMags < 1e-14) return U(N+1) * (N-1) / 12;

                return variance / sumMags;
        }

        U GetStandardDeviation()
        {
                return mStdDev(*mData,GetCentralMomentFunctor<2>(),true);
        }

        U GetSkew()
        {
                return mSkew(*mData,mStdDev,GetCentralMomentFunctor<3>(),true);
        }

        U GetKurtosis()
        {
                return mKurtosis(*mData,GetCentralMomentFunctor<2>(),GetCentralMomentFunctor<4>(),true);
        }

        U GetVariance()
        {
                return GetCentralMoment(SecondOrder);
        }

        T GetEnergy()
        {
                return mEnergy(*mData);
        }

        U GetGeometricMean()
        {
                return mGeometricMean(*mData);
        }

        T GetRMS()
        {
                return mRMS(*mData);
        }

        T GetMax()
        {
                return mMaxElement(*mData);
        }

        T GetMin()
        {
                return mMinElement(*mData);
        }

        U GetSlope()
        {
                // TODO: Sums where Y is used can be taken from Mean and Centroid

                const TSize size  = mData->Size();

                // \sum^{i=0}_{N-1}(x_i)
//              const TData sumY = GetMean()*size;
                // \sum^{i=0}_{N-1}(i x_i)
//              const TData sumXY = GetCentroid()*GetMean()*size;
                // \sum^{i=0}_{N-1}(i)
//              const TData sumX = (size-1)*size/2.0;
                // \sum^{i=0}_{N-1}(i^2)
//              const TData sumXX = (size-1)*(size)*(size+size-1)/6.0;

                //TData num = size*sumXY - sumX*sumY; 
                // = size Centroid Mean size - (size-1)(size)(size)Mean/2
                // = size^2 mean (Centroid - (size-1)/2)
                //num = size*size*GetMean()*(GetCentroid()-(size-1)/2.0);

                // size*sumXX - sumX*sumX =
                // = size (size-1) size (size+size-1)/6 - (size-1)^2(size)^2/4
                // = size^2 ( (size-1)(size+size-1)/6 - (size-1)^2/4 )
                // = size^2 (size-1)( (size+size-1)/6 - (size-1)/4 )
                // = size^2 (size-1)( (4*size-2) - (3*size-3) )/12
                // = size^2 (size-1) (size+1)/12

                //TData denum = (size*sumXX - sumX*sumX)*sumY;
                // = size mean size^2 (size-1) (size+1) / 12
                // = size^3 mean (size-1) (size+1) / 12
                //denum = size*size*size * GetMean() * (size-1) * (size+1) /12.0;

                // return num/denum;
                // = size^2 mean (Centroid - (size-1)/2) / (size^3 mean (size-1) (size+1) / 12)
                // = (Centroid - (size-1)/2) / (size (size-1) (size+1) /12)
                // = ( 12*centroid - 6*size + 6 ) / ( size (size-1) (size+1) )
                // = 6 (2*centroid - size + 1)) / ( size (size-1) (size+1) )
                return 6*(2*GetCentroid() - size + 1) / (size * (size-1) * (size+1));

        }
        U GetFlatness()
        {
                U mean = GetMean();
                U geometricMean = GetGeometricMean();
                if (mean<1e-20) mean=TData(1e-20);
                if (geometricMean<1e-20 ) geometricMean=TData(1e-20);
                return geometricMean/mean;
        }

        void Reset()
        {
                //Note: we keep previously allocated data, we just reset computations
                for (int i=0;i<mMoments.Size();i++)
                        if(mMoments[i]!=NULL) mMoments[i]->Reset();
                for (int i=0;i<mCentralMoments.Size();i++)
                        if(mCentralMoments[i]!=NULL) mCentralMoments[i]->Reset();
                for (int i=0;i<mCenterOfGravities.Size();i++)
                        if(mCenterOfGravities[i]!=NULL) mCenterOfGravities[i]->Reset();

                mKurtosis.Reset();
                mStdDev.Reset();
                mSkew.Reset();
                mEnergy.Reset();
                mRMS.Reset();
                mGeometricMean.Reset();
                mMaxElement.Reset();
                mMinElement.Reset();
                mCentroid.Reset();
        }

private:
        U GetTilt()
        {
                const Array<T>& Y = *mData;
                const TSize size  = mData->Size();
                const U m1 = GetMean();

                TData d1=0;
                TData d2=0;
                for (unsigned i=0;i<size;i++)
                {
                        d1 += i/Y[i];
                        d2 += 1/Y[i];
                }

                // ti = m1/ai *(n - (d1/d2))
                // SpecTilt = m1²/ti² * SUM[1/ai *(i-d1/d2)]

                TData SumTi2 = 0;
                TData Tilt = 0;
                for (unsigned i=0;i<size;i++) 
                {
                        Tilt += (1/Y[i] *(i-d1/d2));
                        TData ti = m1/Y[i]*(i - (d1/d2));
                        SumTi2 += ti*ti;
                }

                Tilt*= (m1*m1/SumTi2);
                return Tilt;
        }

        template<int order> void InitMoment(const O<order>*)
        {
                if(mMoments[order-1]!=NULL)
                        delete mMoments[order-1];
                mMoments[order-1]=new Moment<order,abs,T,U>;
                if(mCentralMoments[order-1]!=NULL)
                        delete mCentralMoments[order-1];
                mCentralMoments[order-1]=new CentralMoment<order,abs,T,U>;
                if(mCenterOfGravities[order-1]!=NULL)
                        delete mCenterOfGravities[order-1];
                mCenterOfGravities[order-1]= new CenterOfGravity<order,abs,T,U>;
                InitMoment((O<order-1>*)(0));
        }

        void InitMoment(O<1>*)
        {
                mMoments[0]=new Moment<1,abs,T,U>;
                mCentralMoments[0]=new CentralMoment<1,abs,T,U>;
                mCenterOfGravities[0]= new CenterOfGravity<1,abs,T,U>;
        }

        template<int order> U GetMoment(const O<order>*,StaticFalse&)
        {
                return (*(dynamic_cast<Moment<order,abs,T,U>*> (mMoments[order-1])))(*mData);
        }

        template<int order> U GetMoment(const O<order>*,StaticTrue&)
        {
                if(order>mMoments.Size())
                {
                        int previousSize=mMoments.Size();
                        mMoments.Resize(order);
                        mMoments.SetSize(order);
                        for(int i=previousSize;i<order;i++) mMoments[i]=NULL;
                }

                if(mMoments[order-1]==NULL)
                {
                        mMoments[order-1]=new Moment<order,abs,T,U>;
                }
                //return GetMoment((const O<order>*)(0),StaticFalse());
                return (*(dynamic_cast<Moment<order,abs,T,U>*> (mMoments[order-1])))(*mData);
        }

        template<int order> void GetChainedMoment(const O<order>* )
        {
                (*pTmpArray)[order-1]=GetMoment((const O<order>*)(0));
                GetChainedMoment((O<order-1>*)(0));
        }

        void GetChainedMoment(O<1>* )
        {
                (*pTmpArray)[0]=GetMoment((O<1>*)(0));
        }

        template<int order> U GetCentralMoment(const O<order>*,StaticFalse&)
        {
                CentralMoment<order,abs,T,U> & tmpMoment = GetCentralMomentFunctor<order>();

                //first we see if we already have corresponding Raw Moments up to the order demanded
                for(int i=0;i<order;i++)
                {
                        //if we don't, we will have to compute them
                        if(mMoments[i]==NULL)
                                return tmpMoment(*mData);
                }

                // if we do, we will use formula that relates Central Moments with Raw Moments
                return tmpMoment(*mData,mMoments);
        }

        template<int order> U GetCentralMoment(const O<order>*,StaticTrue&)
        {
                if(order>mCentralMoments.Size())
                {
                        const int previousSize=mCentralMoments.Size();
                        mCentralMoments.Resize(order+1);
                        mCentralMoments.SetSize(order+1);
                        for(int i=previousSize; i<order; i++) mCentralMoments[i]=NULL;
                }
                if(mCentralMoments[order-1]==NULL)
                {
                        mCentralMoments[order-1] = new CentralMoment<order,abs,T,U>;
                }

                return GetCentralMoment((const O<order>*)(0),StaticFalse());
        }



        template<int order> void GetChainedCentralMoment(const O<order>* )
        {
                (*pTmpArray)[order-1]=GetCentralMoment((const O<order>*)(0));
                GetChainedCentralMoment((O<order-1>*)(0));
        }

        void GetChainedCentralMoment(O<1>* )
        {
                (*pTmpArray)[0]=GetCentralMoment((O<1>*)(0));
        }

        template<int order> U GetCenterOfGravity(const O<order>*,StaticFalse& orderIsGreater)
        {
                return (*dynamic_cast<CenterOfGravity<order,abs,T,U>*> (mCenterOfGravities[order-1]))(*mData);
        }

        template<int order> U GetCenterOfGravity(const O<order>*,StaticTrue& orderIsGreater)
        {
                if(order>mCenterOfGravities.Size())
                {
                        int previousSize=mCenterOfGravities.Size();
                        mCenterOfGravities.Resize(order);
                        mCenterOfGravities.SetSize(order);
                        for(int i=previousSize;i<order;i++) mCenterOfGravities[i]=NULL;
                }
                if(mCenterOfGravities[order-1]=NULL)
                {
                        mCenterOfGravities[order-1]=new CenterOfGravity<order,abs,T,U>;
                }

                return GetCenterOfGravity((const O<order>*)(0),StaticFalse());
        }

        template<int order> void GetChainedCenterOfGravity(const O<order>* )
        {
                (*pTmpArray)[order-1]=GetCenterOfGravity((const O<order>*)(0));
                GetChainedCenterOfGravity((O<order-1>*)(0));
        }

        void GetChainedCenterOfGravity(O<1>* )
        {
                (*pTmpArray)[0]=GetCenterOfGravity((O<1>*)(0));
        }

        template <unsigned order>
        CentralMoment<order,abs,T,U> & GetCentralMomentFunctor()
        {
                CLAM_ASSERT( signed(order-1) < mCentralMoments.Size(),
                        "Calling for a Central Moment order above the configured one");

                typedef CentralMoment<order,abs,T,U> CentralMomentN;
                const unsigned int position = order-1;
                if (!mCentralMoments[position])
                        mCentralMoments[position] = new CentralMomentN;
                return *dynamic_cast<CentralMomentN*>(mCentralMoments[position]);
        }


        Array<BaseMemOp*> mMoments;
        Array<BaseMemOp*> mCentralMoments;
        Array<BaseMemOp*> mCenterOfGravities;
        KurtosisTmpl<abs,T,U> mKurtosis;
        SkewTmpl<abs,T,U> mSkew;
        StandardDeviationTmpl<abs,T,U> mStdDev;
        EnergyTmpl<T> mEnergy;
        RMSTmpl<T> mRMS;
        GeometricMeanTmpl<T,U> mGeometricMean;
        ComplexMaxElement<abs,T> mMaxElement;
        ComplexMinElement<abs,T> mMinElement;
        StatMemory<U> mCentroid;

        const Array<T>* mData;

        Array<T>* pTmpArray;

};





typedef StatsTmpl<> Stats;



};//namespace



#endif