小波 mallat 算法

算法要求：输入序列是大于滤波器长度的偶数列

确实可以通过编程的手段使算法适合所有的情况，但本文章的目的是展示mallat算法的过程，所以就一切从简了

// Mallat.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include "stdio.h"
/*mallat算法分解
* dSIn 输入的序列s，dH0尺度函数展开系数，dH1小波函数展开系数，dSOut输出低频部分，dDOut输出高频部分，
* nSIn_Len 输入序列的长度，nH_Len 滤波器的长度。
*/
int DwtFun(double *pdSIn,double *pdH0,double *pdH1,double *pdSOut,double *pdDOut,int nSIn_Len,int nH_Len)
{
int i,j,k;
//延拓后的Len是一个本体长度加一个滤波器长度
int nLen=nSIn_Len+2*nH_Len;

//建立滤波前的序列pdSArray，滤波后的序列pdSAOut低频部分，pdDAOut高频部分
double *pdSArray=new double[nLen];
double *pdSAOut=new double[nLen];
double *pdDAOut=new double[nLen];
//对称延拓
for(i=0;i<nLen;i++)
{
  if(i<nH_Len)
  {
   pdSArray[i]=pdSIn[nH_Len-i-1];
  }
  else if(i>=nH_Len+nSIn_Len)
  {
   pdSArray[i]=pdSIn[nH_Len+2*nSIn_Len-1-i];
  }
  else
  {
   pdSArray[i]=pdSIn[i-nH_Len];
  }
}

//求输出序列低频部分dSOut,高频部分dDOut.i->nLen,k->nH_Len
double dSTemp,dDTemp;
for(i=0;i<nLen;i++)
{
  dSTemp=0.0;
  dDTemp=0.0;
  for(k=0;k<nH_Len;k++)
  {
   if((i-k)<0)
    continue;
   else
   {
    //低频部分
    dSTemp+=pdH0[nH_Len-k-1]*pdSArray[i-k];

//高频部分
    dDTemp+=pdH1[nH_Len-k-1]*pdSArray[i-k];
   }
  }
  pdSAOut[i]=dSTemp;
  pdDAOut[i]=dDTemp;
}

//二抽取.先将pdSAOut前nH_Len长的一段舍弃,抽取偶数列
for(i=nH_Len,j=0;i<nLen;i+=2,j++)
{
  pdSOut[j]=pdSAOut[i+1];
  pdDOut[j]=pdDAOut[i+1];
}
//返回输出序列的长度
return j;

delete pdSArray;
pdSArray=NULL;
delete pdSAOut;
pdSAOut=NULL;
delete pdDAOut;
pdDAOut=NULL;
}
/*mallat 算法重构
* psSIn 输入的低频序列，pdDIn输入的高频序列，g0，g1重构滤波器，pdOut输出序列，nSInLen输入序列的长度
* nG_Len 滤波器长度
*/
int IDwtFun(double *pdSIn,double *pdDIn,double *pdG0,double *pdG1,double *pdOut,int nSInLen,int nG_Len)
{
int i,j,k;
//建立一个数列存放插入后的数列
int nTemp=2*nSInLen;
double *pdInSertS=new double[nTemp];
double *pdInSertD=new double[nTemp];
//二插入
j=0;
for(i=0;i<nTemp;i++)
{
  if(i%2==0)
  {
   pdInSertS[i]=0;
   pdInSertD[i]=0;
  }
  else
  {
   pdInSertS[i]=pdSIn[j];
   pdInSertD[i]=pdDIn[j];
   j++;
  }
}
//对称拓延
//创建一个nTemp+nG_Len长的数列
int nLen=nTemp+2*nG_Len;
double *pdSAIn=new double[nLen];
double *pdDAIn=new double[nLen];
for(i=0;i<nLen;i++)
{
  if(i<nG_Len)
  {
   pdSAIn[i]=pdInSertS[nG_Len-i-1];
   pdDAIn[i]=pdInSertD[nG_Len-i-1];
  }
  else if(i==nTemp+nG_Len)
  {
   pdSAIn[i]=0.0;
   pdDAIn[i]=0.0;
  }
  else if(i>nTemp+nG_Len)
  {
   pdSAIn[i]=pdInSertS[nG_Len+2*nTemp-i-1];
   pdDAIn[i]=pdInSertD[nG_Len+2*nTemp-i-1];
  }
  else
  {
   pdSAIn[i]=pdInSertS[i-nG_Len];
   pdDAIn[i]=pdInSertD[i-nG_Len];
  }
}
//用滤波器G0和G1对数列进行滤波
double *pdSAOut=new double[nLen];
double *pdDAOut=new double[nLen];
double dSTemp,dDTemp;
for(i=0;i<nLen;i++)
{
  dSTemp=0.0;
  dDTemp=0.0;
  for(k=0;k<nG_Len;k++)
  {
   if((i-k)<0)
    continue;
   else
   {
    //低频部分
    dSTemp+=pdG0[nG_Len-k-1]*pdSAIn[i-k];

//高频部分
    dDTemp+=pdG1[nG_Len-k-1]*pdDAIn[i-k];
   }
  }
  pdSAOut[i]=dSTemp;
  pdDAOut[i]=dDTemp;
}

//合并低频，高频
for(i=2*nG_Len-1,j=0;i<nLen;i++,j++)
{
pdOut[j]=pdSAOut[i]+pdDAOut[i];
}

return j;

delete pdInSertS;
pdInSertS=NULL;
delete pdInSertD;
pdInSertD=NULL;
delete pdSAIn;
pdSAIn=NULL;
delete pdDAIn;
pdDAIn=NULL;
delete pdSAOut;
pdSAOut=NULL;
delete pdDAOut;
pdDAOut=NULL;
}

int main(int argc, char* argv[])
{
int i;
//db4小波，已经取反 h0,h1是分解滤波器，g0,g1是重构滤波器
double dDb4h0[] = { 0.2303778133088964, 0.7148465705529154,
      0.6308807679398587, -0.0279837694168599,
        -0.1870348117190931, 0.0308413818355607,
            0.0328830116668852, -0.0105974017850690 };

double dDb4h1[] = { -0.0105974017850690 , -0.0328830116668852,
      0.0308413818355607 ,   0.1870348117190931,
      -0.0279837694168599 , -0.6308807679398587,
      0.7148465705529154 , -0.2303778133088964};

double dDb4g0[] = { -0.0105974017850690 , 0.0328830116668852,
      0.0308413818355607 , -0.1870348117190931,
      -0.0279837694168599 , 0.6308807679398587,
      0.7148465705529154 , 0.2303778133088964};

double dDb4g1[] = { -0.2303778133088964 , 0.7148465705529154,
      -0.6308807679398587 , -0.0279837694168599,
      0.1870348117190931 , 0.0308413818355607,
      -0.0328830116668852 , -0.0105974017850690};
//生成一个数列，本算法要求输入的数列是比滤波器长的偶数列
double a[]={1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0,11.0,12.0,13.0,14.0,15.0,16.0};
//double a[]={1.0,4.0,5.5,8.2,2.7,5.2,2.0,2.0,2.0,3.0,3.0,4.0,4.0,14.0,17.0,11.0};
//输出
double *pdS=new double[100];
double *pdD=new double[100];
double *pdOut=new double[100];

int l=DwtFun(a,dDb4h0,dDb4h1,pdS,pdD,16,8);
for(i=0;i<l-1;i++)
{
  printf("%f\t",pdS[i]);
  printf("\n");
}
printf("*********************\n");
for(i=0;i<l-1;i++)
{
  printf("%f\t",pdD[i]);
  printf("\n");
}
printf("*********************\n");
int v=IDwtFun(pdS,pdD,dDb4g0,dDb4g1,pdOut,11,8);
//i<v-nG_Len+1
for(i=0;i<v-7;i++)
{
  printf("%f\t",pdOut[i]);
  printf("\n");
}

delete []pdS;
pdS=NULL;
delete []pdD;
pdD=NULL;
delete []pdOut;
pdOut=NULL;
return 0;
}

秒客网

小波 mallat 算法

相关文章