第5章-小波变换的matlab实现汇总课件

第第4章章 小波变换的小波变换的matlab实现实现11.Matlab1.Matlab中小波种类中小波种类 u15种u经典类小波：Harr小波、Morlet小波、Mexican hat小波、Gaussian小波u正交小波：db小波、对称小波、Coiflets小波、Meyer小波u双正交小波u查看命令 wavemngr(read,1)2小波分析示例小波分析示例u一维连续小波一维连续小波 1.coefs=cwt(s,scale,wname)2.coefs=cwt(s,scale,wname,plot)u c=cwt(noissin,1:48,db4,plot);3uC=cwt(noissin,2:2:128,db4,plot)4u图形接口方式（GUI）u命令：wavemenu567一维离散小波分解一维离散小波分解u命令：dwtu格式：cA1,cD1=dwt(X,wname)cA1,cD1=dwt(X,Lo_D,Hi_D)u举例：load leleccum;s=leleccum(1:3920);ls=length(s);cA1,cD1=dwt(s,db1);8u原始信号u低频系数u高频系数9系数重构系数重构u命令：命令：upcoefu格式：格式：1.Y=upcoef(O,X,wname,N)2.Y=upcoef(O,X,wname,N,L)3.Y=upcoef(O,X,Lo_R,Hi_R,N)4.Y=upcoef(O,X,Lo_R,Hi_R,N,L)5.Y=upcoef(O,X,wname)6.Y=upcoef(O,X,Lo_R,Hi_R)uO=a 低频，低频，O=d 高频高频10u举例：A1=upcoef(a,cA1,db1,1,ls);D1=upcoef(d,cD1,db1,1,ls);subplot(1,2,1);plot(A1);title(Approximation A1)subplot(1,2,2);plot(D1);title(Detail D1)11逆变换恢复信号逆变换恢复信号u命令：命令：idwtu格式：格式：1.X=idwt(cA,cD,wname)2.X=idwt(cA,cD,Lo_R,Hi_R)3.X=idwt(cA,cD,wname,L)4.X=idwt(cA,cD,Lo_R,Hi_R,L)12举例：A0=idwt(cA1,cD1,db1,ls);13多尺度一维分解多尺度一维分解u命令：命令：wavedecu格式：格式：C,L=wavedec(X,N,wname)C,L=wavedec(X,N,Lo_D,Hi_D)14uC,L=wavedec(s,3,db1);15低频系数提取低频系数提取u命令：命令：appcoefu格式：格式：1.A=appcoef(C,L,wname,N)2.A=appcoef(C,L,wname)3.A=appcoef(C,L,Lo_R,Hi_R,N)4.A=appcoef(C,L,Lo_R,Hi_R)16高频系数提取高频系数提取u命令：命令：detcoefu格式：格式：1.A=detcoef(C,L,N)2.A=detcoef(C,L)17u举例 cA3=appcoef(C,L,db1,3);cD3=detcoef(C,L,3);cD2=detcoef(C,L,2);cD1=detcoef(C,L,1);18重构系数重构系数u命令：命令：wrcoefu格式：格式：1.X=wrcoef(type,C,L,wname,N)2.X=wrcoef(type,C,L,Lo_R,Hi_R,N)3.X=wrcoef(type,C,L,wname)4.X=wrcoef(type,C,L,Lo_R,Hi_R)utype=a 低频，低频，type=d 高频高频19uA3=wrcoef(a,C,L,db1,3);uD1=wrcoef(d,C,L,db1,1);uD2=wrcoef(d,C,L,db1,2);uD3=wrcoef(d,C,L,db1,3);20重构原始信号重构原始信号u命令：命令：waverecu格式：格式：1.X=waverec(C,L,wname)2.X=waverec(C,L,Lo_R,Hi_R)例子：例子：A0=waverec(C,L,db1);重构最大误差：重构最大误差：Err=max(abs(s-A0)2122图形接口方式（图形接口方式（GUI）2324252627282.二维离散小波二维离散小波u单尺度分解单尺度分解dwt2u格式：格式：1.cA1,cH1,cV1,cD1=dwt2(X,wname)2.cA1,cH1,cV1,cD1=dwt2(X,Lo_D,Hi_D)ucA1,cH1水平;cV1垂直;cD1对角u应用：load wbarb;figure(1);image(X);colormap(map);colorbar;cA1,cH1,cV1,cD1=dwt2(X,bior3.7)29重构系数重构系数u命令：命令：upcoef2u格式：格式：1.Y=upcoef2(O,X,wname,N,S)2.Y=upcoef2(O,X,Lo_R,Hi_R,N,S)3.Y=upcoef2(O,X,wname,N)4.Y=upcoef2(O,X,Lo_R,Hi_R,N)5.Y=upcoef2(O,X,wname)6.Y=upcoef2(O,X,Lo_R,Hi_R)uO：a低频；低频；h水平；水平；v垂直；垂直；d对角对角30uA1=upcoef2(a,cA1,bior3.7,1);uH1=upcoef2(h,cH1,bior3.7,1);uV1=upcoef2(v,cV1,bior3.7,1);uD1=upcoef2(d,cD1,bior3.7,1);ufigure(2);colormap(map);subplot(2,2,1);image(wcodemat(A1,192);title(Approximation A1)subplot(2,2,2);image(wcodemat(H1,192);title(Horizontal Detail H1)subplot(2,2,3);image(wcodemat(V1,192);title(Vertical Detail V1)subplot(2,2,4);image(wcodemat(D1,192);title(Diagonal Detail D1)3132二维逆变换二维逆变换u命令：命令：idwt2u格式：格式：1.X=idwt2(cA1,cH1,cV1,cD1,bior3.7);2.X=idwt2(cA1,cH1,cV1,cD1,bior3.7);3.X=idwt2(cA1,cH1,cV1,cD1,bior3.7);4.X=idwt2(cA1,cH1,cV1,cD1,bior3.7);u应用：应用：Xsyn=idwt2(cA1,cH1,cV1,cD1,bior3.7);33多尺度二维小波多尺度二维小波u命令：命令：wavedec2u格式：格式：1.C,S=wavedec2(X,N,wname)2.C,S=wavedec2(X,N,Lo_D,Hi_D)34uC,S=wavedec2(X,2,bior3.7);%图像的多尺度二维小波分解35提取低频系数提取低频系数u命令：命令：appcoef2u格式：格式：1.A=appcoef2(C,S,wname,N)2.A=appcoef2(C,S,wname)3.A=appcoef2(C,S,Lo_R,Hi_R)4.A=appcoef2(C,S,Lo_R,Hi_R,N)ucA2=appcoef2(C,S,bior3.7,2);%从上面从上面的的C中提取第二层的低频系数中提取第二层的低频系数36提取高频系数提取高频系数u命令：命令：detcoef2u格式：格式：A=detcoef2(type,C,S,wname,N)u说明：说明：Type:h 水平；水平；v垂直；垂直；d对角对角 cH2=detcoef2(h,C,S,2);cV2=detcoef2(v,C,S,2);cD2=detcoef2(d,C,S,2);cH1=detcoef2(h,C,S,1);cV1=detcoef2(v,C,S,1);cD1=detcoef2(d,C,S,1);37重构系数重构系数u命令：命令：wrcoef2u格式：格式：1.X=wrcoef2(type,C,S,wname,N)2.X=wrcoef2(type,C,S,Lo_R,Hi_R,N)3.X=wrcoef2(type,C,S,wname)4.X=wrcoef2(type,C,S,Lo_R,Hi_R,N)u说明：说明：Type:a低频；低频；h 水平；水平；v垂直；垂直；d对角对角38uA2=wrcoef2(a,C,S,bior3.7,2);H1=wrcoef2(h,C,S,bior3.7,1);%重构第重构第1、2层的层的高频信号高频信号 V1=wrcoef2(v,C,S,bior3.7,1);D1=wrcoef2(d,C,S,bior3.7,1);H2=wrcoef2(h,C,S,bior3.7,2);V2=wrcoef2(v,C,S,bior3.7,2);D2=wrcoef2(d,C,S,bior3.7,2);3940重构原始信号重构原始信号u命令：命令：waverec2u格式：格式：X=waverec2(C,S,wname)X=waverec2(C,S,Lo_R,Hi_R)u应用：应用：X0=waverec2(C,S,bior3.7);412D图形接口图形接口42显示显示43小波分析用于信号处理小波分析用于信号处理u常用信号的小波分析常用信号的小波分析u信号的特征提取信号的特征提取u信号处理信号处理uGUI进行信号处理进行信号处理44正弦波的线性组合正弦波的线性组合uS(t)=sin(2t)+sin(20t)+sin(200t)45u间断点检测u波形未来预测u各分信号的频率识别u信号从近似到细节的迁移46分段信号分段信号uS(t)=sin(0.03t)t=1:500 或sin(0.3t)t=500:1000u信号抑制u信号未来预测47信号的特征提取信号的特征提取u信号的突变点往往是它的重要特征信号的突变点往往是它的重要特征u信号的频率谱和它的幅值等表征了信号的许信号的频率谱和它的幅值等表征了信号的许多信息。多信息。u信号的连续性（即信号的奇异性）分析、信信号的连续性（即信号的奇异性）分析、信号的频率谱分析和幅值谱分析不可或缺。号的频率谱分析和幅值谱分析不可或缺。u小波分析进行特征提取时，两种处理方法，小波分析进行特征提取时，两种处理方法，即边界的处理和滤波。即边界的处理和滤波。u利用小波分析得到低频和高频部分。利用小波分析得到低频和高频部分。48检测信号的突变点检测信号的突变点u提取了信号的近似特征a和细节特征d。u在原始信号图像上，无法得知原始信号导数的不连续性。49信号的奇异点检测信号的奇异点检测u【定义定义】在某一尺度在某一尺度 下，如果存在一点下，如果存在一点 使得使得 ，则称点，则称点 是局是局部极值点，且部极值点，且 在在 上有一个模极大值（过零）点。上有一个模极大值（过零）点。如果对如果对 的某一领域内的任意点的某一领域内的任意点 ，则称，则称 为小波变换模极大值（过零）点。为小波变换模极大值（过零）点。尺度空间中所有的模极大值点的连线称为模极大值线。尺度空间中所有的模极大值点的连线称为模极大值线。50u定理：设定理：设 为一严格的整数，为一严格的整数，为具有为具有 阶消失阶消失矩、矩、次连续可微和紧支集的小波，次连续可微和紧支集的小波，（为某一实数区间）为某一实数区间）,若存在尺度若存在尺度 ，使得使得 ，没有局部极大值点，则在区间没有局部极大值点，则在区间 上是一致上是一致Lipschitz （为任一小的正数）。一为任一小的正数）。一般来讲，函数在某一点的般来讲，函数在某一点的Lipschitz指数指数 表征了表征了该点的奇异性大小，该点的奇异性大小，越大，该点的光滑度越高，越大，该点的光滑度越高，越小，该点的奇异性越大。越小，该点的奇异性越大。51u当小波函数可看做某一平滑函数的一阶导数当小波函数可看做某一平滑函数的一阶导数时，信号小波变换模的局部极值点对应于信时，信号小波变换模的局部极值点对应于信号的突变点（或边缘）；当小波函数可看做号的突变点（或边缘）；当小波函数可看做某一平滑函数的二阶导数时，信号小波变换某一平滑函数的二阶导数时，信号小波变换模的过零点，也对应于信号的突变点（或边模的过零点，也对应于信号的突变点（或边缘）。缘）。u因此，采用检测小波变换系数模的过零点和因此，采用检测小波变换系数模的过零点和局部极值点的方法可以检测信号的边缘位置。局部极值点的方法可以检测信号的边缘位置。u比较而言，采用局部边缘进行检测更具有优比较而言，采用局部边缘进行检测更具有优越性。越性。52u信号的奇异性通常可以分为两种情况：信号的奇异性通常可以分为两种情况：u第一种类型的间断点：信号在某一时刻，其第一种类型的间断点：信号在某一时刻，其幅值发生突变，引起信号的不连续，信号的幅值发生突变，引起信号的不连续，信号的突变处是间断点；突变处是间断点；u第二种类型的间断点：信号外观上很光滑，第二种类型的间断点：信号外观上很光滑，其幅值没有突变，但在信号的一阶微分上有其幅值没有突变，但在信号的一阶微分上有突变产生，且一阶微分是不连续的。突变产生，且一阶微分是不连续的。535455信号自相似性的检测信号自相似性的检测u小波系数与自相似性：小波系数与自相似性：u小波分解可通过计算信号和小波的小波分解可通过计算信号和小波的“自相似自相似指数指数”得到。得到。u自相似指数也就是小波系数，如果自相似指自相似指数也就是小波系数，如果自相似指数大，则信号的自相似程度就高，反之亦然。数大，则信号的自相似程度就高，反之亦然。u如果一个信号在不同的尺度上都相似于它本如果一个信号在不同的尺度上都相似于它本身，那么，其身，那么，其“自相似指数自相似指数”，或者小波系，或者小波系数在不同的尺度上也是相似的。数在不同的尺度上也是相似的。5657信号发展趋势的识别信号发展趋势的识别58在某一频率区间上信号的识别在某一频率区间上信号的识别59信号抑制与衰减信号抑制与衰减u消失矩：如果消失矩：如果 （其中（其中 是是小波函数）（小波函数）（）的平均）的平均值为值为0，那么该小波有，那么该小波有n+1个消失矩，并且可以利个消失矩，并且可以利用该小波对用该小波对n次多项式信号进行抑制。次多项式信号进行抑制。6061信号消噪与提取弱信号信号消噪与提取弱信号其中，为含噪信号，为有用信号，为噪声信号。u消噪的三个步骤：1.一维信号的小波分解2.小波分解高频系数的阈值量化3.一维小波重构62信号消噪处理信号消噪处理u命令：命令：wdenu格式：格式：1.sd=wden(s,tptr,sorh,scal,n,wavename)说明：说明：2.tptr指定阈值选取规则；指定阈值选取规则；3.sorh指定选取软阈值指定选取软阈值(sorh=s)或硬阈值或硬阈值(sorh=h)4.scal=one 基本模式基本模式 scal=sln 未知尺度的基本模式，且仅根据第一未知尺度的基本模式，且仅根据第一层的小波分解系数来估计噪声的层次，并只进行一层的小波分解系数来估计噪声的层次，并只进行一次估计，以此来变换阈值的尺度。次估计，以此来变换阈值的尺度。scal=mln 非白噪声的基本模式，且在每个不同非白噪声的基本模式，且在每个不同的小波分解层次上都估计噪声的层次，以此来变换的小波分解层次上都估计噪声的层次，以此来变换阈值的尺度。阈值的尺度。63u命令：命令：wdencmpu格式：格式：xd=wdencmp(opt,x,wavename,n,thr,sorh,keepapp)其中：其中：(1)opt=gbl,thr0,则阈值为全局阈值则阈值为全局阈值 opt=lvd,thr是向量，则阈值是在各层上大小不同的是向量，则阈值是在各层上大小不同的数值。数值。(2)keepapp=1,不对小波分解后的低频系数做处理不对小波分解后的低频系数做处理 keepapp=0,对小波分解后的低频系数也进行阈值化对小波分解后的低频系数也进行阈值化处理处理(3)x是待处理信号。是待处理信号。(4)xd是处理后信号。是处理后信号。(5)其他同其他同wden参数参数64小波分析进行消噪处理小波分析进行消噪处理3种方法种方法u默认阈值消噪处理。该法利用函数默认阈值消噪处理。该法利用函数ddencmp生成生成信号的默认阈值，然后利用函数信号的默认阈值，然后利用函数wdencmp进行消进行消噪处理。噪处理。u给定阈值消噪处理。在实际的消噪处理中，阈值往给定阈值消噪处理。在实际的消噪处理中，阈值往往通过经验公式获得，且这种阈值比默认阈值的可往通过经验公式获得，且这种阈值比默认阈值的可信度高。在进行阈值量化处理时可用函数信度高。在进行阈值量化处理时可用函数wthresh。u强制消噪处理。该法是将小波分解结构中的高频系强制消噪处理。该法是将小波分解结构中的高频系数全置数全置0，即滤掉所有高频部分，然后进行小波重，即滤掉所有高频部分，然后进行小波重构。这种方法比较简单，且消噪后的信号比较平滑，构。这种方法比较简单，且消噪后的信号比较平滑，但容易丢失信号中的有用成分。但容易丢失信号中的有用成分。656667阈值选取规则阈值选取规则u命令：命令：wthreshu格式：格式：yt=wthresh(y,sorh,thr)u说明：该函数根据参数说明：该函数根据参数sorh的取值返回输入的取值返回输入分解系数的软阈值或硬阈值。硬阈值对应于分解系数的软阈值或硬阈值。硬阈值对应于最简单的处理方法，而软阈值具有很好的数最简单的处理方法，而软阈值具有很好的数学特性，并且所得到的理论结果是可用的。学特性，并且所得到的理论结果是可用的。686970信号压缩信号压缩u步骤步骤1.信号的小波分解信号的小波分解2.对高频系数进行阈值量化处理。对第一到第对高频系数进行阈值量化处理。对第一到第N层的层的高频系数，均可选择不同的阈值，并用硬阈值进行高频系数，均可选择不同的阈值，并用硬阈值进行系数的量化。系数的量化。对量化的系数进行小波重构。对量化的系数进行小波重构。u压缩与消噪主要区别：第压缩与消噪主要区别：第2步。步。u有效的信号压缩方法：有效的信号压缩方法：1.对信号进行小波尺度的扩展，并保留绝对值最大的对信号进行小波尺度的扩展，并保留绝对值最大的系数；系数；2.根据分解后各层的效果来确定某一层的阈值，且这根据分解后各层的效果来确定某一层的阈值，且这些阈值是互不相同。些阈值是互不相同。7172小波分析图像处理小波分析图像处理u图像的小波分解与重构图像的小波分解与重构HLLH1LL1HL1HH1HL1HH1LH1LL2HL2LH2HH2图像小波分解示意图73 第1级 L1斜线细节第1级 L1水平细节第1级 L1垂直细节第2级 L2细节近似图象第3级 L3742c(n,m)a2(m)b2(m)2a1(n)b1(n)a1(n)b1(n)2222LLLHHLHH小波分解数据流示意图75小波重构数据流示意图2c(n,m)p2(m)p2(n)q1(n)p2(n)q1(n)222LLLHHLHHq2(m)2276777879图像边缘失真的处理图像边缘失真的处理u补零：假设在原始支撑以外的信号以零补足。补零：假设在原始支撑以外的信号以零补足。其缺点是人为地在边界处制造了不连续。其缺点是人为地在边界处制造了不连续。u边界对称化：假设通过对称的边界值复制可边界对称化：假设通过对称的边界值复制可以恢复信号和图像的原始支撑以外的信号和以恢复信号和图像的原始支撑以外的信号和图像。缺点是在边界处，人为地制造了一阶图像。缺点是在边界处，人为地制造了一阶导数的不连续性，但该方法通常对图像处理导数的不连续性，但该方法通常对图像处理非常有效。非常有效。u一阶平滑填补：假设通过简单的一阶导数外一阶平滑填补：假设通过简单的一阶导数外插（填补时，对前两个值和后两个值使用线插（填补时，对前两个值和后两个值使用线性拟合）能够从信号或图像的原始支撑之外性拟合）能够从信号或图像的原始支撑之外恢复信号或图像。该方法通常对光滑信号较恢复信号或图像。该方法通常对光滑信号较为有效。为有效。80u零阶平滑填补：假设通过简单的常数外插能够从信零阶平滑填补：假设通过简单的常数外插能够从信号或图像的原始支撑之外恢复信号或图像。对于信号或图像的原始支撑之外恢复信号或图像。对于信号延拓来说，该方法是位于左边的第一个值和右边号延拓来说，该方法是位于左边的第一个值和右边的最后一个值的重复。的最后一个值的重复。u周期性填补周期性填补1：这里假设通过周期延拓恢复信号或：这里假设通过周期延拓恢复信号或图像原始支撑以外的信号或图像。其缺点是在边界图像原始支撑以外的信号或图像。其缺点是在边界处人为地制造了不连续性。处人为地制造了不连续性。u周期性填补周期性填补2：如果信号的长度是奇数，首先给信：如果信号的长度是奇数，首先给信号加一个采样点，其值等于最后一个值，接下来对号加一个采样点，其值等于最后一个值，接下来对信号进行周期性填补（信号进行周期性填补（1），即在两端对信号进行），即在两端对信号进行最小周期延拓。对于图像，采取同样的方法。这种最小周期延拓。对于图像，采取同样的方法。这种模式可生成最小长度的小波分解，但为了确保完全模式可生成最小长度的小波分解，但为了确保完全重构，在逆变换中也应采用同样的延拓模式。重构，在逆变换中也应采用同样的延拓模式。u总结：前总结：前5种，会存在一定的冗余，因此在任意一种，会存在一定的冗余，因此在任意一种模式中，在逆变换中都能确保对信号和图像的完种模式中，在逆变换中都能确保对信号和图像的完全重构。全重构。81u补零u对称u平滑填补82图像消噪的步骤图像消噪的步骤u二维图像信号的小波分解二维图像信号的小波分解 选择合适的小波和恰当的分解层次（记为选择合适的小波和恰当的分解层次（记为N），），然后对待分析的二维图像信号然后对待分析的二维图像信号X进行进行N层分解层分解计算。计算。u对分解后的高频系数进行阈值量化。对分解后的高频系数进行阈值量化。对于分解的每一层，选择一个恰当的阈值，并对于分解的每一层，选择一个恰当的阈值，并对该层高频系数进行软阈值量化处理。对该层高频系数进行软阈值量化处理。u二维小波的重构图像信号二维小波的重构图像信号根据小波分解后的第根据小波分解后的第N层近似（低频系数）和层近似（低频系数）和经过阈值量化处理后的各层细节（高频系数）经过阈值量化处理后的各层细节（高频系数），来计算二维信号的小波重构。，来计算二维信号的小波重构。838485图像压缩图像压缩图像扫描DWT视觉加权存储量化8*8方块去冗余86图像扫描小波滤波器量化去冗余存储22小波变换图像压缩87u对图像小波分解后，可得到一系列不同分辨对图像小波分解后，可得到一系列不同分辨率的子图像（它们所对应的频率不相同）。率的子图像（它们所对应的频率不相同）。而对于图像来说，表征它的最主要的部分是而对于图像来说，表征它的最主要的部分是低频部分，而高频部分大部分点的数值均接低频部分，而高频部分大部分点的数值均接近于近于0，而且频率越高，这种现象越明显。，而且频率越高，这种现象越明显。因此，利用小波分解去掉图像的高频部分而因此，利用小波分解去掉图像的高频部分而仅保留图像的低频部分是一种最简单的图像仅保留图像的低频部分是一种最简单的图像压缩方法。压缩方法。8889 压缩前图像的大小：压缩前图像的大小：Name Size Bytes Class Attributes X 256x256 524288 double 第一次压缩图像的大小：第一次压缩图像的大小：Name Size Bytes Class Attributes ca1 135x135 145800 double 第二次压缩图像大小：第二次压缩图像大小：Name Size Bytes Class Attributes ca2 75x75 45000 double 90%装载并显示原始图像装载并显示原始图像load flujet;subplot(1,2,1);image(X);colormap(map);title(原始图像原始图像);axis square;%=%首先利用小波首先利用小波db3对图像对图像X进行进行2层分解层分解c,l=wavedec2(X,2,db3);%=%全局阈值全局阈值thr,sorh,keepapp=ddencmp(cmp,wv,X);%=%进行压缩处理：对所有高频系数进行同样的阈值量化处理进行压缩处理：对所有高频系数进行同样的阈值量化处理Xcmp,cxc,lxc,perf0,perfl2=wdencmp(gbl,c,l,db3,2,thr,sorh,keepapp);%=%图示压缩结果图示压缩结果subplot(1,2,2);image(Xcmp);colormap(map);title(压缩后的图像压缩后的图像);axis square;disp(小波分解系数中置小波分解系数中置0的系数个数百分比：的系数个数百分比：);perf0disp(压缩后保留能量百分比：压缩后保留能量百分比：);perfl291u小波分解系数中置0的系数个数百分比：perf0=89.5226u压缩后保留能量百分比：perfl2=99.983192Wdencmp函数函数%装载并显示原始图像load wbarb;subplot(1,2,1);image(X);colormap(map);title(原始图像);%=%采用默认的全局阈值对图像进行压缩thr,sorh,keepapp,crit=ddencmp(cmp,wp,X);Xc=wpdencmp(X,sorh,3,bior3.1,crit,thr,keepapp);subplot(1,2,2);image(Xc);colormap(map);title(全局阈值压缩图像);93