首页 > 代码库 > 小波 mallat 算法
小波 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 算法