哈尔尺度函数
图象处理与分析三，图象预处理图象变换，小波  补充：图像分析的数据结构  传统的图像数据结构 拓扑数据结构－把图像描述为一组元素集合和它们的关系；区域邻接曲线图是这一类型数据结构的典型,由一组节点V = {v1，v2,. . . ,vn}和一组弧 E={e1，e2,. . . ,em}组成的代数结构。     四叉树－四叉树是T-金字塔的修正，每个分层的水平上图像被分为四个象限；如果一个父节点的四个子节点有相同的值(例如,亮度),不需要分别记录他们。    图像预处理概述  图像预处理方法被分类为四个种类:象素亮度变换，几何学的变 换，使用处理象素的一个局部邻域的预处理方法, 需要有关整个 图像的知识的图像恢复。 本章的主要内容  3.1  概述和分类  3.2  傅里叶变换和性质  3.3  快速傅里叶变换  3.4  其它可分离图像变换  3.5  霍特林变换       3.1 概述和分类  例－傅立叶域的信息表达 例 例 例 分类   背景    法国数学家傅里叶生于1768年，他被世人铭记的最大贡献记载在 1807年的传记中和后来出版于1822年的"La Theorie Analitique  de la chaleur“(热分析理论)一书中。他指出任何周期函数都可以 表示为不同频率的正弦和／或余弦和的形式，每个正弦和／或余弦 乘以不同的系数(现在称这个和为傅里叶级数)。复杂函数可以由简 单的正弦和余弦之和来表示。   非周期的函数(但是这些领域是在曲线是有限的情况下)也可以用 正弦和／或余弦乘以加权函数的积分来表示。在这种情况下的公式 就是傅里叶变换，它的应用在大多数实际应用中比傅里叶级数更广 泛。   用傅里叶级数或变换表示的函数特征可以完全通过傅里叶反过程 来重建，不丢失任何信息。   在过去的一个世纪里，尤其是后50年，傅里叶的思想使整个工业 和学术界都空前繁荣。在20世纪50年代后期，数字计算的出现和快 速傅里叶变换算法的“发明”在信号处理领域产生了巨大变革。            背景 (例) 3.2  傅里叶变换和性质  3.2.2  2-D傅里叶变换  3.2.3  2-D傅里叶变换的性质   2.平移性质   4.旋转性质  离散卷积定理 －统一周期    只有当被采样函数是带限的且周期的(周期为K)，此时使采样间隔满足采样定理才可能完全恢复原来的函数。   2-D采样函数s(x,y)是由一组在X方向间隔△x，在Y方向间隔△y的2-D脉冲构成。 对一个带限函数f(x,y)(即它的傅里叶变换在某个有限区间只外为零)，如果选择如 下频域里的函数与f(x,y)和s(x,y)的乘积的傅里叶变换相乘，就有可能完全恢复 f(x,y)。 3.3.1  快速傅里叶变换算法原理-－库利、图基 只需对正变换的输入作一点小修改就可用于反变换。    如果下式成立：则称正向变换核是可分离的。    3.4.2  沃尔什变换  3.4.3  哈达玛变换  沃尔什变换与哈达玛变换-参考清华大学容观澳书  沃尔什函数的三种排序 3.4.4  离散余弦变换  3.4.5 哈尔变换--哈尔函数定义哈尔小波  哈尔基图像  3.4.6  斜变换－斜(slant)变换也称斯拉特变换。   3.5霍特林变换 －基于图象统计特性的变换  霍特林变换：  小波  1  背景  2  波和小波  3  多分辨率展开 4  一维小波变换  5  二维小波变换  小结       7.1  背景 7.1.1  图像金字塔 图像金字塔  建立图像金字塔  子带编码 (例) 波和小波     波和小波－波与小波之间的差异     时频域分析    生动例子：小波和音乐 变换的理解    变换类型      加窗傅里叶变换     盖伯变换     小波变换－一维连续小波变换      莫莱特(Morlet)小波    莫莱特(Morlet)小波     哈尔变换     离散小波变换的哈尔函数     7.2  多分辨率展开MRA –技术之一                      用V0和W0的展开图7.9(e)函数进行的结果 其中，fa(x)是f(x)使用V0尺度函数 的近似，而fd(x)为f(x)-fa(x)的差， 用W0小波和表示。这两个展开式， 如图7.12(e)和(f)所示，将f(x)用 类似高通和低通滤波器的方法分成 两部分。f(x)的低频部分在fa(x)中 得到——fa(x)给出了f(x)在每个积 分区间上的平均值——而高频细节 则在fd(x)中编码。       7.3一维小波变换 －小波变换三种形式：一般小波序列 展开、离散小波变换和连续小波变换  y=x2的哈尔小波序列展开  7.3.2  离散小波变换  一维离

哈尔尺度函数 - Read
图象处理与分析三，图象预处理图象变换，小波  补充：图像分析的数据结构  传统的图像数据结构 拓扑数据结构－把图像描述为一组元素集合和它们的关系；区域邻接曲线图是这一类型数据结构的典型,由一组节点V = {v1，v2,. . . ,vn}和一组弧 E={e1，e2,. . . ,em}组成的代数结构。     四叉树－四叉树是T-金字塔的修正，每个分层的水平上图像被分为四个象限；如果一个父节点的四个子节点有相同的值(例如,亮度),不需要分别记录他们。    图像预处理概述  图像预处理方法被分类为四个种类:象素亮度变换，几何学的变 换，使用处理象素的一个局部邻域的预处理方法, 需要有关整个 图像的知识的图像恢复。 本章的主要内容  3.1  概述和分类  3.2  傅里叶变换和性质  3.3  快速傅里叶变换  3.4  其它可分离图像变换  3.5  霍特林变换       3.1 概述和分类  例－傅立叶域的信息表达 例 例 例 分类   背景    法国数学家傅里叶生于1768年，他被世人铭记的最大贡献记载在 1807年的传记中和后来出版于1822年的"La Theorie Analitique  de la chaleur“(热分析理论)一书中。他指出任何周期函数都可以 表示为不同频率的正弦和／或余弦和的形式，每个正弦和／或余弦 乘以不同的系数(现在称这个和为傅里叶级数)。复杂函数可以由简 单的正弦和余弦之和来表示。   非周期的函数(但是这些领域是在曲线是有限的情况下)也可以用 正弦和／或余弦乘以加权函数的积分来表示。在这种情况下的公式 就是傅里叶变换，它的应用在大多数实际应用中比傅里叶级数更广 泛。   用傅里叶级数或变换表示的函数特征可以完全通过傅里叶反过程 来重建，不丢失任何信息。   在过去的一个世纪里，尤其是后50年，傅里叶的思想使整个工业 和学术界都空前繁荣。在20世纪50年代后期，数字计算的出现和快 速傅里叶变换算法的“发明”在信号处理领域产生了巨大变革。            背景 (例) 3.2  傅里叶变换和性质  3.2.2  2-D傅里叶变换  3.2.3  2-D傅里叶变换的性质   2.平移性质   4.旋转性质  离散卷积定理 －统一周期    只有当被采样函数是带限的且周期的(周期为K)，此时使采样间隔满足采样定理才可能完全恢复原来的函数。   2-D采样函数s(x,y)是由一组在X方向间隔△x，在Y方向间隔△y的2-D脉冲构成。 对一个带限函数f(x,y)(即它的傅里叶变换在某个有限区间只外为零)，如果选择如 下频域里的函数与f(x,y)和s(x,y)的乘积的傅里叶变换相乘，就有可能完全恢复 f(x,y)。 3.3.1  快速傅里叶变换算法原理-－库利、图基 只需对正变换的输入作一点小修改就可用于反变换。    如果下式成立：则称正向变换核是可分离的。    3.4.2  沃尔什变换  3.4.3  哈达玛变换  沃尔什变换与哈达玛变换-参考清华大学容观澳书  沃尔什函数的三种排序 3.4.4  离散余弦变换  3.4.5 哈尔变换--哈尔函数定义哈尔小波  哈尔基图像  3.4.6  斜变换－斜(slant)变换也称斯拉特变换。   3.5霍特林变换 －基于图象统计特性的变换  霍特林变换：  小波  1  背景  2  波和小波  3  多分辨率展开 4  一维小波变换  5  二维小波变换  小结       7.1  背景 7.1.1  图像金字塔 图像金字塔  建立图像金字塔  子带编码 (例) 波和小波     波和小波－波与小波之间的差异     时频域分析    生动例子：小波和音乐 变换的理解    变换类型      加窗傅里叶变换     盖伯变换     小波变换－一维连续小波变换      莫莱特(Morlet)小波    莫莱特(Morlet)小波     哈尔变换     离散小波变换的哈尔函数     7.2  多分辨率展开MRA –技术之一                      用V0和W0的展开图7.9(e)函数进行的结果 其中，fa(x)是f(x)使用V0尺度函数 的近似，而fd(x)为f(x)-fa(x)的差， 用W0小波和表示。这两个展开式， 如图7.12(e)和(f)所示，将f(x)用 类似高通和低通滤波器的方法分成 两部分。f(x)的低频部分在fa(x)中 得到——fa(x)给出了f(x)在每个积 分区间上的平均值——而高频细节 则在fd(x)中编码。       7.3一维小波变换 －小波变换三种形式：一般小波序列 展开、离散小波变换和连续小波变换  y=x2的哈尔小波序列展开  7.3.2  离散小波变换  一维离

在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小结中学习的公式来做，我们有
因为，这个函数的支撑集为
，所以我们能写成以下有限和的形式
由于这个内积还是十分容易计算的，我们有
我们就逐一计算
现在我们将他们都画在一张图上面
图3.13 函数f0(t)(黑)以及f1(t)(灰)
我们可以看Φ(t-k)的系数，我们不难发现，这些系数都是求平均之后再乘以根号2，我们看(3.21)的f1(t)有
其中a0=5,a1=3，a3=-1,a4=5,a5=7，然后
b0=4√2，b1=2√2,b2=5√2。我们会发现bk遵循以下规律
我们下面对这个例子进行推广，推广至所有的f1(t)∈V1
性质3.9(从V1投影到V0)如果f1(t)∈V1，表示为
同样f1(t)的投影f0(t)=Pf1,0(t)可以写成以下的形式
其中
(3.23)
对所有k∈Z成立
证明：让
现在我们知道
，并且函数的Φ(t-k)的值为1。所以积分就变成了
由于断点在(1/2)Z 处，所以我们把函数分成两部分
，其中
和
。其中，
]
我们下面来确定m的值。
我们可以采用在3.2节的练习3.11来确定。Φ(2t-m)的支撑集为
。我们将左端点m/2设定为左端，也就是说，m/2=k或者m=2k,因此Φ(2t-2k)=Φ1,2k(t)。最左端的重叠部分为[k,k+1]。
从而我们有
因此在区间[k,k+1]
其中，支撑集
，
因此内积变作
在其支撑集上Φ(2t-2k)=1，Φ(2t-(2k+1))=1.内积的结果为
证毕。
剩余函数g(t)=f1(t)-f0(t)
我们将f0(t)视作f1(t)∈V1的逼近，如果我们现在建立一个函数g(t)使得f0(t)+g0(t)=f1(t)，那么f1(t)会是怎么样的呢？
我们能够非常容易的计算出g0(t)=f1(t)-f0(t)，从图3.13我们可以知道g(t)是一个在1/2 ,1 ,3/2 ,2,5/2处有间断点的常分段函数。g(t)是V1的元素，能写成Φ1,k(t)的线性组合。
我们将g0(t)写出来
g0(t)在图3.14中画出
图3.14 函数g0(t)
小波函数Ψ(t)的定义
我们对g0(t)的进行更细一步的分析。如果g0(t)在[0.1)上。那么这个函数g(t)乘以√2，会形成函数Ψ(t)
图3.15 函数 Ψ(t)
我们看g0(t)在[1,2)区间内，可以看做Ψ(t)右移一个单位并乘以3√2。因此在[1.2)处，g0(t)=3√2Ψ(t-1)。同样的，我们也可以写出在区间[2,3)。g0(t)=-Ψ(t-2)。总体来说
我们应该注意到，他们的系数都是有一定的关系的。我们回忆一下(3.21)中f1(t)的系数
。如果我们让
则
所以，不难看出ck遵循
用Ψ(t)书写g0(t)
同样的，我们对其进行推广到所有的g0(t)
性质3.10(剩余函数的公式)假设f1(t)是V1的一个函数，有
f0(t)是在v0空间对f1(t)的逼近。并且有g(t)=f1(t)-f0(t)，g0(t)为
同时，Ψ(t)为(3.25)，以及系数为
证明：证明留作练习3.22
空间W0的定义
回忆我们对于V0的定义，是Φ(t)以及它的整数平移的函数的线性组合。仿照这个定义，我们也规定了Ψ(t)以及它的整数平移的函数的线性组合张成以下空间
定义3.4(剩余空间W0)对于在(3.25)定义的Ψ(t)，我们定义以下空间
我们称这个空间为哈尔小波函数Ψ(t)生成的哈尔小波空间。
我们看我这个小波的这个词是怎么来的，我们从f0(t)来拟合f1(t)，会产生些“小”的波的误差如Ψ(t)以及它的整数平移。
通过定义
张成空间W0，而且在练习3.23中你会知道，他们是线性无关的函数。同样以下的式子很容易得出(练习3.25)
从而有以下的性质
性质3.11(W0的一组正交基)对于(3.28)定义出来的空间W0，以及(3.25)定义出来的Ψ(t)来说。
是它的一组正交基。
证明：这个问题留作练习3.23以及3.25
尺度函数以及膨胀方程
从(3.23)中你可以看见向量
在f1(t)投影到V0空间中扮演了重要的作用，因此，我们能通过(3.29)重写(3.23)。我们对于k∈Z
其中
，我们能通过h来让V1里面的基函数表示Φ(t)，有
同样，我们可以使用向量
将f1(t)投影进W0。正如我们在(3.25)所见的
因此，我们能将Ψ(t)用V1的基函数来表示，滤波系数g0以及g1可以为其系数的线性组合。
等式(3.30)和(3.32)在小波理论中担当重要地位。函数Φ(t)生成了能让我们拟合在L2(R)里面函数的空间Vj。等式(3.30)给予我们从V0和V1空间联系的方式。
同样，Ψ(t)生成的空间W0给我们一种度量从V0拟合V1的失误程度。等式(3.32)也说明了V1里面的函数如何表示Ψ(t)。
我们已经将Ψ(t)称作哈尔小波函数。现在我们给Φ(t)和等式(3.30)和(3.32)来命名。
定义3.5(哈尔尺度函数以及膨胀方程)我们称等式
叫膨胀方程。哈尔函数Φ(t)生成这些用于生成膨胀方程，我们称作尺度函数(作者注：Φ(t)和Ψ(t)有时候我们分别称为父小波和母小波)
性质3.9的另一种证法
我们将Φ(t)表示为Φ1,0(t)和Φ1,1(t)，给予我们另一种证明性质3.9的方法
证明：让
，假设我们将f1(t)投影入V0，使用(3.16)，我们有
而且，
，因此
因此
因为
我们第一个内积可以写为
因为Φ1,k(t)和Φ1,2m(t)是正交的，所以唯一一个非正交的内积就是k=2m 的时候，我们有
同样的，我们有
因此得到(3.23)
在练习3.27你会用(3.34)从而形成另一种方法来证明性质3.10
直和以及正交直和
性质3.12(V0、W0空间的正交性)如果f(t)∈V0，g(t)∈W0。则有=0
证明：因为f(t)∈V0，我们可以将它写作Φ(t)的线性组合
同样，
内积就是
在练习3.25(c)你会要求证明
，从而完成整个的证明。
在后面的章节，我们会对正交的集合非常感兴趣，因此，我们先在这里给出定义。
定义3.6(正交空间)假设V和W是L2(R)的两个子集。如果f(t)∈V，g(t)∈W且=0，则称它们是正交的。写作V⊥W。
在性质3.12中，我们给出了一组符合定义3.6的例子(V0⊥W0)。
我们在性质3.10里面说到f1(t)∈V1，f1(t)=g0(t)+V0(t)，其中g0∈W0，f0(t)∈V0，同时f0(t)是f1(t)在V0的投影，g0(t)是f1(t)在W0的投影。(见练习3.26)我们可以通过W0和V0里面的函数相加构造出V1里面的函数。我们给出下面这个定义
定义3.7(直和以及正交直和)假设V和W都是L2(R)的子空间，我们定义V和W的直和为
如果V⊥W，我们称X为正交直和，并记作
注意到，如果X=W+V，由因为0∈W，因而又有
，所以有
。相似的，我们也可以证明W是X的子空间。
从V0和W0到V1
我们已经知道了V0⊥W0，那么他们的正交直和是什么呢？下面的这个性质给出了答案
性质3.13(V0和W0的正交直和)如果V0和V1分别定义于(3.8)和(3.28)，我们有
证明:从性质3.10我们知道任意f1(t)∈V1能表示为f0(t)∈V0和g0∈W0.我们也知道f0(t)和g0(t)都是V1中的元素，因此原命题等价于f1(t)∈V1与所有的h(t)∈V0正交，则f1(t)一定属于W0(译者注：下面的讨论是考虑V1在V0的投影为0的特殊情况)。假设h(t)是任意一个V0里面的函数，然后
又因为h(t)是任意的，所以必须有f0(t)=0因此f1(t)∈W0，证毕。

另外参见俄罗斯写的http://www.codeproject.com/Articles/22243/Real-Time-Object-Tracker-in-C
Haar小波在图像处理和数字水印等方面应用较多，这里简单的介绍一下哈尔小波的基本原理以及其实现情况。
一、Haar小波的基本原理
数学理论方面的东西我也不是很熟悉，这边主要用简单的例子来介绍下Haar小波的使用情况。
例如：有a=[8,7,6,9]四个数，并使用b[4]数组来保存结果.
则一级Haar小波变换的结果为:
b[0]=(a[0]+a[1])/2,                       b[2]=(a[0]-a[1])/2
b[1]=(a[2]+a[3])/2,                       b[3]=(a[2]-a[3])/2
即依次从数组中取两个数字，计算它们的和以及差，并将和一半和差的一半依次保存在数组的前半部分和后半部分。
例如：有a[8]，要进行一维Haar小波变换，结果保存在b[8]中
则一级Haar小波变换的结果为:
b[0]=(a[0]+a[1])/2,                        b[4]=(a[0]-a[1])/2
b[1]=(a[2]+a[3])/2,                        b[5]=(a[2]-a[3])/2
b[2]=(a[4]+a[5])/2,                        b[6]=(a[4-a[5]])/2
b[3]=(a[6]+a[7])/2,                        b[7]=(a[6]-a[7])/2
如果需要进行二级Haar小波变换的时候，只需要对b[0]-b[3]进行Haar小波变换.
对于二维的矩阵来讲，每一级Haar小波变换需要先后进行水平方向和竖直方向上的两次一维小波变换,行和列的先后次序对结果不影响。
二、Haar小波的实现
使用opencv来读取图片及像素，对图像的第一个8*8的矩阵做了一级小波变换
#include 
#include 
#include 
using namespace std;
int main(){
IplImage* srcImg;
double  imgData[8][8];
int i,j;
srcImg=cvLoadImage("lena.bmp",0);
cout<
for( i=0;i<8;i++)
{
for( j=0;j<8;j++)
{
imgData[i][j]=cvGetReal2D(srcImg,i+256,j+16);
cout<
}
cout<
}
double tempData[8];
//行小波分解
for( i=0;i<8;i++)
{
for( j=0;j<4;j++)
{
double temp1=imgData[i][2*j];
double temp2=imgData[i][2*j+1];
tempData[j]=(temp1+temp2)/2;
tempData[j+4]=(temp1-temp2)/2;
}
for( j=0;j<8;j++)
imgData[i][j]=tempData[j];
}
//列小波分解
for( i=0;i<8;i++)
{
for( j=0;j<4;j++)
{
double temp1=imgData[2*j][i];
double temp2=imgData[2*j+1][i];
tempData[j]=(temp1+temp2)/2;
tempData[j+4]=(temp1-temp2)/2;
}
for( j=0;j<8;j++)
imgData[j][i]=tempData[j];
}
cout<
for( i=0;i<8;i++)
{
for( j=0;j<8;j++)
{
cout<
}
cout<
}
//列小波逆分解
for( i=0;i<8;i++)
{
for( j=0;j<4;j++)
{
double temp1=imgData[j][i];
double temp2=imgData[j+4][i];
tempData[2*j]=temp1+temp2;
tempData[2*j+1]=temp1-temp2;
}
for( j=0;j<8;j++)
{
imgData[j][i]=tempData[j];
}
}
//行小波逆分解
for( i=0;i<8;i++)
{
for( j=0;j<4;j++)
{
double temp1=imgData[i][j];
double temp2=imgData[i][j+4];
tempData[2*j]=temp1+temp2;
tempData[2*j+1]=temp1-temp2;
}
for( j=0;j<2*4;j++)
{
imgData[i][j]=tempData[j];
}
}
cout<
for( i=0;i<8;i++)
{
for( j=0;j<8;j++)
{
cout<
}
cout<
}
return 0;
}
===========================================================================
///  小波变换Mat WDT( const Mat &_src, const string _wname, const int _level )const{int reValue = THID_ERR_NONE;Mat src = Mat_(_src);Mat dst = Mat::zeros( src.rows, src.cols, src.type() );int N = src.rows;int D = src.cols;/// 高通低通滤波器Mat lowFilter;Mat highFilter;wavelet( _wname, lowFilter, highFilter );/// 小波变换int t=1;int row = N;int col = D;while( t<=_level ){///先进行行小波变换for( int i=0; i(0,j) = src.at(i,j);}oneRow = waveletDecompose( oneRow, lowFilter, highFilter );/// 将src这一行置为oneRow中的数据for ( int j=0; j(i,j) = oneRow.at(0,j);}}#if 0//normalize( dst, dst, 0, 255, NORM_MINMAX );IplImage dstImg1 = IplImage(dst);cvSaveImage( "dst.jpg", &dstImg1 );#endif/// 小波列变换for ( int j=0; j(i,0) = dst.at(i,j);}oneCol = ( waveletDecompose( oneCol.t(), lowFilter, highFilter ) ).t();for ( int i=0; i(i,j) = oneCol.at(i,0);}}#if 0//normalize( dst, dst, 0, 255, NORM_MINMAX );IplImage dstImg2 = IplImage(dst);cvSaveImage( "dst.jpg", &dstImg2 );#endif/// 更新row /= 2;col /=2;t++;src = dst;}return dst;}///  小波逆变换Mat IWDT( const Mat &_src, const string _wname, const int _level )const{int reValue = THID_ERR_NONE;Mat src = Mat_(_src);Mat dst = Mat::zeros( src.rows, src.cols, src.type() );int N = src.rows;int D = src.cols;/// 高通低通滤波器Mat lowFilter;Mat highFilter;wavelet( _wname, lowFilter, highFilter );/// 小波变换int t=1;int row = N/std::pow( 2., _level-1);int col = D/std::pow(2., _level-1);while ( row<=N && col<=D ){/// 小波列逆变换for ( int j=0; j(i,0) = src.at(i,j);}oneCol = ( waveletReconstruct( oneCol.t(), lowFilter, highFilter ) ).t();for ( int i=0; i(i,j) = oneCol.at(i,0);}}#if 0//normalize( dst, dst, 0, 255, NORM_MINMAX );IplImage dstImg2 = IplImage(dst);cvSaveImage( "dst.jpg", &dstImg2 );#endif///行小波逆变换for( int i=0; i(0,j) = dst.at(i,j);}oneRow = waveletReconstruct( oneRow, lowFilter, highFilter );/// 将src这一行置为oneRow中的数据for ( int j=0; j(i,j) = oneRow.at(0,j);}}#if 0//normalize( dst, dst, 0, 255, NORM_MINMAX );IplImage dstImg1 = IplImage(dst);cvSaveImage( "dst.jpg", &dstImg1 );#endifrow *= 2;col *= 2;src = dst;}return dst;}/// 调用函数/// 生成不同类型的小波，现在只有haar，sym2void wavelet( const string _wname, Mat &_lowFilter, Mat &_highFilter )const{if ( _wname=="haar" || _wname=="db1" ){int N = 2;_lowFilter = Mat::zeros( 1, N, CV_32F );_highFilter = Mat::zeros( 1, N, CV_32F );_lowFilter.at(0, 0) = 1/sqrtf(N);_lowFilter.at(0, 1) = 1/sqrtf(N);_highFilter.at(0, 0) = -1/sqrtf(N);_highFilter.at(0, 1) = 1/sqrtf(N);}if ( _wname =="sym2" ){int N = 4;float h[] = {-0.483, 0.836, -0.224, -0.129 };float l[] = {-0.129, 0.224,    0.837, 0.483 };_lowFilter = Mat::zeros( 1, N, CV_32F );_highFilter = Mat::zeros( 1, N, CV_32F );for ( int i=0; i(0, i) = l[i];_highFilter.at(0, i) = h[i];}}}/// 小波分解Mat waveletDecompose( const Mat &_src, const Mat &_lowFilter, const Mat &_highFilter )const{assert( _src.rows==1 && _lowFilter.rows==1 && _highFilter.rows==1 );assert( _src.cols>=_lowFilter.cols && _src.cols>=_highFilter.cols );Mat &src = Mat_(_src);int D = src.cols;Mat &lowFilter = Mat_(_lowFilter);Mat &highFilter = Mat_(_highFilter);/// 频域滤波，或时域卷积；ifft( fft(x) * fft(filter)) = cov(x,filter)Mat dst1 = Mat::zeros( 1, D, src.type() );Mat dst2 = Mat::zeros( 1, D, src.type()  );filter2D( src, dst1, -1, lowFilter );filter2D( src, dst2, -1, highFilter );/// 下采样Mat downDst1 = Mat::zeros( 1, D/2, src.type() );Mat downDst2 = Mat::zeros( 1, D/2, src.type() );resize( dst1, downDst1, downDst1.size() );resize( dst2, downDst2, downDst2.size() );/// 数据拼接for ( int i=0; i(0, i) = downDst1.at( 0, i );src.at(0, i+D/2) = downDst2.at( 0, i );}return src;}/// 小波重建Mat waveletReconstruct( const Mat &_src, const Mat &_lowFilter, const Mat &_highFilter )const{assert( _src.rows==1 && _lowFilter.rows==1 && _highFilter.rows==1 );assert( _src.cols>=_lowFilter.cols && _src.cols>=_highFilter.cols );Mat &src = Mat_(_src);int D = src.cols;Mat &lowFilter = Mat_(_lowFilter);Mat &highFilter = Mat_(_highFilter);/// 插值;Mat Up1 = Mat::zeros( 1, D, src.type() );Mat Up2 = Mat::zeros( 1, D, src.type() );/// 插值为0//for ( int i=0, cnt=1; i( 0, cnt ) = src.at( 0, i );     ///< 前一半//    Up2.at( 0, cnt ) = src.at( 0, i+D/2 ); ///< 后一半//}/// 线性插值Mat roi1( src, Rect(0, 0, D/2, 1) );Mat roi2( src, Rect(D/2, 0, D/2, 1) );resize( roi1, Up1, Up1.size(), 0, 0, INTER_CUBIC );resize( roi2, Up2, Up2.size(), 0, 0, INTER_CUBIC );/// 前一半低通，后一半高通Mat dst1 = Mat::zeros( 1, D, src.type() );Mat dst2= Mat::zeros( 1, D, src.type() );filter2D( Up1, dst1, -1, lowFilter );filter2D( Up2, dst2, -1, highFilter );/// 结果相加dst1 = dst1 + dst2;return dst1;}
===========================================================================
*****************************************************************************************************************
其他代码实现
*****************************************************************************************************************
// 哈尔小波.cpp : 定义控制台应用程序的入口点。
//
#include "stdafx.h"
#include 
#include "cv.h"
#include "highgui.h"
#include 
// 二维离散小波变换(单通道浮点图像)
void DWT(IplImage *pImage, int nLayer)
{
// 执行条件
if (pImage)
{
if (pImage->nChannels == 1 &&
pImage->depth == IPL_DEPTH_32F &&
((pImage->width >> nLayer) << nLayer) == pImage->width &&
((pImage->height >> nLayer) << nLayer) == pImage->height)
{
int     i, x, y, n;
float   fValue   = 0;
float   fRadius  = sqrt(2.0f);
int     nWidth   = pImage->width;
int     nHeight  = pImage->height;
int     nHalfW   = nWidth / 2;
int     nHalfH   = nHeight / 2;
float **pData    = new float*[pImage->height];
float  *pRow     = new float[pImage->width];
float  *pColumn  = new float[pImage->height];
for (i = 0; i < pImage->height; i++)
{
pData[i] = (float*) (pImage->imageData + pImage->widthStep * i);
}
// 多层小波变换
for (n = 0; n < nLayer; n++, nWidth /= 2, nHeight /= 2, nHalfW /= 2, nHalfH /= 2)
{
// 水平变换
for (y = 0; y < nHeight; y++)
{
// 奇偶分离
memcpy(pRow, pData[y], sizeof(float) * nWidth);
for (i = 0; i < nHalfW; i++)
{
x = i * 2;
pData[y][i] = pRow[x];
pData[y][nHalfW + i] = pRow[x + 1];
}
// 提升小波变换
for (i = 0; i < nHalfW - 1; i++)
{
fValue = (pData[y][i] + pData[y][i + 1]) / 2;
pData[y][nHalfW + i] -= fValue;
}
fValue = (pData[y][nHalfW - 1] + pData[y][nHalfW - 2]) / 2;
pData[y][nWidth - 1] -= fValue;
fValue = (pData[y][nHalfW] + pData[y][nHalfW + 1]) / 4;
pData[y][0] += fValue;
for (i = 1; i < nHalfW; i++)
{
fValue = (pData[y][nHalfW + i] + pData[y][nHalfW + i - 1]) / 4;
pData[y][i] += fValue;
}
// 频带系数
for (i = 0; i < nHalfW; i++)
{
pData[y][i] *= fRadius;
pData[y][nHalfW + i] /= fRadius;
}
}
// 垂直变换
for (x = 0; x < nWidth; x++)
{
// 奇偶分离
for (i = 0; i < nHalfH; i++)
{
y = i * 2;
pColumn[i] = pData[y][x];
pColumn[nHalfH + i] = pData[y + 1][x];
}
for (i = 0; i < nHeight; i++)
{
pData[i][x] = pColumn[i];
}
// 提升小波变换
for (i = 0; i < nHalfH - 1; i++)
{
fValue = (pData[i][x] + pData[i + 1][x]) / 2;
pData[nHalfH + i][x] -= fValue;
}
fValue = (pData[nHalfH - 1][x] + pData[nHalfH - 2][x]) / 2;
pData[nHeight - 1][x] -= fValue;
fValue = (pData[nHalfH][x] + pData[nHalfH + 1][x]) / 4;
pData[0][x] += fValue;
for (i = 1; i < nHalfH; i++)
{
fValue = (pData[nHalfH + i][x] + pData[nHalfH + i - 1][x]) / 4;
pData[i][x] += fValue;
}
// 频带系数
for (i = 0; i < nHalfH; i++)
{
pData[i][x] *= fRadius;
pData[nHalfH + i][x] /= fRadius;
}
}
}
delete[] pData;
delete[] pRow;
delete[] pColumn;
}
}
}
// 二维离散小波恢复(单通道浮点图像)
void IDWT(IplImage *pImage, int nLayer)
{
// 执行条件
if (pImage)
{
if (pImage->nChannels == 1 &&
pImage->depth == IPL_DEPTH_32F &&
((pImage->width >> nLayer) << nLayer) == pImage->width &&
((pImage->height >> nLayer) << nLayer) == pImage->height)
{
int     i, x, y, n;
float   fValue   = 0;
float   fRadius  = sqrt(2.0f);
int     nWidth   = pImage->width >> (nLayer - 1);
int     nHeight  = pImage->height >> (nLayer - 1);
int     nHalfW   = nWidth / 2;
int     nHalfH   = nHeight / 2;
float **pData    = new float*[pImage->height];
float  *pRow     = new float[pImage->width];
float  *pColumn  = new float[pImage->height];
for (i = 0; i < pImage->height; i++)
{
pData[i] = (float*) (pImage->imageData + pImage->widthStep * i);
}
// 多层小波恢复
for (n = 0; n < nLayer; n++, nWidth *= 2, nHeight *= 2, nHalfW *= 2, nHalfH *= 2)
{
// 垂直恢复
for (x = 0; x < nWidth; x++)
{
// 频带系数
for (i = 0; i < nHalfH; i++)
{
pData[i][x] /= fRadius;
pData[nHalfH + i][x] *= fRadius;
}
// 提升小波恢复
fValue = (pData[nHalfH][x] + pData[nHalfH + 1][x]) / 4;
pData[0][x] -= fValue;
for (i = 1; i < nHalfH; i++)
{
fValue = (pData[nHalfH + i][x] + pData[nHalfH + i - 1][x]) / 4;
pData[i][x] -= fValue;
}
for (i = 0; i < nHalfH - 1; i++)
{
fValue = (pData[i][x] + pData[i + 1][x]) / 2;
pData[nHalfH + i][x] += fValue;
}
fValue = (pData[nHalfH - 1][x] + pData[nHalfH - 2][x]) / 2;
pData[nHeight - 1][x] += fValue;
// 奇偶合并
for (i = 0; i < nHalfH; i++)
{
y = i * 2;
pColumn[y] = pData[i][x];
pColumn[y + 1] = pData[nHalfH + i][x];
}
for (i = 0; i < nHeight; i++)
{
pData[i][x] = pColumn[i];
}
}
// 水平恢复
for (y = 0; y < nHeight; y++)
{
// 频带系数
for (i = 0; i < nHalfW; i++)
{
pData[y][i] /= fRadius;
pData[y][nHalfW + i] *= fRadius;
}
// 提升小波恢复
fValue = (pData[y][nHalfW] + pData[y][nHalfW + 1]) / 4;
pData[y][0] -= fValue;
for (i = 1; i < nHalfW; i++)
{
fValue = (pData[y][nHalfW + i] + pData[y][nHalfW + i - 1]) / 4;
pData[y][i] -= fValue;
}
for (i = 0; i < nHalfW - 1; i++)
{
fValue = (pData[y][i] + pData[y][i + 1]) / 2;
pData[y][nHalfW + i] += fValue;
}
fValue = (pData[y][nHalfW - 1] + pData[y][nHalfW - 2]) / 2;
pData[y][nWidth - 1] += fValue;
// 奇偶合并
for (i = 0; i < nHalfW; i++)
{
x = i * 2;
pRow[x] = pData[y][i];
pRow[x + 1] = pData[y][nHalfW + i];
}
memcpy(pData[y], pRow, sizeof(float) * nWidth);
}
}
delete[] pData;
delete[] pRow;
delete[] pColumn;
}
}
}
int _tmain(int argc, _TCHAR* argv[])
{
// 小波变换层数
int nLayer = 1;
// 输入彩色图像
IplImage *pSrc = cvLoadImage("lena.bmp", CV_LOAD_IMAGE_COLOR);
// 计算小波图象大小
CvSize size = cvGetSize(pSrc);
if ((pSrc->width >> nLayer) << nLayer != pSrc->width)
{
size.width = ((pSrc->width >> nLayer) + 1) << nLayer;
}
if ((pSrc->height >> nLayer) << nLayer != pSrc->height)
{
size.height = ((pSrc->height >> nLayer) + 1) << nLayer;
}
// 创建小波图象
IplImage *pWavelet = cvCreateImage(size, IPL_DEPTH_32F, pSrc->nChannels);
if (pWavelet)
{
// 小波图象赋值
cvSetImageROI(pWavelet, cvRect(0, 0, pSrc->width, pSrc->height));
cvConvertScale(pSrc, pWavelet, 1, -128);
cvResetImageROI(pWavelet);
// 彩色图像小波变换
IplImage *pImage = cvCreateImage(cvGetSize(pWavelet), IPL_DEPTH_32F, 1);
if (pImage)
{
for (int i = 1; i <= pWavelet->nChannels; i++)
{
cvSetImageCOI(pWavelet, i);
cvCopy(pWavelet, pImage, NULL);
// 二维离散小波变换
DWT(pImage, nLayer);
// 二维离散小波恢复
// IDWT(pImage, nLayer);
cvCopy(pImage, pWavelet, NULL);
}
cvSetImageCOI(pWavelet, 0);
cvReleaseImage(&pImage);
}
// 小波变换图象
cvSetImageROI(pWavelet, cvRect(0, 0, pSrc->width, pSrc->height));
cvConvertScale(pWavelet, pSrc, 1, 128);
cvNamedWindow("pWavelet",CV_WINDOW_AUTOSIZE);
cvShowImage("pWavelet",pWavelet);
cvNamedWindow("pSrc",CV_WINDOW_AUTOSIZE);
cvShowImage("pSrc",pSrc);
cvWaitKey(0);
cvResetImageROI(pWavelet); // 本行代码有点多余，但有利用养成良好的编程习惯
cvReleaseImage(&pWavelet);
}
// 显示图像pSrc
// ...
cvReleaseImage(&pSrc);
return 0;
}
************************
-------------------------------------Codeproject Haar -------------------------------------------------------------------------------------------------
void Haar::transrows(char** dest, char** sour, unsigned int w, unsigned int h) const
{
unsigned int w2 = w / 2;
__m64 m00FF;
m00FF.m64_u64 = 0x00FF00FF00FF00FF;
for (unsigned int y = 0; y < h; y++) {
__m64 *mlo = (__m64 *) & dest[y][0];
__m64 *mhi = (__m64 *) & dest[y][w2];
__m64 *msour = (__m64 *) & sour[y][0];
for (unsigned int k = 0; k < w2 / 8; k++) {   //k
__m64 even = _mm_packs_pu16(_mm_and_si64(*msour, m00FF), _mm_and_si64(*(msour + 1), m00FF));       //even coeffs
__m64 odd = _mm_packs_pu16(_mm_srli_pi16(*msour, 8), _mm_srli_pi16(*(msour + 1), 8));              //odd coeffs
addsub(even, odd, mlo++, mhi++);
msour += 2;
}
if (w2 % 8) {
for (unsigned int k = w2 - (w2 % 8); k < w2; k++) {
dest[y][k] = char(((int)sour[y][2*k] + (int)sour[y][2*k+1]) / 2);
dest[y][k+w2] = char(((int)sour[y][2*k] - (int)sour[y][2*k+1]) / 2);
}
}
}
_mm_empty();
}
void Haar::transcols(char** dest, char** sour, unsigned int w, unsigned int h) const
{
unsigned int h2 = h / 2;
for (unsigned int k = 0; k < h2; k++) {
__m64 *mlo = (__m64 *) & dest[k][0];
__m64 *mhi = (__m64 *) & dest[k+h2][0];
__m64 *even = (__m64 *) & sour[2*k][0];
__m64 *odd = (__m64 *) & sour[2*k+1][0];
for (unsigned int x = 0; x < w / 8; x++) {
addsub(*even, *odd, mlo, mhi);
even++;
odd++;
mlo++;
mhi++;
}
}
_mm_empty();
//odd remainder
for (unsigned int x = w - (w % 8); x < w; x++) {
for (unsigned int k = 0; k < h2; k++) {
dest[k][x] = char(((int)sour[2*k][x] + (int)sour[2*k+1][x]) / 2);
dest[k+h2][x] = char(((int)sour[2*k][x] - (int)sour[2*k+1][x]) / 2);
}
}
}
****************************************************************************************************************
1.哈尔基函数
最简单的基函数是哈尔基函数(Haar basis
function)。哈尔基函数在1909年提出，它是由一组分段常值函数组成的函数集。这个函数集定义在半开区间[0,1)上，每一个分段常值函数的数
值在一个小范围里是1，其他地方为0，现以图像为例并使用线性代数中的矢量空间来说明哈尔基函数。
如果一幅图像仅由2^0=1个像素组成，这幅图像在整个[0,1)
区间中就是一个常值函数。
用q_00(x)表示这个常值函数，用V0表示由这个常值函数生成的矢量空间，
即V0：q_00(x)=1(0<=x<1)或0(其它),它是构成矢量空间V0的基。
如果一幅图像由2^1=2个像素组成，这幅图像在[0,1)
区间中有两个等间隔的子区间：[0,1/2) 和[1/2,1) ，每一个区间中各有1个常值函数，分别用q_10(x)
q_11(x)表示。用V1表示由2个子区间中的常值函数生成的矢量空间，即V1: q_10(x)  =1(0<=x<1/2)或0(其它);
q_11(x)=1(1/2<=x<1)或0(其它).这2个常值函数就是构成矢量空间V1的基。
如果一幅图像由2^2=
4个像素组成，这幅图像在[0,1)区间中被分成4个等间隔的子区间：[0,1/4),[1/4,1/2),[1/2,3/4)和[3/4,1)，它们的
常值函数分别用q_20(x) q_21(x) q_22(x)
q_23(x)表示，用V2表示由4个子区间中的常值函数生成的矢量空间，即V2:q_20(x)=1(0<=x<1/4)或0(其它);
q_21(x)=1(1/4<=x<1/2)或0(其它); q_22(x) =1(1/2<=x<3/4)或0(其它);
q_23(x)=1(3/4<=x<1)或0(其它).这4个常值函数就是构成矢量空间V2的基。
我们可以按照这种方法继续定义基函数和由它生成的矢量空间。总而言之，为了表示矢量空间中的矢量，每一个矢量空间V_i都需要定义一个基(basis)。
为生成矢量空间V_i而定义的基函数也叫做尺度函数,这种函数通常用符号q_ij(x)表示。哈尔基函数定义为q(x)=1(0<=x<1)
或0(其它),哈尔尺度函数q_ij(x)定义为q_ij(x)=q(2^i·x-j),j=0,1,…,
2^i-1其中,i为尺度因子，改变i使函数图形缩小或者放大,j为平移参数，改变j使函数沿x轴方向平移。
空间矢量V_i定义为V_i=span{
q_ij(x)},j=0,1,…,
2^i-1.由于定义了基和矢量空间，我们就可以把由2^i个像素组成的一维图像看成为矢量空间V_i中的矢量。由于这些矢量都是在单位区间[0,1)上
定义的函数，所以在V_i矢量空间中的每一个矢量也被包含在V_i+1矢量空间中，即V0∈V1∈...∈V_i∈V_i+1.2.哈尔小波函数
小波函数通常用Ψ_ij(x)表示。与哈尔基函数相对应的小波称为基本哈尔小波函数,并由下式定义:
Ψ(x)=1(0<=x<1/2)或
-1(1/2<=x<1)或0(其它);
哈尔小波尺度函数Ψ_ij(x)定义为
Ψ_ij(x)= Ψ(2^i·x-j),j=0,1,…,
2^i-1.用小波函数构成的矢量空间Wi表示，Wi=span{Ψ_ij(x)}, j=0,1,…, 2^i-1
根据哈尔小波函数的定义，可以写出生成W0，W1和W2等矢量空间的小波函数。
生成矢量空间W0的哈尔小波：Ψ_00(x)=1(0<=x<1/2)或-1(1/2<=x<1)或0(其它)
生成矢量空间W1的哈尔小波：Ψ_10(x)=1(0<=x<1/4)或-1(1/4<=x<1/2)或0(其它);
Ψ_11(x)=1(1/2<=x<3/4)或-1(3/4<=x<1/2)或0(其它).
生成矢量空间W2的哈尔小波：Ψ_00(x)=1(0<=x<1/2)或-1(1/2<=x<1)或0(其它)
Ψ_20(x)=1(0<=x<1/8)或-1(1/8<=x<2/8)或0(其它);
Ψ_21(x)=1(2/8<=x<3/8)或-1(3/8<=x<4/8)或0(其它);
Ψ_22(x)=1(4/8<=x<5/8)或-1(5/8<=x<6/8)或0(其它);
Ψ_00(x)=1(6/8<=x<7/8)或-1(7/8<=x<1)或0(其它).
哈尔基和哈尔小波分别使用下面两个式子进行规范化
q_ij(x)=2^(i/2)·q(2^i·x-j),
Ψ_ij(x)=
2^(i/2)·Ψ(2^i·x-j),
3.哈尔基的结构
使用哈尔基函数
q_ij(x)和哈尔小波函数Ψ_ij(x)生成的矢量空间Vi和Wi具有下面的性质，V_i+1=V_i⊕W_i。这就是说，在矢量空间V_i+1中,
矢量空间W_i中的小波可用来表示一个函数在矢量空间V_i
中不能表示的部分。因此，小波可定义为生成矢量空间W_i的一组线性无关的函数Ψ_ij(x)的集合。
4.哈尔小波变换
小波变换的基本思想是用一组小波函数或者基函数表示一个函数或者信号，例如图像信号。为了理解什么是小波变换，下面用一个具体的例子来说明小波变换的过程。假设有一幅分辨率只有4个像素的一维图像，对应的像素值分别为[9
7 3 5]
计算它的哈尔小波变换系数:
(1).求均值(averaging)。计算相邻像素对的平均值，得到一幅分辨率比较低的新图像，它的像素数目变成了2个，即新的图像的分辨率是原来的1/2，相应的像素值为：[8
4]
(2).
求差值(diqqerencing)。很明显，用2个像素表示这幅图像时，图像的信息已经部分丢失。为了能够从由2个像素组成的图像重构出由4个像素组成
的原始图像，就需要存储一些图像的细节系数,以便在重构时找回丢失的信息。方法是把像素对的第一个像素值减去这个像素对的平均值，或者使用这个像素对的差
值除以2。在这个例子中，第一个细节系数是(9-8)=1，因为计算得到的平均值是8,它比9小1而比7大1，存储这个细节系数就可以恢复原始图像的前两
个像素值。使用同样的方法，第二个细节系数是(3-4)=-1，存储这个细节系数就可以恢复后2个像素值。因此，原始图像就可以用下面的两个平均值和两个
细节系数表示[8 4 1 -1]
(3).重复第1，2步，把由第一步分解得到的图像进一步分解成分辨率更低的图像[6 2 1 -1].
由此可见，通过上述分解就把由4像素组成的一幅图像用一个平均像素值和三个细节系数表示，这个过程就叫做哈尔小波变换,也称哈尔小波分解(Haar wavelet
decomposition)。这个概念可以推广到使用其他小波基的变换。在这个例子中我们可以看到，①变换过程中没有丢失信息，因为能够从所记录的数据
中重构出原始图像。②对这个给定的变换，我们可以从所记录的数据中重构出各种分辨率的图像。例如，在分辨率为1的图像基础上重构出分辨率为2的图像，在分
辨率为2的图像基础上重构出分辨率为4的图像。③通过变换之后产生的细节系数的幅度值比较小，这就为图像压缩提供了一种途径，例如去掉一些微不足道的细节
系数而不影响对重构图像的理解。
求均值和差值的过程实际上就是一维小波变换的过程，现在用数学方法重新描述小波变换的过程。
(1)
用V2中的哈尔基表示
图像I(x)=[9 7 3 5]有2^j =2^2=4像素，因此可以用生成矢量空间V2中的基函数的线性组合表示，
I(x) =9
q_20(x)+
7q_21(x)+3q_22(x)+5q_23(x)
(2)用V1和W1中的函数表示
根据V2=V1⊕W1,I(x)可表示成
I(x)=8q_10(x)+4q_11(x)+Ψ_10(x)-Ψ_11(x)
(3)用V0,W0和W1中的函数表示
V2=V0⊕W0⊕W1,I(x)可表示成
I(x)
=6 q_00(x)+2Ψ_00(x)+Ψ_10(x)-Ψ_11(x)
5.规范化算法
规范化的小波变换与非规范化的小波变换相比，唯一的差别是在规范化的小波变换中需要乘一个规范化的系数。下面用一个例子说明。对函数f(x)=
[2,5,8,9,7,4, -1,
1]作哈尔小波变换。哈尔小波变换实际上是使用求均值和差值的方法进行分解。我们把f(x)看成是矢量空间V3中的一个向量，尺度因子i=
3，因此最多可分解为3层.分解过程如下。
步骤1：c=(2 5,8 9,7,4,-1,1)/√2=(7,17,11,0, -3, -1,3,
-2)/√2
解释：
平均值(7,17,11,0)/2
差值(-3，-1,3，-2)/2(与平均值的差值相对于直接差值的1/2  即 A- (A+B)/2 = (A-B)/2  ,我们将除以2提出来)
所以步骤一的结果是(7,17,11,0)/2  并上(-3，-1,3，-2)/2 =(7,17,11,0, -3, -1,3, -2)/2
步骤2：c=[(7+17)/√2,(11+0)/√2,(7-17)/√2,(11-0)/√2,-3,-1,3,-2]/√2
=(24/√2,11/√2,-10/√2,11/√2,-3,-1,3,-2)/√2
步骤3:c=[(24+11)/2,(24-11)/2
,-10/√2,11/√2,-3,-1,3,-2]/√2=(35/2,13/2,-10/√2,11/√2,-3,-1,3,-2)/√2
这个结果实际上是后面哈尔转换N=8转换矩阵相乘的结果。
N = 8，
。
6.二维哈尔小波变换
前面已经介绍了一维小波变换的基本原理和变换方法。这节将结合具体的图像数据系统地介绍如何使用小波对图像进行变换。一幅图像可看成是由许多像素组成的一
个矩阵，在进行图像压缩时，为降低对存储器的要求,人们通常把它分成许多小块，例如以8×8
个像素为一块，并用矩阵表示，然后分别对每一个图像块进行处理。在小波变换中，由于小波变换中使用的基函数的长度是可变的，虽然无须像以离散余弦变换
(DCT)为基础的JPEG标准算法那样把输入图像进行分块以避免产生JPEG图像那样的“块效应”，但为便于理解小波变换，还是从一个小的图像块入手，
并且继续使用哈尔小波对图像进行变换。
假设有一幅灰度图像，其中的一个图像块用矩阵 A表示为，
[64  2   3   61  60  6
7   57
9   55  54  12  13  51  50  16
17  47  46  20  21  43  42
24
40  26  27  37  36  30  31  33
32  34  35  29  28  38  39  25
41
23  22  44  45  19  18  48
49  15  14  52  53  11  10  56
8   58   59  5
4   62   63
1]
一个图像块是一个二维的数据阵列，进行小波变换时可以对阵列的每一行进行变换，然后对行变换之后的阵列的每一列进行变换，最后对经过变换之后的图像数据阵列进行编码。
从求均值(averaging)与求差值(differencing)开始。在图像块矩阵A中,第一行的像素值为[64 2  3  61  60  6  7
57]
步骤1：在第一行上取每一对像素的平均值，并将结果放到第一行的前4个位置，其余的4个数是R0行每一对像素的第一个数与其相应的平均值之差，并将结果放到第一行的后4个位置
步骤2：对第一行前4个数使用与第一步相同的方法，得到两个平均值和两个差(系数)，并依次放在第一行的前4个位置，其余的4个细节系数位置不动。
步骤3：用与第1和2
步相同的方法，对剩余的一对平均值求平均值和差值。
使用求均值和求差值的方法，对矩阵的每一行进行计算，得到A’
[32.5  0   0.5
0.5   31  -29  27   -25
32.5  0  -0.5  -0.5  -23   21  -19   17
32.5  0
-0.5  -0.5  -15   13  -11    9
32.5  0   0.5   0.5   7   -5    3
-1
32.5  0   0.5   0.5  -1    3   -5    7
32.5  0  -0.5  -0.5   9   -11
13   -15
32.5  0  -0.5  -0.5   17  -19   21   -23
32.5  0   0.5   0.5
-25  27   -29
31]
其中，每一行的第一个元素是该行像素值的平均值，其余的是这行的细节系数。使用同样的方法，对A’的每一列进行计算，得到A’’
[32.5
0  0   0    0     0    0  0
0   0  0   0    0     0    0  0
0   0
0   0    4    -4   4  -4
0   0  0   0    4    -4   4  -4
0   0
0.5  0.5  27   -25  23 -21
0   0 -0.5 -0.5  -11  9   -7  5
0   0
0.5  0.5  -5    7   -9 11
0   0 -0.5 -0.5  21  -23  25 -27]
其
中，左上角的元素表示整个图像块的像素值的平均值，其余是该图像块的细节系数。根据这个事实，如果从矩阵中去掉表示图像的某些细节系数，事实证明重构的图
像质量仍然可以接受。具体做法是设置一个阈值d，例如d>=5的细节系数就把它当作“0”看待，这样经过变换之后的矩阵就变成A’’’
[32.5   0
0  0    0   0    0  0
0   0  0  0    0   0    0  0
0   0  0  0
0   0    0  0
0   0  0  0    0   0    0  0
0   0  0  0   27  -25
23  -21
0   0  0  0   -11  9   -7  0
0   0  0  0   0   7   -9
11
0   0  0  0   21  -23  25
-27]
“0”的数目增加了18个，也就是去掉了18个细节系数。这样做的好处是可提高编码的效率。对A’’’矩阵进行逆变换，得到了重构的近似矩阵AA
[59.5
5.5  7.5  57.5  55.5  9.5  11.5   53.5
5.5   59.5  57.5  7.5  9.5   55.5
53.5  11.5
21.5  43.5  41.5  23.5  25.5  39.5  32.5  32.5
43.5  21.5
23.5  41.5  39.5  25.5  32.5  32.5
32.5  32.5  39.5  25.5  23.5  41.5  43.5
21.5
32.5  32.5  25.5  39.5  41.5  23.5  21.5  43.5
53.5  11.5  9.5  55.5
57.5   7.5   5.5   59.5
11.5  53.5  55.5  9.5  7.5   57.5   59.5
5.5]
如果矩阵A的数据与矩阵AA的数据用图表示，原图和重构图差别不是很大
转：
http://blog.csdn.net/mmmmmmmmnnnnxx/article/details/5169111
Another example:
小波变换的基本思想是用一组小波函数或者基函数表示一个函数或者信号，例如图像信号。为了理解什么是小波变换，下面用一个具体的例子来说明小波变换的过程。1.
求有限信号的均值和差值
[例8. 1]
假设有一幅分辨率只有4个像素 的一维图像，对应的像素值或者叫做图像位置的系数分别为：
[9 7 3
5]
计算它的哈尔小波变换系数。
计算步骤如下：步骤1：求均值(averaging)。计算相邻像素对的平均值，得到一幅分辨率比较低的新图像，它的像素数目变成了2个，即新的图像的分辨率是原来的1/2，相应的像素值为：
[8 4]
步
骤2：求差值(differencing)。很明显，用2个像素表示这幅图像时，图像的信息已经部分丢失。为了能够从由2个像素组成的图像重构出由4个像
素组成的原始图像，就需要存储一些图像的细节系数(detail
coefficient)，以便在重构时找回丢失的信息。方法是把像素对的第一个像素值减去这个像素对的平均值，或者使用这个像素对的差值除以2。在这个
例子中，第一个细节系数是(9-8)=1，因为计算得到的平均值是8,它比9小1而比7大1，存储这个细节系数就可以恢复原始图像的前两个像素值。使用同
样的方法，第二个细节系数是(3-4)=-1，存储这个细节系数就可以恢复后2个像素值。因此，原始图像就可以用下面的两个平均值和两个细节系数表示，
[8 4 1
-1]
步骤3：重复第1，2步，把由第一步分解得到的图像进一步分解成分辨率更低的图像和细节系数。在这个例子中，分解到最后，就用一个像素的平均值6和三个细节系数2,1和-1表示整幅图像。
[6 2 1
-1]
这个分解过程如表8-1所示。
表8-1 哈尔变换过程
分辨率
平均值
细节系数
4
[9 7 3 5]
2
[8 4]
[1 -1]
1
[6]
[2]
由此可见，通过上述分解就把由4像素组成的一幅图像用一个平均像素值和三个细节系数表示，这个过程就叫做哈尔小波变换(Haar
wavelet transform)，也称哈尔小波分解(Haar wavelet
decomposition)。这个概念可以推广到使用其他小波基的变换。
从这个例子中我们可以看到：① 变换过程中没有丢失信息，因为能够从所记录的数据中重构出原始图像。
②
对这个给定的变换，我们可以从所记录的数据中重构出各种分辨率的图像。例如，在分辨率为1的图像基础上重构出分辨率为2的图像，在分辨率为2的图像基础上重构出分辨率为4的图像。
③
通过变换之后产生的细节系数的幅度值比较小，这就为图像压缩提供了一种途径，例如去掉一些微不足道的细节系数并不影响对重构图像的理解。
========================================================================
=原理
==============================
哈尔小波变换是于1909年由Alfréd Haar所提出，是小波变换(Wavelet transform)中最简单的一种变换，也是最早提出的小波变换。他是多贝西小波的于N=2的特例，可称之为D2
哈尔小波的母小波(mother wavelet)可表示为：
且对应的缩放方程式(scaling function)可表示为：
特性
哈尔小波具有如下的特性： (1)任一函数都可以由
以及它们的位移函数所组成
(2)任一函数都可以由常函数，
以及它们的位移函数所组成
(3) 正交性(Orthogonal) 
(4)不同宽度的(不同m)的wavelet/scaling functions之间会有一个关系
φ(t) = φ(2t) + φ(2t ? 1)
ψ(t) = φ(2t) ? φ(2t ? 1)
(5)可以用 m+1的 系数来计算 m 的系数 若
哈尔变换
Haar Transform最早是由A. Haar在1910年“Zur theorie der orthogonalen funktionensysteme”中所提出，是一种最简单又可以反应出时变频谱(time-variant spectrum)的表示方法。其观念与Fourier Transform相近，Fourier Transform的原理是利用弦波sine与cosine来对信号进行调制；而Haar Transform则是利用Haar function来对信号进行调制。Haar function也含有sine、cosine所拥有的正交性，也就是说不同的Haar function是互相orthogonal，其内积为零。
以下面N=8的哈尔变换矩阵为例，我们取第一列和第二列来做内积，得到的结果为零；取第二列和第三列来做内积，得到的结果也是零。依序下去，我们可以发现在哈尔变换矩阵任取两列来进行内积的运算，所得到的内积皆为零。N = 8，
。
在此前提下，利用Fourier Transform的观念，假设所要分析的信号可以使用多个频率与位移不同的Haar function来组合而成，进行Haar Transform时，因为Haar function的正交性，便可求出信号在不同Haar function(不同频率)的情况下所占有的比例。
Haar Transform有以下几点特性：
1.不需要乘法(只有相加或加减)
2.输入与输出个数相同
3.频率只分为低频(直流值)与高频(1和-1)部分
4.可以分析一个信号的Localized feature
5.运算速度极快，但不适合用于信号分析
6.大部分运算为0，不用计算
7.维度小，使用的memory少
8.因为大部分为高频，变换较笼统
对一矩阵A做哈尔小波变换的公式为B = HAHT，其中A为一
的区块且H为N点的哈尔小波变换。而反哈尔小波变换为A = HBHT。以下为H在2、4及8点时的值：N = 2，
。N = 4，
。N = 8，
。此外，当N = 2k时，
。其中H除了第0个row为φ(φ=[1 1 1 ... 1]/
，共N个1)，第2p + q个row为hp,q且。
哈尔小波变换应用于图像压缩
由于数字图片档案过大，因此我们往往会对图片做图像压缩，压缩过后的档案大小不仅存放于电脑中不会占到过大容量，也方便我们于网络上传送。哈尔小波变换其中一种应用便是用来压缩图像。压缩图像的基本概念为将图像存成到一矩阵，矩阵中的每一元素则代表是每一图像的某画素值，介于0到255间。例如256x256大小的图片会存成256x256大小的矩阵。JPEG影像压缩的概念为先将图像切成8x8大小的区块，每一区块为一8x8的矩阵。示意图可见右图。在处理8x8二维矩阵前，先试着对一维矩阵
作哈尔小波变换，公式为
。
范例对8x8的二维矩阵A作哈尔小波变换，由于AH是对A的每一行作哈尔小波变换，作完后还要对A的每一列作哈尔小波变换，因此公式为HTAH。以下为一简单的例子:。列哈尔小波变换(row Haar wavelet transform)。行哈尔小波变换(column Haar wavelet transform)。由以上例子可以看出哈尔小波变换的效果，原本矩阵中变化量不大的元素经过变换后会趋近零，再配合适当量化便可以达到压缩的效果了。此外若一矩阵作完哈尔小波变换后所含的零元素非常多的话，称此矩阵叫稀疏，若一矩阵越稀疏压缩效果越好。因此可对定一临界值ε若矩阵中元素的绝对值小于此临界值ε，可将该元素令成零，可得到更大的压缩率。然而ε取过大的话会造成图像严重失真，因此如何取适当的ε也是值得讨论的议题。
哈尔小波变换运算量比沃尔什变换更少若应用于区域的频谱分析及侦测边缘的话，离散傅立叶变换、Walsh-Hadamard变换及哈尔小波变换的计算量见下表Running Timeterms required for NRMSE 
离散傅里叶变换9.5秒43
沃尔什变换2.2秒65
哈尔小波变换0.3秒128