小波 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;
}

时间： 2024-09-29 13:35:18

小波 mallat 算法的相关文章

图像融合（六）-- 小波融合

基于小波的融合(wavelet) 小波变换的固有特性使其在图像处理中有如下优点:完善的重构能力,保证信号在分解过程中没有信息损失和冗余信息:把图像分解成平均图像和细节图像的组合,分别代表了图像的不同结构,因此容易提取原始图像的结构信息和细节信息:小波分析提供了与人类视觉系统方向相吻合的选择性图像. 离散小波变换(Discrete Wavelet Transform, DWT).DWT的函数基由一个称为母小波或分析小波的单一函数通过膨胀和平移获得.因而,DWT同时具有时域和频域分析能力,与一般的金

小波系数

1. 求小波变化系数时a b怎么取? 小波变换的概念是由法国从事石油信号处理的工程师J.Morlet在1974年首先提出的,通过物理的直观和信号处理的实际需要经验的建立了反演公式,当时未能得到数学家的认可.正如1807年法国的热学工程师J.B.J.Fourier提出任一函数都能展开成三角函数的无穷级数的创新概念未能得到数学家J.L.Lagrange,P.S.Laplace以及A.M.Legendre的认可一样.幸运的是,早在七十年代,A.Calderon表示定理的发现.Hardy空间的原子分解和

第一代程序猿王小波

Managershare:据我所知,最文艺的群体并不是文科生,而是理科生. 多数人知道王小波是小说家,部分人分不清财经作家吴晓波和小说家王小波是不是一回事儿.却很少有人知道王小波可以算的上中国早期的程序员,在 90 年代初的时候因为国内应用软件缺乏,爱捣鼓东西的王小波利用闲暇时间学习了汇编和C语言,编了中文编辑器和输入法.中文编辑器和输入法任何一个都是大牛级的 GEEK 才会去尝试的东西,比如求伯君.王小波通过卖软件还挣了些钱,当时很多中观村的老板要拉他入伙,当然写代码这种来钱快的活对屌丝王小波

好吧，左小波出山了！

我,还是一个不懂世事的毛头小子,第一次写博.万事开头难,没事咱慢慢来.咳,练文笔吗.我觉得写东西最锻炼逻辑思维,我是一个不善于表达的人,可能是程序员的通病,但你看看人家王小波,八九十年代的作家兼职程序员,可不可怕!所以那都不是借口,同学们表达多重要啊!不会表达,你就不会撩妹.看看王小波长成那德行了,几句情诗美的我银河奶奶那么开心.我是一个看到别人优点就想学的人,因为不学习俺就头疼. 最近我不知道为什么陷入无法自拔的自嗨状态,很有可能是自恋又犯了,可那又有啥招,我总被自己帅醒这是真的.当我也纠结一

信号处理——MATLAB小波工具箱使用简介

作者:桂. 时间:2017-02-19  21:47:27 链接:http://www.cnblogs.com/xingshansi/articles/6417638.html 声明:转载请注明出处,谢谢. 前言 本文主要介绍MATLAB小波工具箱的使用.并以一维离散信号为例,简要分析. 一.小波分解 不同于傅里叶变换,小波分解采用小波基的方式对信号进行分解,即通过基信号的平移.伸缩等变换,将信号进行分解.下图给出小波分解的一般特性: 图中可以观察到,a8对应的小波基较大,d8~d1对应的小波基

白银时代-王小波

白银时代-王小波"生活"是天籁,必须凝神静听.对不相信的事情说不在意,这就是我保全体面的方法.在剧痛中死在沙漠里,也比迷失在白银世界里要好得多.文明社会一环扣一环,和谐地运转着,错一环则动全身.你知道什么是天才的秘诀吗?那就是永远只做一件事.假如要做的事情很多,那就排出次序,依次来干.每个人的一生都拥有一些资源,比方说:寿命,智力,健康,身体,性生活.有些人准备把它消费掉,换取新奇,快乐等等.女人上街总是想猎人扛枪进山一样,但是猎取的目标有所不同.他说,我这一生都在等待,等待研究数学,

多尺度二维离散小波重构waverec2

clc,clear all,close all; load woman; [c,s]=wavedec2(X,2,'haar');%进行2尺度二维离散小波分解.分解小波函数haar %多尺度二维离散小波重构(逆变换) Y=waverec2(c,s,'haar'); figure; subplot(1,2,1),imshow(X,map),title('原始图像'); subplot(1,2,2),imshow(Y,map),title('重构图像');

3.3 哈尔小波空间W0

在3.2节我们学习了关于(3.8)定义的Vj的性质.特别的,我们可以乘以系数从一个Vj空间变换到另一个.我们这节学习V0和V1的关系. 将f1(t)∈V1投影至V0 我们考虑一个属于V1的函数f1(t),有 这个函数在图3.12中画出 图3.12 函数f1(t) 从性质2.8我们知道f1(t)属于Vj,只要j≥1.然而这个函数不属于V0,因为他的间断点在.如果我们现在想找一个在V0里面的函数f0(t)来逼近f1(t),那我们可以采用在3.2小结中学习的公式来做,我们有 因为,这个函数的支撑集为,

单尺度二维离散小波重构（逆变换）idwt2

clc,clear all,close all; load woman; %单尺度二维离散小波分解.分解小波函数haar [cA,cH,cV,cD]=dwt2(X,'haar'); %单尺度二维离散小波重构(逆变换) Y=idwt2(cA,cH,cV,cD,'haar'); figure; subplot(1,2,1),imshow(X,map),title('原始图像'); subplot(1,2,2),imshow(Y,map),title('重构图像');