libwavelet.c

linux系统下实现数据压缩的小波方法

#include "stdlib.h"
#include "stdio.h"
#include "math.h"

//---------------------------------------------------------------------------
//单层小波变换
//输入：
//      num:    数据点的个数
//      value:  数据点
//输出：
//      coef:   系数
idwt(int num,float value[20000],float coef[20000])
{
        //规范化haar小波基系数
        float appr[]={1.0/sqrt(2.0),1.0/sqrt(2.0)},detail[]={1.0/sqrt(2.0),-1.0/sqrt(2.0)};

        int i=0,half=0;

        half=num/2;
        i=0;
        while(i<half)
        {
             coef[i]=value[2*i]*appr[0]+value[2*i+1]*appr[1];
             coef[half+i]=value[2*i]*detail[0]+value[2*i+1]*detail[1];
             i++;
        }
}

//单层小波重构
//输入：
//      num:    系数的个数
//      coef:   系数
//输出：
//      value:   重构数据
idwt_rec(int num,float coef[20000],float value[20000])
{
        //规范化haar小波基系数
        float appr[]={1.0/sqrt(2.0),1.0/sqrt(2.0)},detail[]={1.0/sqrt(2.0),-1.0/sqrt(2.0)};

        int i=0,half=0;

        half=num/2;
        i=0;
        while(i<half)
        {
             value[2*i]=coef[i]*appr[0]+coef[half+i]*detail[0];
             value[2*i+1]=coef[i]*appr[1]+coef[half+i]*detail[1];
             i++;
        }
}

//小波变换
//输入：
//      lev:    变换层数(若指定层数超过最大层数将会被修改为最大层数)
//      num:    数据点个数
//      value:  数据点
//输出：
//      low_num:最低频部分的系数个数
//      total_num:系数个数
//      coef:   系数
dwt(int *lev,int num,float value[20000],int *low_num,int *total_num,float coef[20000])
{
        int max_lev=0,real_num=0,inum=0;
        int i=0,j=0;
        float tmp[20000];

        //求取最大分解层数
        max_lev=(int)(log10(num)/log10(2)+0.000001);
        if(num>(int)(pow(2,max_lev)+0.000001))
        {
                max_lev=max_lev+1;
        }

        //分解时的信号个数
        real_num=(int)(pow(2,max_lev)+0.000001);

        //以最后的信号值补齐信号
        for(i=num;i<real_num;i++)
        {
                value[i]=value[num-1];
        }

        //检查层数
        if(*lev>max_lev)
        {
                *lev=max_lev;
        }

        for(i=0;i<real_num;i++)
        {
                coef[i]=value[i];
        }

        //逐层进行小波变换
        for(i=0;i<*lev;i++)
        {
                inum=real_num/((int)(pow(2,i)+0.000001));
                for(j=0;j<inum;j++)
                {
                        tmp[j]=coef[j];                        
                }                
                idwt(inum,tmp,coef); //单层小波变换
        }
        *low_num=inum;
        *total_num=real_num;
}

//小波重构
//输入：
//
//      coef:   系数
//      lev:    层数
//      low_num:最低频部分的系数个数
//输出：
//      value:  重构数据
dwt_rec(float coef[20000],int lev,int low_num,float value[20000])
{
        int inum=0;
        int i=0,j=0;
        float tmp[20000];

        for(i=0;i<low_num;i++)
        {
                value[i]=coef[i];
        }

        //逐层进行小波重构
        for(i=0;i<lev;i++)
        {
                inum=low_num*((int)(pow(2,i)+0.000001));                
                for(j=0;j<inum;j++)
                {
                        if(j<inum/2)
                        {
                                tmp[j]=value[j];
                        }
                        else
                        {
                                tmp[j]=coef[j];
                        }                        
                }                
                idwt_rec(inum,tmp,value); //单层小波重构
        }
}

//均值
mymean(int num,float value[20000],float *mean)
{
        int i=0;

        *mean=0;
        for(i=0;i<num;i++)
        {
               *mean=*mean+value[i]/num;
        }
}

//标准差
mystdvari(int num,float value[20000],float *stdvari)
{
        int i=0;
        float mean=0,vari=0;

        for(i=0;i<num;i++)
        {
               mean=mean+value[i]/num;
        } 

        for(i=0;i<num;i++)
        {
               vari=vari+pow((value[i]-mean),2)/(num-1);
        }
        *stdvari=sqrt(vari);
}

//获取value绝对值排序对应的位置
//flag=0 升序；flag=1 降序；
//pos存储原始数据的位置
pos_abssort(int num,float value[20000],int flag,int pos[20000])
{
        float tmpvalue[20000],tmp;
        int i=0,j=0,k=0;

        for(i=0;i<num;i++)
        {
                pos[i]=i;
                tmpvalue[i]=value[i];
        }

        for(i=0;i<num;i++)
        {
                for(j=i+1;j<num;j++)
                {
                        if(flag==0)
                        {
                                if(fabs(tmpvalue[j])<fabs(tmpvalue[i]))
                                {
                                        k=pos[i];
                                        pos[i]=pos[j];
                                        pos[j]=k;

                                        tmp=tmpvalue[i];
                                        tmpvalue[i]=tmpvalue[j];
                                        tmpvalue[j]=tmp;
                                }
                        }
                        else
                        {
                                if(fabs(tmpvalue[j])>fabs(tmpvalue[i]))
                                {
                                        k=pos[i];
                                        pos[i]=pos[j];
                                        pos[j]=k;

                                        tmp=tmpvalue[i];
                                        tmpvalue[i]=tmpvalue[j];
                                        tmpvalue[j]=tmp;
                                }
                        }
                }
        }
}

//小波数据压缩
//压缩结果存储在resu中
compress(int num,float time[10000],float value[10000],char resu[64000],float *cpr)
{
        int lev;
        float e_ratio;
        char strlev[20],stre_ratio[20];

        float mean,stdvari;
        float tmpvalue[20000],coef[20000];
        double total_energe=0,comp_energe=0;
        int i=0,low_num,total_num,comp_num=0;
        int pos[20000];
        int sample=(int)((time[num-1]-time[0])/(num-1));
        char str[1000];

		
		if(rdini("WAVELET","LEVEL",strlev)==0)
		{
			WriteLog(-301);
		}
		lev=atoi(strlev);
		
		if(rdini("WAVELET","ENERGY",stre_ratio)==0)
		{
			WriteLog(-301);
		}
		e_ratio=atof(stre_ratio);		
		
        //数据标准化
        mymean(num,value,&mean);
        //mystdvari(num,value,&stdvari);
        for(i=0;i<num;i++)
        {
                tmpvalue[i]=(value[i]-mean);
        }

        //小波变换
        dwt(&lev,num,tmpvalue,&low_num,&total_num,coef);

        //基于能量保持率进行压缩
        pos_abssort(total_num,coef,1,pos);
        for(i=0;i<total_num;i++)
        {
                total_energe=total_energe+pow(coef[i],2);                
        }        
        for(i=0;i<total_num;i++)
        {
                if(comp_energe>=total_energe*e_ratio)
                {
                        break;
                }
                comp_energe=comp_energe+pow(coef[pos[i]],2);                
        }
        comp_num=i;

        //返回压缩结果
        strcpy(resu,"");
        sprintf(str,"%d:%d:%.4f:%.4f:%d:%d:%d:%d:",num,sample,mean,stdvari,comp_num,lev,total_num,low_num);
        strcat(resu,str);
        for(i=0;i<comp_num;i++)
        {
                sprintf(str,"%d:%.4f:",pos[i],coef[pos[i]]);
                strcat(resu,str);
        }

        *cpr=num/((float)comp_num);
}

//小波数据解压缩
//返回数据的个数、时间和值
decompress(char resu[64000],int *num,float time[10000],float value[10000])
{
        int sample=0,lev=0,comp_num=0,total_num=0,low_num=0;
        int i=0,j=0,k=0,flag=1;
        float mean=0,stdvari=0;
        float coef[20000],tmpcoef[20000];
        int pos[20000];

        //从压缩结果中获取信息
        sscanf(resu,"%d:%d:%f:%f:%d:%d:%d:%d",num,&sample,&mean,&stdvari,&comp_num,&lev,&total_num,&low_num);
        
        for(i=0;i<strlen(resu);i++)
        {
                if(resu[i]==':')
                {
                        k++;
                        if(k>=8)
                        {
                                flag++;
                                continue;
                        }
                }
                if(flag==2)
                {
                        sscanf(resu+i,"%d:%f",&(pos[j]),&(coef[j]));
                        flag=0;
                        j++;
                }
        }

        //还原小波系数
        for(i=0;i<total_num;i++)
        {
                tmpcoef[i]=0;
        }
        
        for(i=0;i<comp_num;i++)
        {
                tmpcoef[pos[i]]=coef[i];
        } 

        //小波重构        
        dwt_rec(tmpcoef,lev,low_num,value);

        //返回解压缩结果
        for(i=0;i<total_num;i++)
        {
                value[i]=value[i]+mean; 
                time[i]=sample*i;
        }
}