tswavelet.c

将提升小波用于目标跟踪。 小波提升的核心就是更新算法和预测算法,通过预测算法可以得到高频信息,而通过更新算子可以得到正确的低频信息.提升样式可以实现原位计算和整数提升,并且变换的中间结果是交织排列的

#include <stdio.h>
#include <math.h>
#include <string.h>
#define MAX_LAYER_NUM      3

void LWTrans( double *pData   , int nDataLength , 
             double *pWeight , int nWeightLength ,
             double *pBuff , int nLayer)
{
    int nC, nD , nE , i ,k,m;
    double *pHighFreqCoef, dfTemp;
    double *pWei;
    
    for( m=0; m<nLayer; m++ )
    {
        nC = (nDataLength+1)/2;
        nD = (nDataLength  )/2;
        nE = nC/2;
        
        for( i=0; i<nE; i++ )
            pBuff[i] = pData[i+i+1];
        for( i=1; i<nC; i++ )
            pData[i] = pData[i+i];
        
        pHighFreqCoef = &pData[nC];
        for( i=nD-1; i>=nE; i-- )
            pHighFreqCoef[i] = pData[i+i+1];
        for( i=0; i<nE; i++ )
            pHighFreqCoef[i] = pBuff[i];
        
        //transform
        pWei=pWeight;
        for( i=0; i<nWeightLength; i++ )
        {
            for( k=0; k<nC-1; k++ )
                pData[k] +=(pWei[0]*pHighFreqCoef[k] + pWei[1]*pHighFreqCoef[k+1]);
            pData[nC-1] +=pWei[0]*( pHighFreqCoef[nC-1] );
            
            pHighFreqCoef[0] +=pWei[2]*pData[0];		
            for( k=1; k<nD; k++ )
                pHighFreqCoef[k] +=(pWei[2]*pData[k] + pWei[3]*pData[k-1]);
            
            pWei = pWei+4;
        }
        
        for( k=0; k<nC; k++ )
            pData[k] *=pWei[0];
        dfTemp = 1.0/pWei[0];
        for( k=0; k<nD; k++ )
            pHighFreqCoef[k] *=dfTemp;

//        for( k=0; k<nC; k++)
//        {
//            printf("%4d %10.6f %10.6f\n", k, pData[k] , pHighFreqCoef[k]);
//        }

        nDataLength=nD;
        pData=pHighFreqCoef;
    }
}


void ILWTrans( double *pInputData   , int nDataLength , 
              double *pWeight , int nWeightLength ,
              double *pBuff , int nLayer)
{
    int nC, nD , k , i , j , m;
    double *pHighFreqCoef , *pCoef , *pData , dfTemp;
    int scalC[MAX_LAYER_NUM+1];
    int scalD[MAX_LAYER_NUM+1];
    
    k=0;
    m=nDataLength;
    for( i=1; i<=nLayer; i++ )
    {
        j=(m+1)/2;
        scalC[i]=j;
        m=m/2;
        scalD[i]=m;
        k=k+j;
    }

    for( m=nLayer; m>0; m--)
    {
        nD = scalD[m];
        nC = scalC[m];
        nDataLength=nC+nD;
        k  = k-nC;
        pData = pInputData+k;
        
        pHighFreqCoef = &pData[nC];
        dfTemp = pWeight[nWeightLength*4];
//        for( i=0; i<nD; i++ )
//        {
//            printf("%d ===%d %10.6f %10.6f\n", m, i, pHighFreqCoef[i], pData[i] );
//        }

        for( i=0; i<nD; i++ )
        {
            pHighFreqCoef[i] *=dfTemp;
        }
        dfTemp = 1.0/dfTemp;
        for( i=0; i<nC; i++ )
        {
            pData[i] *=dfTemp;
        }
        
        for( i=0; i<nWeightLength; i++ )
        {
            pCoef = &pWeight[(nWeightLength-1-i)*4];
            
            pHighFreqCoef[0] -=pCoef[2]*pData[0];		
            for( j=1; j<nD; j++ )
                pHighFreqCoef[j] -=(pCoef[2]*pData[j] + pCoef[3]*pData[j-1]);
            
            for( j=0; j<nC-1; j++ )
                pData[j] -=(pCoef[0]*pHighFreqCoef[j] + pCoef[1]*pHighFreqCoef[j+1]);
            pData[nC-1] -=( pCoef[0]* pHighFreqCoef[nC-1] );	
        }
        
        for( i=0; i<nC; i++ )
            pBuff[i+i] = pData[i];
        for( i=0; i<nD; i++ )
            pBuff[i+i+1] = pHighFreqCoef[i];
        memcpy( pData , pBuff , nDataLength*sizeof(double));
    }
//    for(i=0; i<nDataLength; i++ )
//    {
//        printf("??? %3d %f %f\n", i, pData[i] , pHighFreqCoef[i] );
//    }
}

void LWDenoise( double *pData   , int nDataLength , 
                int nLayer)
{
    int nC, i,m, nD;
    double dfSigma;

    nC = (nDataLength+1)/2;
    nD = nDataLength/2;
    for( i=0; i<nC; i++)
    {
        pData[i] = 0;
    }
    
    pData=pData+nC;
    nC= (nD+1)/2;
    nD= nD/2;
    m=nLayer-2; 
    while( m>0 )
    {
        dfSigma=0;
        for( i=0; i<nC; i++)
        {
            dfSigma=dfSigma+pData[i]*pData[i];
        }
        dfSigma = sqrt(dfSigma);
        dfSigma = (dfSigma/nC)*sqrt(2*log(1+nC));///log(2.0);
        for( i=0; i<nC; i++)
        {
            if( fabs(pData[i])>dfSigma )
            {
                //pData[i]=0;
            }
/*
            if( pData[i]>dfSigma )
            {
                pData[i]-=dfSigma;
            }
            else if( pData[i]<(-dfSigma) )
            {
                pData[i]+=dfSigma;
            }
*/
        }
        pData=pData+nC;
        nC= (nD+1)/2;
        nD=nD/2;
        m--;
    }
}