离散小波变换算法剖析及其通用程序实现

第28卷第20期Vol. 28

No. 20

计算机工程与设计

Computer Engineering and Design

2007年10月Oct. 2007

离散小波变换算法剖析及其通用程序实现

熊智新，胡慕伊，陈朝霞，胡

摘

明

(南京林业大学化学工程学院，江苏南京210037)

要：针对小波变换工程应用软件开发需要，结合Mallat 算法原理分析介绍了离散小波变换的主要功能步骤以及程序设计技术的关键问题。算法采用Delphi 语言实现，大量数据验算表明，程序运算中间及最终结果和用Matlab 小波工具箱编写的小波变换程序执行情况完全一致。主要算法均以子函数形式给出，便于研究人员把在Matlab 中开发的小波变换应用算法成果转化为其它高级语言程序，构建独立的专用软件系统。

关键词：离散小波变换; Mallat 算法; 算法描述; 程序设计; 软件开发中图法分类号：TP311.11

文献标识码：A

文章编号：1000-7024(2007) 20-4856-04

Analysis of discrete wavelet transform algorithm and its general

implementation for programming

XIONG Zhi-xin,

HU Mu-yi,

CHEN Zhao-xia,

HU Ming

(College of Chemical Engineering, Nanjing Forest University, Nanjing 210037, China )

Abstract ：For the application software development of wavelet transform (WT ) , some key steps and techniques for programming are introduced based on Mallat algorithm analysis and descriptions. Many data are employed to check the proposed algorithm programmed by the language of Delphi. The interval and final results shows no differences from the outputs produced by WT programs with the Matlab wavelet toolbox. The main algorithms are written as sub-functions that are helpful for the researchers to convert their achievements in Matlab to other computer languages and construct the WT software for specialized application independently.

Key words ：discrete wavelet transform; Mallat algorithm; algorithm description; programming; software development

borland C++中将其编译成DLL 。这种方法虽然提高了通用性，但由于Matlab 函数及数据类型与C/C++不能兼容变换，简洁的M 文件翻译成C/C++程序后可能晦涩繁杂，不易维护改进，调用时可能会出现非法操作错误。

有鉴如此，在利用Visual C++、Delphi 等高级语言开发小波变换应用专业软件时，还是希望软件从界面到算法能以相同语言设计，实现系统的无缝集成，便于调试、修改和维护，这就需要自己编写小波变换程序。由于一般书籍没有给出小波变换通用程序实现方法，本文根据作者多年从事小波分析应用研究与软件开发实践，从程序设计角度介绍了离散小波变换的算法原理及编程中的关键技术。

0引言

小波变换作为一种新的信号处理工具，近年来在许多工程领域中已被广泛应用[1-2]。目前科技工作者使用最广泛的小波变换工具软件是工作在Matlab 环境下的Wavelet Toolbox 和Wavelab 软件包。由于Matlab 的M 文件采用解释执行，运行速度慢，界面编写复杂，算法也只能用源码的方式发行，知识产权得不到保障，不利于成品软件开发。很多研发者为了能充分利用Matlab 中研究的算法成果，提高工作效率，采用了各种Matlab 与其它应用程序接口技术[3-4]，开发可脱离Matlab 环境运行的应用软件。采用的接口技术通常有文件传递和动态链接库(DLL ) 调用两种方式。前者利用工具软件(如Mideva ) 将M 文件编译成可执行程序，然后被其它应用程序调用，处理结果通过中间文件传递。该方法虽然简便，但程序间交互性差，在数据量大时处理速度依然较慢，较长的程序执行时间可能引起读写文件错误。后者则利用Matlab 或Mideva 的C 语言编译器将M 文件转换为C/C++程序，再在Visual C++或

1离散小波变换Mallat 算法分析

利用小波变换处理实际工程中经采样系统得到的离散信

号，一般并不使用尺度函数或小波函数来运算，而是直接利用双通道滤波器组进行，小波基的选择也因此转换为相应滤波器系数选择。离散系列小波变换通常采用Mallat 算法，包括

收稿日期：2007-04-30E-mail ：[email protected]

基金项目：南京林业大学高学历人才基金项目(163030010) 。

作者简介：熊智新(1973－) ，男，湖北应城人，博士，讲师，研究方向为测试技术与信息智能处理；胡慕伊(1958－) ，男，教授，硕士生导师，研究方向为控制系统仿真与优化；陈朝霞(1972－) ，女，讲师，研究方向为控制系统优化；胡明(1959－) ，男，高级实验师，研究方向为控制仪表智能化。

－4856－

分解和重构两部分，该算法在小波分析中的地位相当于快速Fourier 变换算法在经典Fourier 分析中的地位[5]

。Mallat 分解

算法公式为

=(1

)

=

(2)

式中：——低频和高频分解滤波器，，——第

为滤波器中元素个数，为计算方便，

令

，

将公式变换如下

=+

=

+

£¬

²ãµÍͨ·Ö½â£¬Ôòʽ(3) 实

现过程如图1所示。

滤波器

1

3

5

+

1

+

2

+

3

+

4

+5

+1

层

图1

Mallat 分解算法计算

{

}

图1的计算过程相当于以

{

1

,

3

,

5

}作为靠模，

把

＋２后作加权平均，

得到第。

计算时只需换用滤波器{1, 3, 5}以同样算法进行。Mallat 算法还可从工程技术人员熟悉的双通道滤波器组角度来理解，此时式(1) 、(2) 所描述的分解过程如图2所示，

其中

分别是与

层偶数项才在

H

2

G

2

图2离散系列的一次小波分解

，

～

=

+

(5)

为看清楚运算过程，

同样令，并整理成下面的回复公式[7]

｛

～

=/2

=/

2

(6)

=

～

+

式中：～

分别是低频和高频重构滤波器，图4表示按式(6) 计算

～

的过程(假定

～

,

2

,

4

,

的

计算可分两步进行：

第1步计算

，此时可对

～

～

+1

层

４

１

５

／

２

＋１

图4

Mallat 重构算法计算

{

,

4

}，按式(7) 计算出=

2/2

++1

(7)

第2步计算

+1

，其循序推进

方式同第1步，但要用奇数靠模{

3

,

+1

=

3

/2

+

+1(8)

式(6)

中用

也按以上两步进行，但此时要采用

{

～

1

,

3

,

5

}滤波器。把计算出的

和

对应相加即得第

。从重构过程可看出，图3所示的二插值

和滤波运算通过以上两步得以实现。

2

小波变换Mallat 算法描述

2.1

程序的初始化

通过Mallat 算法分析可以看出，算法启动需要原始离散

信号系列0

～

ºÍ

4

个滤波器组，，记作：Lo_D，Hi_D，

Lo_R，Hi_R，分别用于分解与重构运算。对于正交小波(如Daubechies 系列) ，只要已知低频分解滤波器Lo_D，则其它3个

滤波器可推导出来[5]，即：Lo_R:=Wrev(Lo_D) ；Hi_R:=Qmf(Lo_R ) ；Hi_D:=Wrev(Hi_R) ，其中Wrev 为滤波器组反转函数，Qmf 为镜像二次滤波器(quadrature mirror filters, QMF ) 函数，算法描述(以类Pascal 语言[8]) 分别如下：

//滤波器反转函数

输入：n 个元素的一维滤波器数组p (0..n-1) 。输出：p (0..n-1) 首尾反转数组wrev_P(0..n-1) 。

Procedure WREV (p, wrev_P) For i=0To n-1Do

wrev_P[n-i-1]:=p [i ]; Return

//滤波器镜像函数

输入：n 个元素的一维滤波器数组p (0..n-1) 。

输出：p (0..n-1) 镜像数组qmf_P(0..n-1) 。Procedure QMF (p, qmf_P) For i=0To n-1Do {sgn:=1;

If (i Mod 2) =0Then //判断i 为偶数否。sgn:=-1;

qmf_P[n-i-1]:=sgn*p [i ];

－4857－

}Return

如果采用的小波函数不是正交的，则不满足上述关系，必须分别输入各滤波器的值。为应用方便，可在Matlab 中利用函数

[Lo_D,Hi_D,Lo_R,Hi_R]=wfilters ('wname' )

生成相应的4个滤波器组并以文件形式保存，应用时可根据用户小波函数选择设置由程序自动读取文件数据，初始化4个滤波器。

2.2信号的边界延拓

由于信号是有限长的，为避免产生边界效应，以保证原始

信号的分解和重构的精确实现，可在卷积计算前对原始信号边界延拓处理，常用方法有

[9]

：①零延拓；②光滑常数延拓；③对称延拓；④周期延拓。实际应用中，对称延拓和周期延拓使用较多。Matlab 的Wavelet Toolbox 默认采用对称延拓，为便于对照，本文介绍的算法也采用该方式。延拓算法要注意以下问题：

(1) 每层分解运算前都要对信号(原始或低频信号) 进行边界延拓，重构时不需要。

(2) 假定滤波器长度为L ，则信号两端各延拓L-2个元素。如果延拓前信号长度是奇数，则需要在右端再延拓一个元素，以保证延拓后信号始终为偶数，便于随后的二抽取运算。

算法描述如下：//对称边界延拓

输入：信号c (0..n-1) ，滤波器长度len 。

输出：延拓信号exp_c(0..k-1) ，k 由n 与len 的具体值确定。Procedure EDGE_EXPEND(c, len, exp_c) len:=len-2; k:=0;

For i:=0To ln-1Do //左边界延拓{exp_c[k ]:=c[len-i-1]; k:=k+1;}For i:=0To n-1Do //获取原始信号{exp_c[k ]:=c[i ]; k:=k+1;}

If (n Mod 2) =1Then //如奇数长加1len:=len +1;

For i:=0To len-1Do //右边界延拓{exp_c[k ]:=c[n-i-1]; k:=k+1;}Return;

2.3小波分解与重构卷积程序

Mallat 算法从理论上讲都归结为信号与相应的低通和高

通滤波器的迭代卷积，因此卷积运算是小波变换程序的核心。但从前面的分析及图1与图4可以看出，分解和重构卷积运算实际上并不相同，因此需各自编写卷积算法。

算法描述如下：//分解卷积算法

输入：延拓后原始信号或第

+1层分解结果dec_c(0..k ) ，k 由m 与n 的具体值确定。Procedure DEC_CONVOL(c, h, dec_c) ; EDGE_EXPEND(c, m, exp_c) ; //边界延拓k:=0;

－4858－

For i:=0To n-m Do If (i mod 2) =0Then {For j:=0To m-1Do If ((i+j)

dec_c[k ]:=dec_c[k ]+exp_c[i+j]*h[m-1-j ]; k:=k+1;}Return;

调用该函数时，只要把Lo_D或Hi_D分别赋值给h ，就可得到c 的

+1层低频平滑系数或小波系数dec_c(0..n ) ，重构滤

波器数组h (0..m-1) 。

输出：第

层重构分量

～

(或

) 。

2.4小波分解与重构卷积迭代运算

给定一个初始系列

｛｝

，Mallat 算法通过分解卷积迭代就可计算

{}，

{}(0

//分解卷积计算高频系数DEC_CONVOL(c,Hi_D,cd ) ; //记录分解系数ca,cd

RECORD_DCoef(ca,cd,ca_set,cd_set) ; c:=ca;}Return;

算法中函数RECORD_DCoef用来记录小波分解各层低频和高频系数ca 和cd ，以供信号分析或重构时应用，其算法可根据应用要求具体定义。

本文程序设计时，该函数的两个输出变量ca_set,cd_set定义为二维数组，按行存储结果，即每行第0号元素记录ca (或cd ) 的长度，其余元素记录ca (或cd ) 的值，并且第0行仅存储原始信号的长度。

//小波变换重构迭代运算

//输入：第L 层分解低频系数ca 和全部小波系数cd_set，重构滤波器组Lo_R与Hi_R

//输出：重构信号c

Procedure REC_WT(ca, cd_set,L,Lo_R,Hi_R,c ) ; For i:=LDownto 1Do

{REC_CONVOL(ca, Lo_R,rec_c) ; //计算

～

//二维动态数组，记录各层分解结果TMatrx2=arrayof array of double; End;

Delphi 采用的Object Pascal 在4.0版本之后就引入了动态数组概念，该数据结构在使用时不必事先说明数组的大小，只是在程序设计中为程序动态地开辟存储空间，尤其是允许创建行长度不等的动态数组，这对于在程序中按行存储各层小波变换结果极为方便。因为信号小波分解在实际运算中每层得到的系数个数并不相等，而且也不是有规律地对半递减，采用行不等的动态数组存储，则可根据每层分解得到系数的具体数目分配相应的存储空间，节约内层。

前述RECORD_DCoef函数中输出变量ca_set与d_set即被定义为TMatrx2类型数据。在工程应用时，还可根据实际需要对数组ca_set和cd_set进行操作，完成诸如小波滤波、数据压缩、信号特征提取等处理。

为便于对照，本文利用Matlab 小波工具箱编写了小波变换程序。大量验算结果表明，对同一组数据分解和重构，两个程序得到的无论是中间结果还是最终结果均完全一致，可实现对任意长度信号的小波变换，并且Delphi 程序的执行速度也远快于在Matlab 中的执行过程。这对于把已经在Matlab 上研究成熟的算法转移为Delphi 或其它可直接编译执行的语言，以便开发独立的小波变换应用软件极有意义。

len:=cd_set[i, 0]; //取出第i 层小波系数个数For j:=1To len Do //取出第i 层小波系数cd [i ]:=cd_set(i, j ) ;

REC_CONVOL(cd, Hi_R,rec_d) ; //计算//

计算

=

～

4结束语

本文从程序设计角度介绍了一维小波变换Mallat 算法的

+

len:=length (rec_c) ; For j:=0To len-1Do

rec_c[j ]:=rec_c[j ]+rec_d[j ];

len:=cd_set[i-1, 0]; //取出第i-1层系数个数//截取rec_c中左边len 个元素For j:=0To len-1Do ca [j ]:=rec_c[j ]; }c:=ca;Return;

程序中对第i 次迭代重构低频信号rec_c截取处理，以保证其和第i-1层高频系数cd 长度一致，

便于实现

=

～

基本原理和程序实现技术，以供程序设计人员开发小波变换应用软件时参考。其中涉及的关键算法均以子函数形式给出，便于研发者在应用时根据需要构建更为复杂的小波变换信号处理系统，或按面向对象程序设计思想，开发小波变换类。由于二维信号小波变换是一维小波变换的推广，可看作是对信号按行和列分别进行一维小波变换，因此利用本文所述算法构建小波变换图像处理处理程序也是不难的。

参考文献:

[1]

张和君, 张跃. 基于小波变换的心电信号综合检测算法研究[J ]. 计算机工程与设计,2006,27(20) :3831-3834.[2][3][4][5][6][7][8][9]

熊智新, 路文初, 胡上序. 小波变换和RBF 网络用于模式法分解重叠色谱峰[J ]. 浙江大学学报(工学版) ,2005,39(4) :516-521.王周益, 刘继兴, 柳长安.VC ++与Matlab 混合编程研究及开发实例[J ]. 计算机应用研究,2006,23(5) :154-155.

李善姬, 芦成刚. 采用VB 与Matlab 混合编程的数字滤波器设计[J ]. 计算机工程与设计,2006,27(18) :3486-3490.

胡昌华, 张军波, 夏军, 等. 基于Matlab 的系统分析与设计-小波分析[M ]. 西安:西安电子科技大学出版社,1999.

彭玉华. 小波变换与工程应用[M ]. 北京:科学出版社,1999. 徐长发, 李国宽. 实用小波方法[M ]. 武汉:华中科技大学出版, 2001.

黄刘生, 唐策善. 数据结构[M ]. 合肥:中国科技大学出版社, 2002. 唐远炎, 王玲. 小波分析与文本文字识别[M ]. 北京:科学出版社, 2004.

+计

算。由于卷积从信号左端开始，因此右端边界附近卷积结果误差较大，加上边界延拓效果，截去右端部分信号不会对重构精度有明显影响[9]。

3小波变换算法程序实现

上述算法采用类Pascal 语言描述，因此可以很方便地转

化成某种具体计算机语言编写的程序。本文选用Delphi 可视化编程环境实现了上述算法，因程序中涉及较多的数组运算，为便于处理以及函数之间的参数传递，预先定义两个动态数组，程序如下：

Type

//一维动态数组，记录分解信号或滤波器系数TRealArray1=arrayof double;

－4859－