《小波基础》PPT课件.ppt

数字信号处理,11 小波基础,什么是小波 小波发展历史 一维连续小波变换 一维离散小波变换 MATLAB小波工具箱,基本概念,引言,Amplitude,Time,Frequency,(,a),Amplitude,Time,Time Domain,(,c),信号时频域关系图,傅里叶变换：将时域 频域，使信号的频率特性一目了然。,引言,实质：任何信号在一定程度上均可表示为一些列正弦波之和。虽然傅里叶分析自诞生以来，在科学与工程领域发挥了巨大的作用，但也有明显的不足。,引言,从本质上讲，傅里叶变换就是一个棱镜，它把一个信号函数分解为众多的频率成分，这些频率又可以重构原来的信号函数。这种变换是可逆的，且保持能量不变。,引言,在音乐信号中人们关心的是什么时刻演奏什么样的音符；对地震波的记录人们关心的是什么位置出现什么样的反射波；图像识别中的边缘检测关心的是信号突变部分的位置，即纹理结构。,这些，傅里叶变换都不能完成，需要引入时-频局部化分析。,在不少实际问题中，我们关心的是信号在局部范围中的特征,例如：,短时傅里叶变换,引言,【例如】歌声可看作为声音震荡的波函数，而傅里叶变换即是将这个波函数转化成某种乐谱。不过，傅里叶变换无法反映信号在哪一时刻有高音，在哪一时刻有低音，因而所有的音符都挤在了一起。,傅里叶变换的实质，是一种全局的变换，要么完全在时间域，要么完全在频率域，因此无法表述信号的时频局部性质，而时-频局部性质恰好是非平稳信号最基本和最关键的性质。,一叶障目不见泰山,引言,需要注意的是：线性系统理论中的傅里叶变换是以在两个方向上都无限伸展的正弦曲线波作为基函数的。对于瞬态信号或高度局部化的信号（例如边缘），由于这些成分并不类似于任何一个傅里叶基函数，它们的变换系数（频谱）不是紧凑的，频谱上呈现出一幅相当混乱的构成。在这种情况下，傅里叶变换是通过复杂的安排，以抵消一些正弦波的方式构造出在大部分区间都为零的函数而实现的。因此：,傅里叶分析不能刻画时域信号的局部特性；傅里叶分析对非平稳信号的处理效果不是很好。,为了克服上述缺陷，使用有限宽度基函数的变换方法逐步发展起来了。这些基函数不仅在频率上而且在位置上是变化的，它们是有限宽度的波小波。,引言,小波：克服了傅里叶变换的缺点，不仅能检测到高音与低音，而且还能将高音与低音发生的位置与原始信号相对应。,管中窥豹略见一斑,引言,引言,11.1 什么是小波,什么是小波？所谓小波，即小区域的波，是一种长度有限、均值为零的波。,均值为0,假设 的傅里叶变换为，若满足下式：,例如：,则称 为基小波或母小波。,完全重构条件或恒等分辨条件,例 1,【例】判断Haar 小波是否为基小波。,【解】,例 1,是基小波,设有基小波，可通过尺度（伸缩）因子和位移因子来产生一系列小波基函数：,11.1 什么是小波,其中：a，是尺度因子，b，是位移因子。,11.1 什么是小波,小波的基本特点：“小”：在时域具有紧支集或近似紧支集，具有很强的衰减性；“波动性”：正负交替的“波动性”，即直流分量为零；既能在时域刻画信号的局部性，也能在频域反映信号的局部性。而且由于对高频成分采用逐渐精细的时域或频域取样步长，从而可以聚焦到对象的任何细节，所以被称为“数学显微镜”。基本思想：从函数分解的角度看，主要目标是希望能找到另外一种基函数(t)来代替傅里叶变换中的sin(t)函数，使得任何信号f(t)都能由以函数(t)为基底、经过伸缩或平移产生的一些列函数的线性组合来表示。,11.1 什么是小波,Haar小波,在小波分析中最早用到的一个具有紧支撑的正交小波函数，同时也是最简单的一个函数，它是非连续的，类似一个阶梯函数。,wname=haar;phi,psi,xval=wavefun(wname,20);plot(xval,psi);title(haar);,11.1 什么是小波,2.Daubechies小波,由世界著名的小波分析学者Inrid Daubechies构造。除了db1(即haar小波)外，其他小波没有明确的表达式，但转换函数h的平方模是很明确的。Daubechies小波函数提供了比Haar组更有效的分析和综合。Daubechies系中的小波基记为dbN，N为序号，且N1，2，10。,11.1 什么是小波,3.Coiflet小波系Coiflet函数也是由Daubechies构造的一个小波函数，它具有coifN(N1，2，3，4，5)这一系列。Coiflet具有比dbN更好的对称性。从支撑长度的角度看，coifN具有和db3N和sym3N相同的支撑长度；从消失矩的数目来看，coifN具有和db2N和sym2N相同的消失矩数目。,11.1 什么是小波,4.Morlet小波,Morlet小波不存在尺度函数;快速衰减但非紧支撑.,Gabor 小波,Morlet小波,11.1 什么是小波,5.高斯小波,这是高斯函数的一阶导数，在信号与图象的边缘提取中具有重要的应用。主要应用于阶梯型边界的提取。,特性：指数级衰减，非紧支撑；具有非常好的时间频率局部化；关于0轴反对称。,6.Marr小波,这是高斯函数的二阶导数，在信号与图象的边缘提取中具有重要的应用。主要应用于屋脊型边界和Dirac边缘的提取。,（也叫墨西哥草帽小波）,特性：指数级衰减，非紧支撑；具有非常好的时间频率局部化；关于0轴对称。,11.1 什么是小波,11.1 什么是小波,7.Meyer小波,它的小波函数与尺度函数都是在频域中进行定义的。具体定义如下：,11.1 什么是小波,8.Shannon小波,在时域，Shannon小波是无限次可微的，具有无穷阶消失矩，不是紧支的，具有渐近衰减性但较缓慢；在频域，Shannon小波是频率带限函数，具有好的局部化特性。,11.1 什么是小波,9.Battle-Lemarie样条小波,Battle-Lemarie线性样条小波及其频域函数的图形,下面不是小波,11.1 什么是小波,11.2 小波技术发展史,1807:傅里叶 傅里叶理论指出，一个信号可表示成一系列正弦或余弦函数之和，叫做傅里叶展开式,小波变换(wavelet transform)是什么老课题：函数的表示方法 新方法：傅里叶Haarwavelet transform,11.2 小波技术发展史,只有频率分辨率，而没有时间分辨率；可确定信号中包含哪些频率的信号，但不能确定具有这些频率的信号出现在什么时间。,其中，,缺憾：,11.2 小波技术发展史,1909:Alfred HaarAlfred Haar对在函数空间中寻找一个与傅里叶类似的基非常感兴趣。1909年，发现并使用了小波，后来被命名为哈尔小波(Haar wavelets)。,11.2 小波技术发展史,1945:Gabor开发了STFT(short time 傅里叶变换),这里：,11.2 小波技术发展史,1980：Morlet20世纪70年代，在法国石油公司工作的年轻地球物理学家Jean Morlet提出小波变换(wavelet transform，WT)的概念；20世纪80年代,开发了连续小波变换(continuous wavelet transform，CWT)1986：Y.Meyer法国科学家Y.Meyer与其同事创造性地构造出具有一定衰减性的光滑函数，用于分析函数；用缩放(dilations)与平移(translations)均为2 j(j0的整数)的倍数构造了L2(R)空间的规范正交基，使小波分析得到发展。,11.2 小波技术发展史,1988：Mallat算法法国科学家Stephane Mallat提出多分辨率概念，从空间上形象说明小波的多分辨率的特性，并提出了正交小波的构造方法和快速算法，称为Mallat算法；该算法统一了在此之前构造正交小波基的所有方法，包括语音识别中的镜向滤波，图象处理中的金字塔方法，地震分析中短时波形处理等。其地位相当于快速傅里叶变换在经典傅里叶分析中的地位。,11.2 小波技术发展史,小波理论与工程应用Inrid Daubechies于1988年最先揭示了小波变换和滤波器组之间的内在关系，使离散小波分析变成为现实；Ronald Coifman和Victor Wickerhauser等著名科学家在把小波理论引入到工程应用方面做出了极其重要贡献；自从发现滤波器组与小波基函数有密切关系之后，小波分析在信号处理中得到极其广泛的应用。,小波变换：逆小波变换：,11.3 连续小波变换,其中：,一、什么是连续小波变换,是卷积吗？,含义：类似于用镜头观察信号。代表镜头。则b相当于使镜头相对于目标的平行移动，而a相当于镜头向目标推进或远离。,物理意义：时域上的意义：数学显微镜，一组有效宽度不同的窗口傅里叶变换的汇集。,11.3 连续小波变换,频域上的意义：若f(t)的傅里叶变换为F()，根据Parseval定理。,11.3 连续小波变换,11.3 连续小波变换,若设 的中心为，半径为，则有：,同理：的中心为，半径为，,的中心为：半径为：,恒Q性质（品质因数）,11.3 连续小波变换,小波母函数在频域具有带通特性，其伸缩和平移系列就可以看作是一组带通滤波器。,所研究频带的的中心为：,通带的宽度为：,通常，将通带宽度与中心频率的比值称为带通滤波器的品质因数（Q），即,即带宽与中心频率的比与中心频率的位置无关,Q恒定,11.3 连续小波变换,分析高频成分,分析低频成分,（1）伸缩在时间轴上对信号进行压缩和伸展，如图所示。,（2）平移：小波函数在时间轴上的波形平行移动。如图所示。,11.3 连续小波变换,时刻1,时刻2,原始信号,二、基本步骤选择一个小波函数，并将其与要分析的信号起始点对齐；计算在这一时刻信号与小波函数的逼近程度，即计算变换系数C。C越大，就意味着此刻信号与所选择的小波函数波形越相近。,11.3 连续小波变换,变换系数依赖于所选择的小波。因此，为了检测某些特定波形的信号，应该选择波形相近的小波进行分析。,调整参数b，调整信号的分析时间段，向右平移小波，重复步骤，直到分析时段已经覆盖了信号的整个区间。,11.3 连续小波变换,b是位移因子。,调整参数a，尺度伸缩，重复步骤。,重复步骤，计算完所有尺度的连续小波变换系数。,11.3 连续小波变换,a是尺度因子。,11.3 连续小波变换,以上5步表示了下式的求解过程：,11.3 连续小波变换,示意图：,三、分析（1）尺度与频率的关系,小尺度：压缩的小波快速变化的细节 高频部分大尺度：拉伸的小波缓慢变化的粗部 低频部分,11.3 连续小波变换,11.3 连续小波变换,（2）计算复杂性小波变换比快速傅里叶变换还要快一个数量级。信号长度为M时，傅里叶变换和小波变换的计算复杂性分别如下公式：,MATLAB函数,COEFS=cwt(S,SCALES,wname)：在尺度SCALES下计算向量一维小波系数；COEFS=cwt(S,SCALES,wname,plot)：计算小波系数，并以图形显示；COEFS=cwt(S,SCALES,wname,PLOTMODE)：计算并画出连续小波变换的系数，并使用PLOTMODE对图形着色。PLOTMODE值的含义：lvl:scale-by-scale着色模式glb:考虑所有尺度的着色模式abslvl或lvlabs:使用系数绝对值的scale-by-scale着色模式absglb或glbabs:使用系数绝对值并考虑所有尺度的着色模式,MATLAB函数,【例如】COEFS=cwt(S,SCALES,wname,plot)相当于：COEFS=cwt(S,SCALES,wname,absglb)COEFS=cwt(S,SCALES,wname,PLOTMODE,XLIM)能够计算并画出连续小波变换的系数。系数使用PLOTMODE和XLIM进行着色。其中：XLIM=x1,x2，并且有如下关系：1=x1=x2=length(S),例 2,【例】已知一信号f(t)3sin(100pt)2sin(68pt)5cos(72pt)，且该信号混有白噪声，请对该信号进行连续小波变换。【解】clear all;t=0:0,01:1;f=3*sin(100*pi*t)+2*sin(68*pi*t)+5*cos(72*pi*t)+randn(1,length(t);coefs=cwt(f,1:0.2:3,db3,plot);title(对不同的尺度小波变换系数值);Ylabel(尺度);Xlabel(时间);,例 2,横坐标：变换的系数纵坐标：尺度灰度颜色越深，表示系数的值越大。,11.4 离散小波变换,连续小波变换存在的不足连续小波变换中含有很多冗余信息，冗余信息不利于对信号的分析和处理；另外，由于连续小波变换中有冗余信息，可能对尺度和平移参数进行离散化后仍可重构信号；连续小波变换的计算量也大。“离散化”要解决的核心问题尺度和平移参数要怎样离散化？尺度和平移参数离散化后，要想重构信号对小波函数应有什么样的要求？,11.4 离散小波变换,参数的离散化尺度参数a的离散化a=a0j,j Z 通常取a0的值为2，称为二进小波平移参数b的离散化取决于尺度参数b=ka0j,j,k Z,11.4 离散小波变换,离散小波变换（Discrete Wavelet Transform，DWT）缩放因子和平移参数都选择（j 0，整数）的倍数，这种变换称为双尺度小波变换。,离散小波变换,重构公式：,变换系数：,11.4 离散小波变换,执行DWT的有效方法：Mallat算法Mallat在1988年开发出一种滤波器，如下图所示，其中，S表示原始的输入信号，通过两个互补的滤波器产生A和D两个信号。,11.4 离散小波变换,特别注意：在使用滤波器对真实的数字信号进行变换时，得到的数据将是原始数据的两倍。例如，如果原始信号的数据样本为1000个，通过滤波之后每一个通道的数据均为1000个，总共为2000个。于是，根据 Nyquist采样定理提出降采样方法，即在每个通道中每两个样本数据中取一个，得到的离散小波变换的系数分别用cD和cA表示。,（a）“近似”基函数 A（b）“细节”基函数 D低频滤波系数 高频滤波系数其中：,一、Haar小波基函数,补充 1：一维Mallat算法,信号s可表示为小波近似分解 a 与小波细节分解d之和：s=a+d其中：近似分解 a：=小波近似系数 wa 基函数 A细节分解 d：=小波细节系数 wd 基函数 D,二、信号分解,补充 1：一维Mallat算法,正变换：原始信号在小波基上，获得“小波系数”分量反变换：所有“小波分解”合成原始信号例如：小波分解 a=小波系数 wa 小波基A,小波变换就是将原始信号s变换成小波系数W：w=wa,wd 近似系数wa：平均成分（低频）细节系数wd：变化成分（高频）,补充 1：一维Mallat算法,信号 s 有M个样本，J 级小波变换：,小波分解,小波系数,三、离散小波变换公式,正变换,反变换,其中：是小波基函数,补充 1：一维Mallat算法,设AJ(n)长度为L，从s(n)中取等长的L个样本值与AJ(n)对应元素相乘，再求部分和。结果为M/L个元素。,补充 1：一维Mallat算法,A 函数,D 函数,由于Haar是正交变换，除以常数的目的是使变换后平方和不变。例如：,Haar小波，例子：16点信号:6 5 9 8 3 7 8 5 6 5 9 8 1 3 3 9小波正变换：小波系数小波近似系数（加）小波细节系数（减）小波反变换：可以由分解信号恢复原始信号近似分解细节分解,四、一维信号小波变换的例子,补充 1：一维Mallat算法,原始信号：,补充 1：一维Mallat算法,小波变换,补充 1：一维Mallat算法,wd1,wa1,wa2,wd2,wd1,补充 1：一维Mallat算法,代码实现：x=6 5 9 8 3 7 8 5 6 5 9 8 1 3 3 9;c,l=wavedec(x,2,haar);%多层小波分解ca3=appcoef(c,l,haar,2)%提取第2级近似系数cd2=detcoef(c,l,2)%提取第2级细节系数cd1=detcoef(c,l,1)%提取第1级细节系数结果如下：ca3=14.0000 11.5000 14.0000 8.0000cd2=-3.0000-1.5000-3.0000-4.0000cd1=0.7071 0.7071-2.8284 2.1213 0.7071 0.7071-1.4142-4.2426,结果基本相同。,2级近似分解（原始信号每4个平均值）2级细节分解（原始信号每2个平均的差值）1级细节分解（原始信号单数和双数的差值）恢复信号,补充 1：一维Mallat算法,信号分解,函数 离散小波变换,function cA,cD=mydwt(x,lpd,hpd,dim)%对x进行一维离散小波分解%lpd：低通滤波器%hpd：高通滤波器%dim：小波分解级数%cA：平均部分的小波分解系数%cD：细节部分的小波分解系数cA=x;%初始化cD=;for i=1:dim cvl=conv(cA,lpd);%低通滤波 dnl=downspl(cvl);%降抽样，求出平均部分的分解系数 cvh=conv(cA,hpd);%高通滤波 dnh=downspl(cvh);%降抽样，求出本层分解后的细节部分系数 cA=dnl;%降抽样后的平均部分系数进入下一层分解 cD=cD,dnh;%将本层分解所得的细节部分系数存入序列cDend,函数 降采样,function y=downspl(x)%对输入序列进行降采样%降抽样是对输入序列取其偶数位，舍弃奇数位。例如%x=x1,x2,x3,x4,x5，则 y=x2,x4。N=length(x);%读取输入序列长度M=floor(N/2);%输出序列的长度是输入序列长度的一半i=1:M;y(i)=x(2*i);,MATLAB函数,降采样函数dyaddown：Y=dyaddown(X,EVENODD)Y=dyaddown(X)Y=dyaddown(X,EVENODD,type)Y=dyaddown(X,type,EVENODD),MATLAB函数,cA,cD=dwt(X,wname)使用小波wname对信号X进行单层分解，求得的近似系数存放在数组cA中，细节系数存放在数组cD中。C,L=wavedec(X,N,wname)利用小波wname对信号X进行多层分解cA,cD=dwt(X,wname)和C,L=wavedec(X,1,wname)功能是等效的，但是返回的结果不一样。cA,cD=dwt(X,wname)中返回的cA，cD分别存放是信号的近似和细节，C,L=wavedec(X,1,wname)中返回的近似和细节都存放在C中，L存放是近似和各阶细节系数对应的长度 A=appcoef(C,L,wname,N)利用小波wname从分解系数C,L中提取第N层近似系数,11.5 小波分解树,原始信号通过一对滤波器进行的分解叫做一级分解。信号的分解过程可以迭代，即可进行多级分解。一、小波分解树用下述方法分解形成的树：对信号的高频分量不再继续分解，而对低频分量连续进行分解，得到许多分辨率较低的低频分量。,11.5 小波分解树,11.5 小波分解树,二、小波包分解树不仅对信号的低频分量连续进行分解，而且对高频分量也进行连续分解，这样不仅可得到许多分辨率较低的低频分量，而且也可得到许多分辨率较低的高频分量。,例 3：小波包树,load noisbump;%读信号x=noisbump;t=wpdec(x,3,db2);%3层小波包分解fig=plot(t);%显示小波树结构,在菜单“Node Label”下选择“Index”，然后点击左边的某节点，可显示其波形。在菜单“Node Action”下选择“Split/Merge”，然后点击某节点可合并节点，再在菜单“Node Action”下选择“Visualize”，再点击节点，以显示其信号波形。,例 3：小波包树,newt=plot(t,read,fig);%获取新的树newt=wpjoin(newt,3);%从命令行修改新树fig2=plot(newt);%显示,例 3：小波包树,查看颜色表示的小波包系数wpviewcf(t,1);,11.6 小波重构,所谓重构，是指把分解的系数还原成原始信号，也称“合成”，实质是逆离散小波变换(IDWT)。包含升采样和滤波两个过程，见下图。,11.6 小波重构,升采样是在两个样本数据之间插入“0”，目的是把信号的分量加长，其过程见下图。,11.6 小波重构,重构滤波器在信号的分解期间，降采样引进了畸变，这就需要在分解和重构阶段精心选择关系紧密但不一定一致的滤波器以尽量消除这种畸变。低通分解滤波器(L)和高通分解滤波器(H)以及重构滤波器构成正交镜像滤波器系统。,函数 逆离散小波变换,function y=myidwt(cA,cD,lpr,hpr)%逆离散小波变换%lpr、hpr：重构所用的低通、高通滤波器lca=length(cA);%求出平均、细节部分分解系数的长度lcd=length(cD);while(lcd)=(lca)%若lcd小于lca，则重构停止upl=upspl(cA);%对平均部分系数进行升抽样cvl=conv(upl,lpr);%低通卷积cD_up=cD(lcd-lca+1:lcd);%取出本层重构所需的细节部分系数uph=upspl(cD_up);%对细节部分系数进行升抽样cvh=conv(uph,hpr);%高通卷积cA=cvl+cvh;%用本层重构的序列更新cAcD=cD(1:lcd-lca);%舍弃本层重构用到的细节部分系数lca=length(cA);lcd=length(cD);end%lcd lca，重构完成，结束循环y=cA;%y 等于平均部分系数序列 cA,函数 升采样,function y=upspl(x)%对输入序列x进行升抽样，即对序列x每个元素之间插零，例如 x=x1,x2,x3,x4,上抽样后为 y=x1,0,x2,0,x3,0,x4;N=length(x);%读取输入序列长度M=2*N-1;%输出序列的长度是输入序列长度的2倍再减一for i=1:M%输出序列的偶数位为0，奇数位按次序等于相应位置的输入序列元素 if mod(i,2)y(i)=x(i+1)/2);else y(i)=0;endend,MATLAB函数,升采样函数dyadup：Y=dyadup(X,EVENODD)Y=dyadup(X)Y=dyadup(X,EVENODD,type)Y=dyadup(X,type,EVENODD)重构函数waverec：X=waverec(C,L,wname)X=waverec(C,L,Lo_R,Hi_R),例 4,原始信号clear all;t=0:pi/100:4*pi;s=sin(t+pi/4);subplot(211);plot(s);title(原始信号);,例 4,小波分解c,l=wavedec(s,3,db1);%多层小波分解ca3=appcoef(c,l,db1,3);%提取低频系数（第三层近似系数）cd3=detcoef(c,l,3);%提取第3层高频（近似）系数cd2=detcoef(c,l,2);%提取第2层高频（近似）系数cd1=detcoef(c,l,1);%提取第1层高频（近似）系数figure(2);subplot(221);plot(ca3);title(第三层低频系数);subplot(222);plot(cd1);title(第一层高频系数);subplot(223);plot(cd2);title(第二层高频系数);subplot(224);plot(cd3);title(第三层高频系数);,例 4,小波重构s1=waverec(c,l,db1);figure(1);subplot(2，1，2);plot(s1);grid;title(db1小波重构信号);,重构效果非常好！,11.7 小波分析工具箱,常用的小波基函数,11.7 小波分析工具箱,计算小波滤波器系数的函数,11.7 小波分析工具箱,用于验证算法的数据文件,例 5：小波滤波器,clear all;1.正弦波f1=50;%频率1f2=100;%频率2fs=2*(f1+f2);%采样频率%fs=256Ts=1/fs;%采样间隔N=120;%采样点数n=1:N;y=sin(2*pi*f1*n*Ts)+sin(2*pi*f2*n*Ts);%正弦波混合figure(1)subplot(2,1,1);plot(y);title(两个叠加的正弦信号);subplot(2,1,2);plot(abs(fft(y);title(频谱);,例 5：小波滤波器,2.小波滤波器谱h=wfilters(db30,l);%低通g=wfilters(db30,h);%高通h=h,zeros(1,N-length(h);%补零g=g,zeros(1,N-length(g);figure(2);subplot(2,2,1);plot(g);title(低通h(n);subplot(2,2,2);plot(h);title(低通g(n);subplot(2,2,3);plot(abs(fft(h);title(低通频响H(w);subplot(2,2,4);plot(abs(fft(g);title(高通频响G(w);,例 5：小波滤波器,3.MALLAT分解算法sig1=ifft(fft(y).*fft(h);%低通(低频分量)sig2=ifft(fft(y).*fft(g);%高通(高频分量)figure(3);%信号图subplot(2,2,1);plot(real(sig1);title(分解信号1：低通);subplot(2,2,2);plot(real(sig2);title(分解信号2：高通);subplot(2,2,3);plot(abs(fft(sig1);title(分解信号1频谱);subplot(2,2,4);plot(abs(fft(sig2);title(分解信号2频谱);,例 5：小波滤波器,4.MALLAT重构算法sig1=dyaddown(sig1);%2抽取：降采样sig2=dyaddown(sig2);%2抽取：降采样sig1=dyadup(sig1);%2插值：升采样sig2=dyadup(sig2);%2插值：升采样sig1=sig1(1,1:N);%去掉最后一个零sig2=sig2(1,1:N);%去掉最后一个零hr=h(end:-1:1);%重构低通gr=g(end:-1:1);%重构高通hr=circshift(hr,1);%位置调整圆周右移一位gr=circshift(gr,1);%位置调整圆周右移一位sig1=ifft(fft(hr).*fft(sig1);%低频sig2=ifft(fft(gr).*fft(sig2);%高频sig=sig1+sig2;%源信号,例 5：小波滤波器,5.结果figure(4);subplot(2,2,1);plot(real(sig1);title(重构低频信号);subplot(2,2,2);plot(real(sig2);title(重构高频信号);subplot(2,2,3);plot(abs(fft(sig1);title(重构低频信号频谱);subplot(2,2,4);plot(abs(fft(sig2);title(重构高频信号频谱);,11.8 消噪处理,小波消噪处理方法一般有三种。强制消噪处理。把小波分解结构中的高频系数全部变为0，即把高频部分全部滤除掉，然后冉对信号进行重构处理。这种方法比较简单，重构后的消噪信号也比较平滑，但容易丢失信号的有用成份。默认阈值消噪处理。该方法利用ddencmp函数产生信号的默认阈值，然后利用wdencmp函数进行消噪处理。给定软、硬阈值消噪处理。在实际的消噪处理过程中，阈值往往可以通过经验公式获得。而且这种阈值比默认阈值更具有可信度。在进行阈值量化处理中可用wthresh函数进行。,基于小波的信号消噪处理和压缩处理,11.8 消噪处理,例 6：小波去噪,clear all;snr=4;%设置信噪比init=2055415866;%设置随机数xref,x=wnoise(1,11,snr,init);%产生原始信号xref=xref(1:2000);x=x(1:2000);用sym8小波进行五层分解，实现去噪xd=wden(x,heursure,s,one,5,sym8);subplot(311);plot(xref);title(原始信号);subplot(312);plot(x);title(含噪信号);subplot(313);plot(xd);title(去噪信号);,二维小波变换是由一维小波变换扩展而来的，二维尺度函数和二维小波函数可由一维尺度函数和小波函数张量积得到，即：,补充 2：二维小波变换,图像是二维信号，其小波变换相当于二次一维信号的小波变换：（1）第一次一维信号的小波变换：相当于图像的行变换（2）第二次一维信号的小波变换：相当于图像的列变换,一、基本概念,【例】图像的二维小波变换包括沿行向(水平方向)和列向(垂直方向)滤波和2-降采样。,补充 2：二维小波变换,对图像沿行进行变换，得到系数矩阵IL(x,y)和IH(x,y)；再对IL(x,y)和IH(x,y)分别沿列变换；最后得到一层小波分解的4个子图：ILL(x,y)I(x,y)的近似子图 IHL(x,y)I(x,y)的水平方向细节子图 ILH(x,y)I(x,y)的垂直方向细节子图 IHH(x,y)I(x,y)的对角线方向细节子图,对近似子图重复此过程，直到确定的分解水平。,(a)一层分解，(b)二层分解,补充 2：二维小波变换,一维行小波变换,一维列小波变换,补充 2：二维小波变换,二、二维Mallat算法,（1）基本概念,补充 2：二维小波变换,其中：,可得：,可分离的二维小波变换,补充 2：二维小波变换,【例】可分离的二维小波变换,（1）确定LL,（2）另一种计算方法,（2）确定HL,补充 2：二维小波变换,补充 2：二维小波变换,简单的压缩方案:,方案1:只保留低频部分.,方案2:全局阈值 法.,方案3:保留绝对值较大的若干小波系数,二维小波变换的塔式结构图,图像块的三级小波分解系数,【例子】,（3）二维重构,补充 2：二维小波变换,补充 2：二维小波变换,一维列小波变换,一维行小波变换,双正交滤波器的情况,例 7：二维图像,T=256;%图像维数SUB_T=T/2;%子图维数1.调原始图像矩阵load wbarb;%下载图像f=X;%原始图像2.二维小波分解l=wfilters(db10,l);%低通分解,长度为20L=T-length(l);l_zeros=l,zeros(1,L);%矩阵行数与输入图像一致，为2的整数幂h=wfilters(db10,h);%高通分解,长度为20h_zeros=h,zeros(1,L);%矩阵行数与输入图像一致，为2的整数幂,例 7：二维图像,for i=1:T;%列变换 row(1:SUB_T,i)=dyaddown(ifft(fft(l_zeros).*fft(f(:,i).;%row(SUB_T+1:T,i)=dyaddown(ifft(fft(h_zeros).*fft(f(:,i).;%end;for j=1:T;%行变换 line(j,1:SUB_T)=dyaddown(ifft(fft(l_zeros).*fft(row(j,:);line(j,SUB_T+1:T)=dyaddown(ifft(fft(h_zeros).*fft(row(j,:);end;decompose_pic=line;%分解矩阵lt_pic=decompose_pic(1:SUB_T,1:SUB_T);%左上方：fi(x)*fi(y)rt_pic=decompose_pic(1:SUB_T,SUB_T+1:T);%右上方：fi(x)*psi(y)lb_pic=decompose_pic(SUB_T+1:T,1:SUB_T);%左下方：psi(x)*fi(y)rb_pic=decompose_pic(SUB_T+1:T,SUB_T+1:T);%右下方：psi(x)*psi(y),例 7：二维图像,3.分解结果显示figure(1);colormap(map);subplot(2,1,1);image(f);title(原始图像);subplot(2,1,2);image(abs(decompose_pic);%title(分解图像);figure(2);colormap(map);subplot(2,2,1);image(abs(lt_pic);title(Phi(x)*Phi(y);subplot(2,2,2);image(abs(rt_pic);title(Phi(x)*Psi(y);subplot(2,2,3);image(abs(lb_pic);title(Psi(x)*Phi(y);subplot(2,2,4);image(abs(rb_pic);title(Psi(x)*Psi(y);,例 7：二维图像,4.重构图像l_re=l_zeros(end:-1:1);%重构低通滤波l_r=circshift(l_re,1);%位置调整h_re=h_zeros(end:-1:1);%重构高通滤波h_r=circshift(h_re,1);%位置调整 top_pic=lt_pic,rt_pic;%图像上半部分t=0;for i=1:T;%行插值低频 if(mod(i,2)=0)topll(i,:)=top_pic(t,:);%偶数行保持 else t=t+1;topll(i,:)=zeros(1,T);%奇数行为零 endend;for i=1:T;%列变换 topcl_re(:,i)=ifft(fft(l_r).*fft(topll(:,i);%圆周卷积FFTend;,例 7：二维图像,bottom_pic=lb_pic,rb_pic;%图像下半部分t=0;for i=1:T;%行插值高频 if(mod(i,2)=0)bottomlh(i,:)=bottom_pic(t,:);%偶数行保持 else bottomlh(i,:)=zeros(1,T);%奇数行为零 t=t+1;endend;for i=1:T;%列变换 bottomch_re(:,i)=ifft(fft(h_r).*fft(bottomlh(:,i);end;construct1=bottomch_re+topcl_re;%列变换重构完毕,例 7：二维图像,left_pic=construct1(:,1:SUB_T);%图像左半部分t=0;for i=1:T;%列插值低频 if(mod(i,2)=0)leftll(:,i)=left_pic(:,t);%偶数列保持 else t=t+1;leftll(:,i)=zeros(T,1);%奇数列为零 endend;for i=1:T;%行变换 leftcl_re(i,:)=ifft(fft(l_r).*fft(leftll(i,:);%圆周卷积FFTend;,例 7：二维图像,right_pic=construct1(:,SUB_T+1:T);%图像右半部分t=0;for i=1:T;%列插值高频 if(mod(i,2)=0)rightlh(:,i)=right_pic(:,t);%偶数列保持 else rightlh(:,i)=zeros(T,1);%奇数列为零 t=t+1;endend;for i=1:T;%行变换 rightch_re(i,:)=ifft(fft(h_r).*fft(rightlh(i,:);%圆周卷积FFTend;construct_pic=rightch_re+leftcl_re;%重建全部图像,例