Mallat算法

6.9.1.1 尺度空间的有限分解

MRA框架表明f（t）∈L2（R）可分解为无穷个小波分量的直和，但在实际应用中，仅知道f（t）的近似函数。为不失一般性，可假设原信号是在a=1或j=0的分辨率下测得的，用f0（t）表示，它属于子空间V0。而子空间V0又可分解成两个子空间。因此，在MRA框架下理解为f0（t）∈V0，这样就有如下的尺度空间的有限分解表现

V0=V1⊕W1=（V1⊕W1）⊕W0=（V2⊕W2）⊕W1

=（Vj⊕Wj）⊕Wj-1⊕Wj-2⊕…⊕W1

︙

=VJ+WJ⊕WJ-1⊕WJ-2⊕…⊕W1

其中子空间及分量分别为

地球物理信息处理基础

V0的有限分解关系仅是MRA无穷分解中的一部分。因此，就子空间而言、就函数分量而言以及就频率范围而言，有限分解的含义都与MRA相同。例如，V0中的元素f0（t）是有限频率范围的，fJ（t）∈VJ是f0（t）的最低频表现，Wj中元素δj（t）是具有特定带宽的，它们互不重叠，这些频带的总和就是f0（t）的频率范围。

为了数字计算和分析处理的目的，需要将fj（t）和δj（t）用离散数据来表示。显然，｛cj，k｝∈l2是合适的，因为 表明fj（t）∈Vj和｛cj，k｝∈l2是一一对应的，或者说利用｛cj，k｝就一定可恢复fj（t）。根据同样的理由，可用数据｛dj，k｝∈l2来表示δj（t）∈Wj，这类数据｛cj，k｝和｛dj，k｝对应着特定的离散表现，j表示尺度方面的2进制规则离散，k表示整节点标记的时移离散。

6.9.1.2 分解算法

分解算法要实现的目标是：在｛φ（t-k）｝是标准正交基条件下，已知｛cj-1，k｝、｛hk｝和｛gk｝，求出｛cj，k｝和｛dj，k｝。

ψ（t）关于φ（t）的两尺度关系式（6-95）提供了一条由尺度函数φ（t）构造母小波ψ（t）的途径。根据尺度函数

地球物理信息处理基础

和小波函数

地球物理信息处理基础

得到

地球物理信息处理基础

用l取代2k+l，上式成为

地球物理信息处理基础

计算f（t）∈L2（R）与上式两端的内积，便得到如下计算小波级数系数的一个公式

地球物理信息处理基础

式中cj-1，l=＜fj，φj-1，l＞=＜f，φj-1，l＞是信号f（t）与φj-1，l（t）的内积，即cj，l=＜f，φj，l＞。信号在Vj的正交投影 ，将两尺度关系式的φj，k（t）代入定义式φj，k（t）=2-j/2φ（2-jt-k），有

地球物理信息处理基础

由此可以计算f（t）与φj，k（t）的内积，即

地球物理信息处理基础

在数学上，fj（t）是平方可积连续函数，它是L2（R）的元素；cj，l（l∈Z）是平方可和序列，它是平方可和序列矢量空间l2（Z）中的元素。 说明fj（t）与cj，k（k∈Z）有等效对应关系，在数学上称l2（Z）与L2（R）之间存在着同构关系。显然，cj，k比fj（t）更适合于在计算机上运算。cj，k称为在分辨2j下的f（t）离散逼近，因此，可以把f0（tk）=f（tk）的采样值作为｛cj，k｝的初始数据，即｛c0，k｝=｛f（tk）｝。

综合式（6-109）和（6-110），得到一般的分解公式

地球物理信息处理基础

式（6-111）的实现过程如图6-26所示，每个尺度j所存储的数据｛cj，k｝都是按整数编号的。以j尺度层为基础来观察j-1尺度层，j尺度层的采样节点编号k对应着j-1尺度层上编号为2k的采样节点，或者说。j-1尺度的采样节点是在j尺度采样节点基础上均匀加密的结果；若以j-1尺度层为基础来观察j尺度层，则j-1尺度层上隔2取样（“隔1取1”）的节点正好对应着j尺度层上的采样节点。图6-26表明了由j-1尺度向j尺度的变换过程，｛hk｝可看做滤波器的单位冲激响应（权系数）。假设｛hk｝仅有6个元素，则式（6-111）所表明的变换过程相当于把｛h-2，h-1，h0，h1，h2，h3，｝作为权值，其中心点h0对准｛cj，k｝后再作加权平均，即

地球物理信息处理基础

这一实现过程是非常快捷的，相当于先计算序列cj-1，k（k∈Z）与 的卷积，然后抽取偶数下标（隔1抽1）的卷积结果，即得到逼近级数的系数cj，k（k∈Z）。

图6-26 分解算法计算｛cj，k｝示意图

图6-26虽然仅表明由｛cj-1，k｝计算｛cj，k｝的变换过程，它同样也表明了由｛cj-1，k｝计算｛dj，k｝的变换过程，只不过要将｛hk｝换成｛gk｝而已。对正交小波而言，gk=（-1）kh-k+1是由｛hk｝决定的。

由c0，k开始，利用式（6-111）进行迭代运算，陆续计算出c1，k、c2，k等等，与此同时，利用c0，k、c1，k、c2，k等值，同样不断计算出d1，k、d2，k等小波级数系数值。

式（6-111）所表明的计算过程由算子表示会更简单些。记Aj=｛cj，k｝，Dj=｛dj，k｝。记算子H：l2➝l2，其运算意义如式（6-111）的第1式所示，即

地球物理信息处理基础

同样，记算子G：l2➝l2，其运算意义如式（6-111）的第2式所示，即

地球物理信息处理基础

采用算子表示后，式（6-111）所表明的分解算法结构如图6-27所示。它表明

Aj=HjA0；Dj=GHj-1A0

即A0经H算子j次作用后即可获得Aj，A0经H算子j-1次作用后再经G算子作用1次即可获得Dj。

图6-27 分解算法结构示意图

图6-27还表明，只要在细密采样间隔的尺度层次上给定A0，就可利用分解算法快速地获得较粗采样间隔尺度层上的有关数据Aj和Dj（0＜j≤J）。假设实际问题中A0有N个数据，则2尺度层上的A1和D1各有N/2个数据，依此类推；还假设｛hk｝和｛gk｝分别有M个数据；那么，用A0计算A1和D1共需2MN/2次运算，从2到3尺度层共需2MN/4次运算，依此类推可知，要得到Aj和Dj（0＜j≤J）共需2MN（2-1+2-2+…2-J）次运算。由此可见，分解算法是快速的。

原信号f0（t）具有下列唯一分解

地球物理信息处理基础

δj是信号在Wj（j=1，2，…，J）各子空间上的正交投影，它们是从一个较精细的逼近变成较粗略的逼近（两个逼近的分辨率相邻近）时所丢失的信息，即

地球物理信息处理基础

可将这一不断降低逼近分辨率的过程看成是“一层又一层地把信号进行剥皮”的过程。当J选得足够大时，“剥”下来的信息总和

地球物理信息处理基础

足够多，将足以精确表示原信号f0（t），而最终的逼近信号fJ（t）的分辨率已经非常低，这样反而可以把式（6-113）当成原信号的估计，而把fJ（t）看成是估计误差。这就是说，用小波级数在所有分辨率下的全部系数（j=1，2，…，J）来代替原信号，其误差fJ（t）可以任意小。按照这种解释，式（6-111）算法就是将f0（t）的信息（c0，k是它的离散表示）表示成c1，k～cJ，k等信息和一个估计误差（实际上它是在分辨率最低即2J下的逼近）cJ，k。这一过程实际上是在一次又一次地改变着正交基（或子空间）。

6.9.1.3 信号重构算法

重构算法是分解算法的逆过程。此时已知数据｛cj，k｝和｛dj，k｝（0≤j≤J），希望利用这些数据快速准确地重构出原始数据。

一般而言，相邻两分辨率下的逼近信号存在着下列关系：

地球物理信息处理基础

计算上式左右两端与φj-l，k（t）的内积，得

地球物理信息处理基础

为了便于讨论运算过程，在上式中令m=p-2k，所以信号重构公式为

地球物理信息处理基础

根据式（6-115）计算｛ ｝的过程见图6-28所示。为了理解式（6-115）和图6-28所示的计算过程，应注意下面两点：一是j-1尺度层偶数编号的采样点对应着j尺度层上的采样点，二是示意图中假设｛hk｝的有限数据为｛h-2，h-1，h0，h1，h2，h3，｝。

可以看出，计算｛ ｝可分为两步进行：

地球物理信息处理基础

第一步计算j-1尺度层上的偶数编号采样点处的｛ ｝，此时可对j尺度层采样节点循序推进，所用数据为｛h-2，h0，h2，｝，按公式

地球物理信息处理基础

计算出j-l尺度层上偶数编号节点处的相应值。

第二步计算j-1尺度层上的奇数编号的采样点处的｛ ｝，其循序推进过程仍同于第一步，但需另用｛h3，h1，h-1｝，且按公式

地球物理信息处理基础

图6-28也表明了式（6-115）中｛dj，k｝计算｛ ｝的过程，此时要采用｛gk｝，计算式应为

地球物理信息处理基础

为了更清楚地表明｛cj，p｝和｛dj，p｝，0≤ j≤J-1重构｛c0，p｝的关系，与前述的分解情况一样，记Aj=｛cj，p｝，Dj=｛dj，p｝，采用算子的记号H*：l2➝l2，其运算意义为

地球物理信息处理基础

记G*：l2➝l2，其运算意义为

地球物理信息处理基础

在这些重构算子意义下，式（6-115）可表示为

地球物理信息处理基础

于是从尺度层j=J到尺度层j=0的重构算法过程可用图6-29表示。重构算法也是快速的，现估计其计算量。设｛hm｝和｛gm｝分别有M个数据，0尺度层｛cj，0｝有N个数据，则尺度层1上｛cj，1｝有N/2个数据，用｛c1，p｝计算｛ ｝的偶节点值和奇节点值分别需要（M/2）N/2次运算；由｛c1，p｝计算｛ ｝共需MN/2次运算；｛d1，p｝计算｛ ｝也需MN/2次运算；由1尺度层数据重构0尺度层数据共需MN次运算。可依此类推出其它相邻尺度层之间采用重构算法的运算量，最后，得总运算量为

2MN（2-1+2-2+…2-J）

图6-29 信号重构算法示意图

图6-30所示的是信号分解与重构算法的计算流程示意图，图中的“ ”表示“按因子2抽取”或“隔1抽1”运算。符号“ ”表示，“按因子2内插”或“两相邻样本间内插一个零样本”。输入是c0，k，相当于原信号的离散表示；输出是不同尺度下的小波级数系数cj，k（j∈Z）以及信号在子空间VJ的正交投影cJ，k。这是在做信号分析，其分解公式见（6-111）。与此相反，当输入是dj，k（j∈Z）和cJ，k，输出是c0，k时，为信号重建或合成，此时有

地球物理信息处理基础

将合成与分析算法合起来画成图6-30所示的形式，称之为Mallat算法。

图6-30 信号分解与重构算法示意图

6.9.1.4 Mallat算法实现中的一些问题

Mallat算法是一种纯数字的快速递推算法，在使用Mallat算法时，有一些具体问题需引起注意。

（1）对正交尺度函数φ（t）而言，Mallat算法中仅需数据｛c0，k｝和｛hk｝可进行快速的分解和重构递推运算。要存储的数据为｛cJ，k｝和｛dj，k｝，0＜j≤J，这些有用的数据的存储量等于｛c0，k｝的数据存储量。特别值得强调的是，Mallat算法中隐含着两类关系，一类是关于多分辨分析方面的，例如对0＜j≤J，有

地球物理信息处理基础

fj-1（t）=fj（t）+δj（t）

地球物理信息处理基础

另一类关系是由于尺度函数φ（t）平移正交性产生的，例如

cj，k=＜f，φj，k＞；dj，k=＜δj，ψj，k＞

gk=（-1）kh-k+1，0＜j≤J

Mallat算法正是利用了这些关系，在算法实施过程中不需尺度函数φ（t）和小波函数ψ（t）的具体形式，只要求它们存在并找出｛hk｝，就可以顺利地进行分解和重构处理了。因此，只要查得正交尺度函数双尺度方程的传递系数｛hk｝，就可以应用Mallat算法了。

顺便指出，如果尺度函数φ（t）（例如样条函数）不是平移正交的，它虽然可以生成MRA，但由此构造出的样条小波ψ（t）仅关于尺度正交，没有平移正交性。此时

地球物理信息处理基础

所以Mallat算法不再适用，必须另行推导相应的分解和重构算法。

（2）初始数据｛c0，k｝的选用，在正交小波分解中，φ（t）是正交尺度函数，0尺度层上的展180开系数｛c0，k｝=＜f（t），φ0，k（t）＞，用复杂的计算来确定初始数据｛c0，k｝是不合算的，应采用变通处理办法，即简单采用｛c0，k｝=｛f（tk）｝，用细尺度层上的采样值作为初始数据｛c0，k｝。这种做法似乎有些不严密，但可以证明，虽然｛f（tk）｝作为｛c0，k｝的近似值时有微小误差，但数据｛f（tk）｝同样有效地表现了f0（t）的变化波动状况和有效的频率范围，这种替代不会影响对f0（t）的时频分析；同时还应看到，用｛f（tk）｝作初始数值，不仅可使问题简单化，而且也可使Mallat算法准确地分解和重构初始数据。总之，用｛f（tk）｝替代｛c0，k｝是实用方便的。

（3）分解层数和采样间隔的关系，这个问题主要从以下几方面考虑便可得出结论。

第一，因为最细的0尺度层的采样间隔T决定了f0（t）的频率范围，由取样效应可知，最大的频率范围为｜ω（f0）｜≤1/（2T）；同样，最粗的J尺度层的采样间隔为J=2JT，最低的频率范围为｜ω（fJ）｜≤1/（2T）=2-J｜ω（f0）｜。于是可从需要分辨的最高频率和需要分辨的最低频率这两个指标来决定最细尺度层的采样间隔和数据分解的层数。

第二，在最细的0尺度层上，应取用多少个数据才能满足J个层次的数据分解呢？在Mallat算法中，｛c1，k｝的数据量仅为｛c0，k｝数据量的一半，依此类推。同样，在J尺度层至少要取用NJ个数据才能表现低频量，于是推知，在0尺度层，至少取用N=2JNJ个数据，才能满足J个分解层次的需要。

（4）在Mallat算法的运算中，需用到所存储的数据外面的数据（见图6-31），图中实线框内数据是要存储的。在分解算法中，若｛hk｝仅有6个数据，由图6-26可知，要用到实线框外的数据c0，-2和c0，-1才能计算出c1，-1，要用到c0，N0+1、c0，N0+2和c0，N0+3才能计算出c1，N1，其它层次的情形类似。由于细密采样层（图6-31中对应着j=0的尺度层次）中的近似函数f0（t）∈L2（R），当t➝∞时，f0（t）➝0，实线框内足够多的采样数据｛c0，k｝，k=0，1，…，N0，已反应了f0（t）的基本特性，f0（t）在图6-31实线框外的数据几乎消失，因此在实际分解过程中可简单地令实线框外的数据为零。同样，重构过程也会用到实线框以外的数据（见图6-28），也可以简单地令实线框外的数据为零。

图6-31 Mallat算法涉及没有存储的数据示意图

另一种办法也是常用的。在最细尺度层上较多地取用数据，在计算过程中适当地多存储些数据，如图6-31中虚线框所示。此时应以实线框存储数据作为分解重构算法和进行数字分析的依据。

6.9.1.5 Mallat 算法所表现的频域分解特点

有限尺度空间的正交小波子空间的直和分解关系，例如，

V0=V1⊕W1=V2⊕W2⊕W1=V3⊕W3⊕W2⊕W1

在Mallat算法中是通过算子H和G来表现的。数据Aj表征fj（t）∈Vj，数据Dj表征δj（t）∈Wj，那么在Mallat算法中，Aj=HAj-1；Dj=GAj-1，从而实现了子空间的分解。这种子空间分解和算子H与G之间的关系示意于图6-32中（以V0分解成3层为例）。

图6-32 Mallat算法所确定的数据分解和子空间分解的对应关系

因为fj-1（t）和fj（t）都是有限频率范围的，fj（t）的频率范围仅是fj-1（t）的相对低频部分，δj（t）的频率范围是fj-1（t）的相对高频部分，所以δj（t）是fj-1（t）的关于高频带的分量。换成子空间来描述，Vj-1所表现的频率范围被分解为两部分，一部分是由Vj表现的低频部分，另一部分是由Wj表现的频带部分。因此，V0尺度空间的直和分解，表明把f0（t）的频率范围分解为由W1、W2、W3和V3所表现的频带，且这些频带是互不重叠的。这里所说的W2所表现的频带宽度，也就是其基函数ψ2，k（t）（窗函数）的频窗宽度。由aDω可知，W2所表现的频带宽是W1的一半，尺度指标增加，则小波子空间所表现的频带宽分半减小。由于对任意尺度而言，小波子空间所表现的时频窗面积是恒定不变的常数，所以对时频窗的时窗宽度而言，尺度指标增加，则相应的时窗宽度成倍增加。图6-33表示了正交小波分解用于分析t*处局部时域信号时，在各个频带的时频窗表现。

图6-33 正交小波分解在各频带的时频窗示意图

mallat算法是mallat提出的用于某一函数F(t)的二进小波分解与重构的快速算法，其地位相当于傅立叶变换中的FFT。

即相当于构造一定的函数空间，将信号F(t)分解到函数空间中进行一定的计算，获取你想要得到的成分，然后再重构返回原始信号。

具体的mallat算法原理很复杂，你可以在看看书或者相关文献。。希望能有所帮助，如果有高手，请帮忙指正！

---