Mallat小波的s手工编程算法说明(原创)
MATLAB小波变换指令及其功能介绍（超级有用）
MATLAB小波变换指令及其功能介绍（超级有用）MATLAB小波变换指令及其功能介绍1 一维小波变换的 Matlab 实现(1) dwt函数功能：一维离散小波变换格式：[cA,cD]=dwt(X,'wname')[cA,cD]=dwt(X,Lo_D,Hi_D)别可以实现一维、二维和 N 维DFT 说明：[cA,cD]=dwt(X,'wname') 使用指定的小波基函数'wname'对信号X 进行分解，cA、cD 分别为近似分量和细节分量；[cA,cD]=dwt(X,Lo_D,Hi_D) 使用指定的滤波器组 Lo_D、Hi_D 对信号进行分解。
(2) idwt 函数功能：一维离散小波反变换格式：X=idwt(cA,cD,'wname')X=idwt(cA,cD,Lo_R,Hi_R)X=idwt(cA,cD,'wname',L)函数 fft、fft2 和 fftn 分X=idwt(cA,cD,Lo_R,Hi_R,L)说明：X=idwt(cA,cD,'wname') 由近似分量 cA 和细节分量 cD 经小波反变换重构原始信号 X 。
'wname' 为所选的小波函数X=idwt(cA,cD,Lo_R,Hi_R) 用指定的重构滤波器 Lo_R 和Hi_R 经小波反变换重构原始信号 X 。
X=idwt(cA,cD,'wname',L) 和X=idwt(cA,cD,Lo_R,Hi_R,L) 指定返回信号 X 中心附近的 L 个点。
2 二维小波变换的 Matlab 实现二维小波变换的函数别可以实现一维、二维和 N 维 DFT函数名函数功能---------------------------------------------------dwt2 二维离散小波变换wavedec2 二维信号的多层小波分解idwt2 二维离散小波反变换waverec2 二维信号的多层小波重构wrcoef2 由多层小波分解重构某一层的分解信号upcoef2 由多层小波分解重构近似分量或细节分量detcoef2 提取二维信号小波分解的细节分量appcoef2 提取二维信号小波分解的近似分量upwlev2 二维小波分解的单层重构dwtpet2 二维周期小波变换idwtper2 二维周期小波反变换----------------------------------------------------------- (1) wcodemat 函数功能：对数据矩阵进行伪彩色编码函数 fft、fft2 和 fftn 分格式：Y=wcodemat(X,NB,OPT,ABSOL)Y=wcodemat(X,NB,OPT)Y=wcodemat(X,NB)Y=wcodemat(X)说明：Y=wcodemat(X,NB,OPT,ABSOL) 返回数据矩阵 X 的编码矩阵Y ；NB 伪编码的最大值，即编码范围为 0～NB，缺省值 NB＝16；OPT 指定了编码的方式（缺省值为 'mat'），即：别可以实现一维、二维和 N 维 DFTOPT＝'row' ，按行编码OPT＝'col' ，按列编码OPT＝'mat' ，按整个矩阵编码函数 fft、fft2 和 fftn 分ABSOL 是函数的控制参数（缺省值为 '1'），即：ABSOL＝0 时，返回编码矩阵ABSOL＝1 时，返回数据矩阵的绝对值 ABS(X)1. 离散傅立叶变换的 Matlab实现(2) dwt2 函数功能：二维离散小波变换格式：[cA,cH,cV,cD]=dwt2(X,'wname')[cA,cH,cV,cD]=dwt2(X,Lo_D,Hi_D)说明：[cA,cH,cV,cD]=dwt2(X,'wname')使用指定的小波基函数'wname' 对二维信号 X 进行二维离散小波变幻；cA，cH,cV,cD 分别为近似分量、水平细节分量、垂直细节分量和对角细节分量；[cA,cH,cV,cD]=dwt2(X,Lo_D,Hi_D) 使用指定的分解低通和高通滤波器 Lo_D 和 Hi_D 分解信号 X 。
小波变换基本原理
第五章 小波变换基本原理问题①小波变换如何实现时频分析？其频率轴刻度如何标定？ —尺度 ②小波发展史③小波变换与短时傅里叶变换比较a ．适用领域不同 b.STFT 任意窗函数 WT （要容许性条件） ④小波相关概念，数值实现算法多分辨率分析（哈尔小波为例） Daubechies 正交小波构造 MRA 的滤波器实现⑤小波的历史地位仍不如FT ，并不是万能的5.1 连续小波变换一．CWT 与时频分析 1.概念：⎰+∞∞--ψ=dt abt t S ab a CWT )(*)(1),( 2.小波变换与STFT 用于时频分析的区别小波 构造？1910 Harr 小波80年代初兴起 Meyer —小波解析形式80年代末 Mallat 多分辨率分析—WT 无须尺度和小波函数—滤波器组实现90年代初 Daubechies 正交小波变换90年代中后期 Sweblews 第二代小波变换3.WT 与STFT 对比举例（Fig 5–6, Fig 5–7） 二．WT 几个注意的问题1.WT 与)(t ψ选择有关 — 应用信号分析还是信号复原2.母小波)(t ψ必须满足容许性条件 ∞<ψ=⎰∞+∞-ψdw ww C 2)(①隐含要求 )(,0)0(t ψ=ψ即具有带通特性 ②利用ψC 可推出反变换表达式⎰⎰+∞∞-+∞∞-ψ-ψ=dadb ab t b a CWT a C t S )(),(11)(23.CWT 高度冗余（与CSTFT 相似）4.二进小波变换（对平移量b 和尺度进行离散化） )2(2)()(1)(2,22,,n t t a b t at n b a m mn m b a mm-ψ=ψ⇒-ψ=⇒•==--ψdt t t S n CWT d n m m m n m )(*)()2,2(,,⎰+∞∞---ψ=•=5.小波变换具有时移不变性),()(),()(00b b a CWT b t S b a CWT t S -↔-↔6.用小波重构信号 ∑∑∑∑+∞-∞=+∞-∞=+∞-∞=+∞-∞=ψψ=m n m n nm nm nm n m t dt d t S )(ˆ)(ˆ)(,,,,正交小波 中心问题：如何构建对偶框架{}n m ,ˆψ如何构建正交小波？5.2 分段逼近P1. =)(t φ逼近函数)2(2)(n t n t -→-φφ)2(2)()()(S ,1,0n t C t S n t C t nn nn -≈⇒-≈∑∑φφ 尺度21=a ⇒一般式：∑-=-≈nm m nm m a n t Ct S 2)2(2)(,2尺度φ)(,0,τS a m 逼近收敛于→∞→ 0,,0→∞→→逼近a m2.两尺度函数间关系 )12()2()(-+=t t t φφφ①张成空间满足10V V ⊂ ②两尺度空间差异在哪？ 3.表征细节的小波变换的引入很显然采样率越高，s T 越小， 逼近误差越小，采样率∞→无误差发现2)()()12(2)()()2(t t t t t t ϕφφϕφφ-=-+=⇒∑-≈⇒nn n t C S )2(2)t (,1φ 12,2+=m m n⎥⎦⎤⎢⎣⎡--+-∑∑+m m m m m t C m t C )122()22(212,12,1φφ⎥⎦⎤⎢⎣⎡---+-+-=∑∑+m m m m m t m t C m t m t C 2)()(2)()(212,12,1ϕφϕφ ∑∑-•-+-•+→++nn n mn n n t C C n t C C n m )(2)(212,12,112,12,1ϕφ001W V V ⊕=⇒ 4.推广⇓⊕⊕⊕⊕⊕=⊕⊕=⊕=⇒----012011011W W W W V W W V W V V m m0121W W W V V ⊕⊕⊕=--∞- ↑⊕⊕⊕=---m W W W V m m m m ,123,lim ,1012=↓↓⊕⊕⊕⊕⊕==↑↑∞---∞→∞V m W W W W V V m m m 逼近精度逼近精度⎭⎬⎫⎩-)2(22n t m m ϕ包含信息量决定 →形成最简单的MRA尺 度2V二．分段逼近与小波变换（哈尔小波） 1.信号的尺度逼近与小波表示 尺度逼近 ∑→-nm nm m t S n t C)()2(2,2φ 小波表示 ∑∑+∞-∞=+∞-∞=-=m n m mnm n t dt S )2(2)(2,ϕ Harr 小波2.Harr 小波特性①同一尺度平移正交性：⎰+∞∞-'-='--)()(*)(n n dt n t n t δϕϕ②尺度，平移均正交 ⎰∞+∞-''''+''='-->=<n n m m m m m m n m n m dt n t n t t t ,,2)(,,)2(*)2(2)(),(δδϕϕϕϕ ⇒⎭⎬⎫⎩⎨⎧-⇒形成正交基)2(22n t mm ϕ⎰∞+∞--=dt n t t S d mm n m )2(*)(22,ϕ影即为小波系数信号在正交基函数上投 分段逼近的推广—MRA 一．多分辨率分析含义①由内空间 ⊂⊂⊂⊂+-110m m m V V V 组成②若0V 空间尺度函数)(t ϕ平移正交：⎰+∞∞-=-)()(*)(n n t t δφφ则)(t ϕ为0V 空间尺度函数,任一函数S(t)可用表示)(t φ③成立当且仅当1)2()(+∈∈m m V t S V t S ④{}00=m mm V V 交集为⑤平方可积空间即为并集逼近m V )(lim 2R L V m m =∞→ 问题：Harr 小波构成最简单MRA⇓同尺度m 也满足⎰+∞∞-''-=)()(*)(,,n n dt t t n m n m δϕϕ 作变量替换即可证明⎰∑∞+∞--=-=dtn t t S C n t C t S n nn )(*)()()(φφ如何构造选其它具体的MRA 体系 二．正交小波函数的系统构造 1.两尺度方程引入 ①低通滤波器与尺度关系Harr 小波满足 ⎥⎦⎤⎢⎣⎡-+=-+=)12(21)2(212)12()2()(t t t t t φφφφφ∑-=⎥⎦⎤⎢⎣⎡=nn t n h th 卷积关系满足)()(2)2(212100φφ②频域反映令 )2(2)2()()()()(00w tw t w H n h φφφφ↔⇒↔↔)()(00w w H h φφ↔*⇒)()()2()()(2)2(200w w H w w w H w φφφφ==⇒即③含义a. LPF n h H 为)(,1)0(00=b ．根据MRA ，∏∞==Φ=Φ100)0()2()2()2()(k k wH w w H w φc.1)0(=Φ 2.QMF 的引入①)(t φ的尺度正交关系的频域反映⎰+∞∞-=-)()(*)(n n t t δφφ⇒↔--)()(w e n t jnw φφ 频域也正交⎰∑+∞∞-=njnw n dw e w w )()(*)(21δφφπ两边对n 求和 ⎰∑+∞∞-=⇒ninw dw e w w 1)(*)(21φφπ利用泊松求和公式∑∑+=-nnjnwn w F en f )2()(π（令)(2)(,1)(w w F n f πδ==则） 有 ∑∑+=-nnjnwn w e)2(2πδπ∑∑-=⇒nnjnwn w e)2(21πδπ⎰∑+∞∞-=-⇒ndw n w w w 1)2()(*)(πδφφ∑⎰+∞∞-=-ndw n w w 1)2()(2πδφ即：∑∑=+⇒=-knk w n w 1)2(1)2(22πφπφ② QMF 正交镜像滤波器组的导出 利用两尺度关系∑=++k k wH k w 1)2()2(20ππφ对k 分奇偶讨论1))12(2())12(2()22()22(2020=+++++++⇒∑∑nn n wn w H n w n w H πφππφπ1))12(2()2()22()2(22220=+++++∑∑nnn ww H n w wH πφππφ 1)2()2(2020=++⇒πwH w H1)2(*)()(*)()()(00002020=+++=++⇒πππw H w H w H w H w H w H ③含义a.镜像为)()(,1)(1)0(0000w H w H H H ππ+=⇒=b.功率互补条件—半带条件 )(*)()(00w H w H w P =20)(π+w H1π20)(w H3.正交小波滤波器满足的条件 ①频域关系根据0)(),(=-k x x φϕ可推出0)(*)()(*)(1010=+++ππw H w H w H w H 上式的解为 )(*)(01π+-=-w H e w H jw ②时域关系 令 ∑-=↔↔njnw e n h w H w H n h w H n h )()()()()()(0011根据)(*)1()1()()(*)1()1()(*)()1()(*)(0010010000πππ+↔--=+↔--+↔--↔-⇒---w H e n h n h w H en h w H n h w H n h jw n jwn n③易证 QMF w H 也为)(1④小波滤波器同样满足两尺度关系∏∑∞==Φ=-=20111)2()2()2()2()()2()(2)(k k kwH w H w w H w k t k h t ϕφϕ4.尺度与小波滤波器频域关系的矩阵表示⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡++⎥⎦⎤⎢⎣⎡++1001)()()()()(*)()()(11001010ππππw H w H w H w H W H w H w H w H 5.{}{}解释的与MRA t t n m n m )()(,,φϕ{}{}m nm mnm V t W t →→)()(,,φϕ 正交补 112+-⊕⊕⊕=⇒m m m W W W L⎰∑∑∞+∞-+∞-∞=+∞-∞===dtt t S d t dt S n m n m m n m n nm )(*)()()(,,,,ϕϕ例：求Harr 小波的频域尺度函数和小波函数⎥⎦⎤⎢⎣⎡-=⎥⎦⎤⎢⎣⎡=2121212110h h 解： 2)2()2()2()(11210w w Sin e w Cos e w H w k k w j k w j k •===Φ∏∏∞=∞=-+- ∑⋅⋅=-==---nwj jwjnww Sin e j e e n h w H )2()1(21)()(211 4)4()()2()2()(21w w Sin w w w H w =⇒=Φ=ϕϕ 其频域幅值图如Fig 5–13所示可发现其缺陷在于波纹太大 （原因—时域紧支撑） 例：理想LPF 也构成正交小波⎪⎩⎪⎨⎧≤=其它021)(0πw w H解：[]())1()1(2)()(00n n Sin w H IFT n h --==ππ 小波函数Sinc Sinc →•)( 三．有关小波函数的一些概念 1.小波消失矩 （vanishing moment ） 满足 阶消失矩具有则称N t N k dt t t k m k )(1,1,0,0)()(1ϕϕ-===⎰+∞∞-①母小波)(t ϕ平滑度由消失矩决定，消失矩越大，则)(w ϕ频域衰减越快)(t ϕ越平滑②消失矩越大，小波振荡程度越高 2.小波正则度（regularity ） ①定义：小波)(t ϕ的连续可导次数②正则度为n 的小波)(t ϕ具有（n +1）阶消失矩（必要条件） 四．问题讨论1.根据MRA 理论①小波和尺度函数均可由无穷频域次乘积得出，最终由)(0n h 决定 ②不关心其解析表达式2.MRA 理论 离散小波的数值实现5.4 小波变换与数字滤波器组一．时间离散小波变换的实现途径 1.不能直接对定义式离散化实现)2(2),()(),(2,,n t t S t t S d m mn m n m -==ϕϕ 令 )(采样周期→=T kT l 当m 较小时，n t m -2不为整数2.第一代小波变换：根据MRA 理论，由数字滤波器组实现3.第二代小波变换：Swelden 算法 由预测和更新滤波器进行交替提升实现 二．Mallat 算法 1.两个近似假设①∑∑∑-=+=nn m k nkn nk n m n m t dt C t S t S 1,000)()()()(ϕφ似由某一尺度空间函数近②n m C ,由采样数据直接近似 ⎰∞+∞--=dt n t t S C m m n m )2(*)(22,φm m w jnm jnw w e n t w e n t w t m----•↔-⇒↔-⇒↔2)2()2()()()()(2φφφφφφ滤波器组（Mallat 算法） （根据尺度函数和小波函数）)2(2)2(2222w en t m wjn m mm m -⋅⋅---↔-⇒φφ⎰∞+∞---⋅=⇒dw e w w S C w nj m mnm m 22,)2(*)(221φπ当分辨率m 足够高时 0)2(*→-w m φnt m m m nwj mn m m mt S n S dwe w S C --=---∞+∞--==⋅≈⇒⎰22222,)(2)2(2)(212π故可直接用样本数据取代 2.Mallat 算法 ①分解算法 a.推导⎰⎰⎰∞+∞--∞+∞-∞+∞-----=-==-dtn t t S dtn t t S dt t t S C m m m m n m n m )222(*)(2)2(*)(2)()(1121*,1,1φφφ两尺度关系 ⎰∑∞+∞--+-⋅im m dt i n t i h t S ))2(2(*)(2)(2021φ∑∑∑⎰++∞+∞->=<⋅=+-=iiin m i n m im m C i h t t S i h dti n t t S i h 2,02,020)(2)(),()(2))2(2(*2)()(φφ∑-+='i i m C n i h in i ,0)2(22同理-=-i m n m C n i h d ,1,1)2(2②重构算法a.推导（由两尺度关系，正交关系，及奇偶讨论可导出）⎪⎭⎫⎝⎛-+-=∑∑--i i i m i m n m d i n h C i n h C ,11,10,)2()2(2b.滤波器组实现（上采样+滤波）5.5 小波变换的应用一．小波地位小波曾火热一时，但小波不是万能的，在某些应用场合特别适用 小波无法求解微分方程纯数字和物理地位不如FT 二．信号检测方面应用 发动机声音中的撞击声检测傅里叶分析：时间平均作用模糊了信号局部特性 Gabor 变换 ：仍需长窗去包含振荡波形 小波变换 ： 小波基可任意窄 三．降噪应用 1.适用场合经典滤波：要求信号与噪声频率足够窄且不重合 高斯类噪声和脉冲噪声 → 宽带噪声 → 小波去噪 2.滤波效果①经典滤波：丢失波形尖锐处信息②小波降噪：基本保留波形尖锐处信息（与小波基选择有关） 3.滤波手段①传统方法：Prony 参数建模法②小波降噪b.可证明其统计最优性c.阈值比较（阈值T 可基于信号标准差得出） 硬阈值：比较n m d ,软阈值：考虑n m d ,符号，及其其它系数相关性 4.小波基选择：小波基应与主体信号量相近相似度越高，主小波系数越大，噪声系数则越小 NI 信号处理工具箱分解重构。
小波3.4
§4双正交小波及Mallat 算法1992，Cohen, Daubechies & Feaureau 提出双正交小波,具有很强的对称性（图像处理中，一般用的是双正交小波）。
详见：“Biorthogonal Bases of Compactly Supported Wavelets ”，Comm. Pure & Applied Math, 45, 485-560, 1992框架k j ,ψ，对偶框架kj ,~ψ，则： ∑∑∑∑><=><=jk k j k j jkkj k j x f x f x f )(~,)(~,)(,,,,ψψψψ 如果小波ψ及其对偶小波ψ~满足： nk m j n m k j ,,,,~,δδψψ>=< 则称为><ψψ~,为双正交小波. 双正交小波基及其Mallat 算法：ψ ψ~ k j ,ψ k j ,~ψ，各自不正交，但彼此正交，即： nk m j n m k j ,,,,~,δδψψ>=< },|)({,Z k j x k j ∈ψ，},|)(~{,Z k j x kj ∈ψ为基 ∑∑∑∑><=><=jk k j k j jkkj k j x f x f x f )(~,)(~,)(,,,,ψψψψ两个多分辨分析：⊂⊂⊂⊂-101V V V⊂⊂⊂⊂-101~~~V V V尺度函数ϕ，ϕ~ }|)({Z k k x ∈-ϕ构成0V 的Riesz 基}|)(~{Z k k x ∈-ϕ构成0~V 的Riesz 基作直和分解：11++⊕=j j j W V V 11~~~++⊕=j j j W V Vψ↔0W ，ψ~~0↔W 要求 11~++⊥j j W V ，11~++⊥j j W V假设存在ψ（ψ~）使得()}|)(~)({Z k k x k x ∈--ψψ构成()00~W W 的Riesz 基 则：正交性意味着0,~,,>=<'k j k j ψϕ 0~,,,>=<'k j k j ψϕ， k k '∀, 双尺度方程10-⊂∈V V ϕ， ∑-=kk k x h x )2(2)(ϕϕ，)2(ˆ)2()(ˆϖϕϖϖϕH =,∑-=kik ke h H ϖϖ21)(; 10-⊂∈V W ψ， ∑-=kk k x g x )2(2)(ϕψ,)2(ˆ)2()(ˆϖϕϖϖψG =,ikwkk e g w G -∑=21)(10~~~-⊂∈V V ϕ ∑-=kk k x h x )2(~~2)(~ϕϕ； )2(ˆ~)2(~)(ˆ~ϖϕϖϖϕH =,∑-=kik k eh H ϖϖ~21)(~; 10~~~-⊂∈V W ψ， ∑-=kk k x g x )2(~~2)(~ϕψ； )2(ˆ~)2(~)(ˆ~ϖϕϖϖψG =，∑-=kik ke g G ϖϖ~21)(~ 由0W ∈ψ，0~~W ∈ψ以及正交性条件有)~,(ψψ构成双正交小波的必要条件: ⎪⎪⎪⎩⎪⎪⎪⎨⎧=+++=+++=+++=+++0)()(~)()(~0)(~)()(~)(1)(~)()(~)(1)(~)()(~)(πϖπϖϖϖπϖπϖϖϖπϖπϖϖϖπϖπϖϖϖG H G H G H G H G G G G H H H H当取⎪⎩⎪⎨⎧∏+=∏+=--)()(~)(~)(ωωωωωH e G H e w G i i是必要条件为：1)(~)()(ˆ)(=+++z H z H H H ωωωω.------------------（1）引理1：如果ϕϕ~,满足⎪⎩⎪⎨⎧+≤+≤----εεωωϕωωϕ2121)1()(ˆ~)1()(ˆC C ，---------------------（2） C 为常数，0>ε，则：{}k j ,ψ与{}kj ,~ψ构成对偶框架。
(整理)小波变换课件第4章小波变换的实现技术.
第4章 小波变换的实现技术4.1 Mallat 算法双正交小波变换的Mallat 算法：设{}n h h =、{}n g g =、{}n h h =、{}n g g =为实系数双正交小波滤波器。
h ，g 是小波分析滤波器，h ，g 是小波综合滤波器。
h 表示h 的逆序，即n n h h -=。
若输入信号为n a ，它的低频部分和高频部分以此为1n a -和1n d -，小波分解与重构的卷积算法：11()()n n n na D a h d D a g --⎧⎪=*⎨=*⎪⎩ n11()()n n a Uah Ud g --=*+*先进行输入信号和分析滤波器的巻积，再隔点采样，以形成低频和高频信号。
对于有限的数据量，经过多次小波变化后数据量大减，因此需对输入数据进行处理。
4.1.1 边界延拓方法下面给出几种经验方法。
1. 补零延拓是假定边界以外的信号全部为零，这种延拓方式的缺点是，如果输入信号在边界点的值与零相差很大，则零延拓意味着在边界处加入了高频成分，造成很大误差。
实际应用中很少采用。
0121,0,,,,...,,0,0,......n s s s s -2.简单周期延拓将信号看作一个周期信号,即k n k s s +=。
简单周期延拓后的信号变为这种延拓方式的不足之处在于，当信号两端边界值相差很大时，延拓后的信号将存在周期性的突变，也就是说简单周期延拓可在边界引入大量高频成分，从而产生较大误差。
3. 周期对称延拓这种方法是将原信号在边界上作对称折叠，一般分二1）当与之做卷积的滤波器为奇数时，周期延拓信号为2）当与之做卷积的滤波器为偶数时，周期延拓信号为4. 光滑常数延拓在原信号两端添加与端点数据相同的常数。
0121,,,...,,n s s s s -0121,,,...,,n s s s s -0121,,,...,,n s s s s -0,...s 1,...,n s -01221,,,...,,,n n s s s s s --0121,,,...,,n s s s s -21012,...,,,,,...n s s s s s -321212,,,...,,,,...n n n s s s s s s ---10012,,...,,,,...n n s s s s s --10112,,,...,,,n n n s s s s s ---5. 平滑延拓在原信号两端用线性外插法补充采样值，即沿着信号两端包络线的一阶导数方向增加采样值。
小波学习之一（单层一维离散小波变换DWT的Mallat算法C++和MATLAB实现）---转载
⼩波学习之⼀（单层⼀维离散⼩波变换DWT的Mallat算法C++和MATLAB实现）---转载1 Mallat算法离散序列的Mallat算法分解公式如下：其中，H(n)、G(n)分别表⽰所选取的⼩波函数对应的低通和⾼通滤波器的抽头系数序列。
从Mallat算法的分解原理可知，分解后的序列就是原序列与滤波器序列的卷积再进⾏隔点抽取⽽来。
离散序列的Mallat算法重构公式如下：其中，h(n)、g(n)分别表⽰所选取的⼩波函数对应的低通和⾼通滤波器的抽头系数序列。
2 ⼩波变换实现过程(C/C++)2.1 ⼩波变换结果序列长度⼩波的Mallat算法分解后的序列长度由原序列长SoureLen和滤波器长FilterLen决定。
从Mallat算法的分解原理可知，分解后的序列就是原序列与滤波器序列的卷积再进⾏隔点抽取⽽来。
即分解抽取的结果长度为(SoureLen+FilterLen-1)/2。
2.2 获取滤波器组对于⼀些通⽤的⼩波函数，简单起见，可以通过Matlab的wfilters(‘wavename’)获取4个滤波器；特殊的⼩波函数需要⾃⾏构造获得。
下⾯以db1⼩波函数(Haar⼩波)为例，其变换与重构滤波器组的结果如下：[cpp]1. //matlab输⼊获取命令2. >> [Lo_D,Hi_D,Lo_R,Hi_R] = wfilters('db1')3.4. //获取的结果5. Lo_D =6. 0.7071 0.70717. Hi_D =8. -0.7071 0.70719. Lo_R =10. 0.7071 0.707111. Hi_R =12. 0.7071 -0.7071//matlab输⼊获取命令>> [Lo_D,Hi_D,Lo_R,Hi_R] = wfilters('db1')//获取的结果Lo_D =0.7071 0.7071Hi_D =-0.7071 0.7071Lo_R =0.7071 0.7071Hi_R =0.7071 -0.70712.3 信号边界延拓在Mallat算法中，假定输⼊序列是⽆限长的，⽽实际应⽤中输⼊的信号是有限的采样序列，这就会出现信号边界处理问题。
Mallat小波的s手工编程算法说明(原创)
Mallat 算法及问题Mallat 算法在小波多分辨率分析中具有极其重要的地位。
Mallat 算法中，与尺度函数)(t φ相联系的是低通滤波器)(n h ，与小波函数)(t ψ相联系的是高通滤波器)(n g 。
分解后得到离散逼近信号)(n a j [又称：尺度系数]，和离散细节信号)(n d j [又称：离散小波系数]。
本文以《小波分析及其应用》为主要参考书。
未指明情况下，均指该书。
1． 由滤波器系数h 计算滤波器系数g尺度滤波器（低通滤波器）)(n h 是核心，小波滤波器)(n g 可由)(n h 计算得到。
计算公式为)1()1()(1n h n g n --=-。
该公式的含义为：将)(n h 以0=n 翻转，得到)(n h -，再将序列右移1位，即得到)1(n h -。
再乘上符号n --1)1(，即得到)(n g 。
如图所示：可见，实现方法为：对)(n h 倒序，然后在对该序列的适当位置添加负号。
如最后（从左到右）1位乘以1-（因为其0=n ，1)1(1-=--n ），倒数第2位保持原数（因为其1-=n ，1)1(1+=--n ），以此类推。
特别注意，以上序列索引从1=n 开始（Matlab 中就是如此），而以0=n 点为中心翻转。
如果序列从0=n 开始索引，则需要作调整，原乘以1-改为乘以1+，原乘以1+的改为乘以1-。
在Matlab 中，用Orthfilt( )函数可得到各种正交小波的滤波器系数Lo_D （低频分解滤波器）和Hi_D （高频分解滤波器），另外，Lo_R 和Hi_R 分别为低频重构滤波器和高频重构滤波器。
特别注意的是，Lo_R 为Lo_D 的倒序：Lo_D = wrev(Lo_R)；Hi_R 为Hi_D 的倒序：Hi_D = wrev(Hi_R)。
wrev 即矢量倒序。
另，Wfilters( )函数可得到各种小波的滤波器系数，不限于正交小波。
Matlab 中，分解是按照公式 3.2.6（即 3.4.9）和 3.2.18（即 3.4.10），即：∑∑-=-=++k j j kj j k a n k g n d k a n k h n a )()2()()()2()(11进行的，其对应的抽取图为图3.1，即滤波器为分解滤波器h 和g 。
小波变换课件ch4 Mallat算法及二维小波47页PPT
Hale Waihona Puke 46、我们若已接受最坏的，就再没有什么损失。——卡耐基 47、书到用时方恨少、事非经过不知难。——陆游 48、书籍把我们引入最美好的社会，使我们认识各个时代的伟大智者。——史美尔斯 49、熟读唐诗三百首，不会作诗也会吟。——孙洙 50、谁和我一样用功，谁就会和我一样成功。——莫扎特
小波变换课件ch4 Mallat算法及二维 小波
1、纪律是管理关系的形式。——阿法 纳西耶 夫 2、改革如果不讲纪律，就难以成功。
3、道德行为训练，不是通过语言影响 ，而是 让儿童 练习良 好道德 行为， 克服懒 惰、轻 率、不 守纪律 、颓废 等不良 行为。 4、学校没有纪律便如磨房里没有水。 ——夸 美纽斯
小波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_Lendouble 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+1for(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;}。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
Mallat 算法及问题Mallat 算法在小波多分辨率分析中具有极其重要的地位。
Mallat 算法中，与尺度函数)(t φ相联系的是低通滤波器)(n h ，与小波函数)(t ψ相联系的是高通滤波器)(n g 。
分解后得到离散逼近信号)(n a j [又称：尺度系数]，和离散细节信号)(n d j [又称：离散小波系数]。
本文以《小波分析及其应用》为主要参考书。
未指明情况下，均指该书。
1． 由滤波器系数h 计算滤波器系数g尺度滤波器（低通滤波器）)(n h 是核心，小波滤波器)(n g 可由)(n h 计算得到。
计算公式为)1()1()(1n h n g n --=-。
该公式的含义为：将)(n h 以0=n 翻转，得到)(n h -，再将序列右移1位，即得到)1(n h -。
再乘上符号n --1)1(，即得到)(n g 。
如图所示：可见，实现方法为：对)(n h 倒序，然后在对该序列的适当位置添加负号。
如最后（从左到右）1位乘以1-（因为其0=n ，1)1(1-=--n ），倒数第2位保持原数（因为其1-=n ，1)1(1+=--n ），以此类推。
特别注意，以上序列索引从1=n 开始（Matlab 中就是如此），而以0=n 点为中心翻转。
如果序列从0=n 开始索引，则需要作调整，原乘以1-改为乘以1+，原乘以1+的改为乘以1-。
在Matlab 中，用Orthfilt( )函数可得到各种正交小波的滤波器系数Lo_D （低频分解滤波器）和Hi_D （高频分解滤波器），另外，Lo_R 和Hi_R 分别为低频重构滤波器和高频重构滤波器。
特别注意的是，Lo_R 为Lo_D 的倒序：Lo_D = wrev(Lo_R)；Hi_R 为Hi_D 的倒序：Hi_D = wrev(Hi_R)。
wrev 即矢量倒序。
另，Wfilters( )函数可得到各种小波的滤波器系数，不限于正交小波。
Matlab 中，分解是按照公式 3.2.6（即 3.4.9）和 3.2.18（即 3.4.10），即：∑∑-=-=++k j j kj j k a n k g n d k a n k h n a )()2()()()2()(11进行的，其对应的抽取图为图3.1，即滤波器为分解滤波器h 和g 。
其中，)()()()(n g n g n h n h -=-=。
公式3.4.9/10可转换为公式3.4.11/12，令n k k 2'-=，得n k k 2'+=，即可。
所以，本设计与Matlab 中的滤波器组对应关系应当为：D Hi g DLo h RHi g R Lo h ____⇒⇒⇒⇒（P46）2． 一维小波分解（DWT ）在程序中采用的是∑∑+=+=++k j j kj j k n a k g n d k n a k h n a )2()()()2()()(11进行分解，即书上公式（3.4.11）和（3.4.12）。
因此本次的设计中仅出现重构滤波器对)(n h 和)(n g 。
在小波一级分解时，0=j 。
这里只讨论)(1n a j +部分。
设滤波器h 的长度为L ，即：)1(),1(),0(-L h h h ；输入序列j a 的长度为N ，即：)1(,),1(),0(-n x x x 。
)14()1()14()1()04()0()2()12()1()12()1()02()0()1()1()1()1()1()0()0()0()12()1()0()2()1()1()2()0()1()14()1()2()2()3()1()4()0()2()2()()(000100010001000010000101-+-+++++=-+-+++++=--+++=-+--+++-+-=--+--++-+-+-=-+=∑L a L h a h a h a L a L h a h a h a L a L h a h a h a L a L h a h a h a h a L a L h a h a h a h a k n a k h n a k上式中，0a 的负标号，和0a 的大于等于N 的标号均为超出边界的元素。
以上每一式可以表示为矩阵相乘，即：⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡---⨯- )2()3()4()]1(),...,1(),0([000a a a L h h h h 以L=5，N=8为例⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡----⨯=-)0()1()2()3()4()2(000001a a a a a h a ， ⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡--⨯=-)2()1()0()1()2()1(000001a a a a a h a ， ⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⨯=)4()3()2()1()0()0(000001a a a a a h a ， ⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⨯=)6()5()4()3()2()1(000001a a a a a h a ， ⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⨯=)8()7()6()5()4()2(000001a a a a a h a ， ⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⨯=)10()9()8()7()6()3(000001a a a a a h a 以上可见：)4(0-a ~)1(0-a 和)8(0a ~)10(0a 为边界的延拓。
边界延拓请看下一节。
Matlab 中对应的函数为：[cA,cD] = dwt(X, 'wname') —— 'wname'为小波函数，确定滤波器系数。
[cA,cD] = dwt(X, 'wname', 'mode', MODE) —— 'mode'为延拓方式，见下一节。
[cA,cD] = dwt(X, Lo_D, Hi_D) —— 直接给定分解滤波器系数。
3． 边界的延拓首先确定需要延拓的位数。
包括前补部分，和后补部分。
滤波器h 的长度为L ，即：)1(),1(),0(-L h h h ；输入序列j a 的长度为N ，即：)1(,),1(),0(-n x x x 。
补充序列，到满足上述矩阵相乘即可。
关注0a 的索引号。
与)0(h 相乘的0a 索引号为偶数，如 ,4,2,0,2,4,--形式。
在每个矩阵乘法中，列矩阵0a 从该索引号开始，顺序1+，列矩阵长度为h 的长度L 。
明显可以看出以下延拓规律：L 为奇数：前补L-1位L 为偶数：前补L-2位N 为奇数：后补L-1位N 为偶数：后补L-2位延拓方式主要有：补零延拓（Matlab 中的ZPD ），周期延拓（PPD ），对称延拓（SYM ）。
（1）如采用补零延拓，则延拓后的序列为：000)7()6()5()4()3()2()1()0(000000000000a a a a a a a a算法：直接构造两个指定长度的、所有元素值为0的前补和后补序列，分别插入（Insert Array ）到原始序列的头和尾。
（2）如采用周期延拓，则延拓后的序列为：)2()1()0()7()6()5()4()3()2()1()0()7()6()5()4(000000000000000a a a a a a a a a a a a a a a算法：对于后延拓，先构造一个指定长度的、所有元素值为0的后补序列，然后利用For 循环（次数为延拓序列的长度）逐位修改该序列中的数值，规律为循环次数i （即延拓序列的索引数）除以原始序列长度，以其余数（R ）索引原始序列，用得到的原始序列中的数值，修改延拓序列中第i 个索引元素的值。
对于前延拓，采用以上类似的方法。
为方便处理，程序中用余数（R ）索引原始序列倒序后的序列，用循环对每个元素修改完成后，再对输出序列倒序，得到正确的后延拓序列。
前延拓序列和后延拓序列分别插入到原始序列的头和尾。
（3）如采用对称延拓，则延拓后的序列为：若边界对称点不参与对称行为（Half - Point ）则为：)5()6()7()7()6()5()4()3()2()1()0()0()1()2()3(000000000000000a a a a a a a a a a a a a a a 算法：采用周期延拓类似的方式。
不同之处在于，是用余数R 去索引原始序列，还是原始序列的倒序序列。
这时就要用上商（IQ ）的奇偶值了。
对于后延拓序列，若IQ 为奇数，则索引原始序列；若IQ 为偶数，则索引原始序列的倒序序列。
对于前延拓序列，类似考虑。
若边界对称点参与对称行为（Whole - Point ）则为：)4()5()6()7()6()5()4()3()2()1()0()1()2()3()4(000000000000000a a a a a a a a a a a a a a a 算法：采用以上类似的方法。
不同之处在于，是用余数R 去索引的原始序列必须减少一位，)0(x 或)1(-N x ，根据IQ 的奇、偶不同，以及前、后延拓不同。
设：延拓后序列长度为EXN ，如上例为15。
在编程中，延拓后的序列索引从0开始，顺序从偶数索引元素（0，2，4，…，L EXN -）开始，在延拓后的序列中取出L 长的向量，与h 向量或g 向量相乘。
向量相乘用点积（Dot Product ）。
在用Matlab 进行DWT 对比试验时，一定要用dwtmode(‘mode ’)命令更改DWT Extension Mode 为VI 中的对应模式。
4． 一维小波的重构（IDWT ）重构公式：)()2()()2()(11k d k n g k ak n h n a j k j kj ++∑∑-+-=，为便于编程实现，用k n k 2'-=代换，得：∑∑-+-=++'1'1)2'()'()2'()'()(k j k j j k n d k g k n a k h n a 。
若仅考虑：j = 0，即1级重构。
有：∑∑-+-='1'10)2'()'()2'()'()(k k k n d k g k n a k h n a )255()5()245()4()235()3()225()2()215()1()205()0()5()244()4()234()3()224()2()214()1()204()0()4()233()3()223()2()213()1()203()0()3()222()2()212()1()202()0()2()211()1()201()0()1()200()0()0(111111011111011110111011010-+-+-+-+-+-=+-+-+-+-+-=+-+-+-+-=+-+-+-=+-+-=+-=a h a h a h a h a h a h a a h a h a h a h a h a a h a h a h a h a a h a h a h a a h a h a a h a上式中，只考虑h 项。