《基于MATLAB的小波分析应用》课件第10章

第10章小波变换与矩阵

方程求解10.1小波变换快速求解矩阵方程

10.2快速求解矩阵方程实例分析10.1小波变换快速求解矩阵方程采用积分方程法后，偏微分方程可以转化为下列矩阵方程：ZJ=V

对于一般迭代的解法，例如共轭梯度法，需要k·O(N 2)计算复杂度，其中，k是迭代需要的步数。(10-1)用小波变换求解(10-1)矩阵方程的步骤如下：

(1)构造正交小波变换矩阵W(10-2)(10-3)式中，H和G分别是由消失矩为P的分解低通滤波器h和分解高通滤波器g构成的矩阵。它们的每行是由长度为N/2n-1的向量[h(0)，h(1)，…h(2P-2)，h(2P-1)，0，0，…，0]和向量[g(0)，g(1)，…g(2P-2)，g(2P-1)，0，0，…0]分别圆周2移位得到。根据滤波器的正交性质，可知WWT=WTW=I。

(2)由于U=WZ，对矩阵Z和向量V分别作二维和一维小波变换(10-4)(10-5)经变换后，矩阵方程(10-1)等价于(10-6)式中，。称式(10-6)为小波域矩阵方程。

(3)通过硬阈值方法，稀疏化矩阵。(10-7)式中，阈值的选择基于L1范数标准(10-8)小波域矩阵经稀疏化后转化为(10-9)

(4)采用共轭梯度(CG)或者广义最小余量法(GMRES)求解小波域的解。

(5)通过逆小波变换，得到矩阵方程的近似解。(10-10)对于利用小波变化加速矩阵方程的求解，需注意以下几点：

(1)若式(10-2)中，时，该变换等价于一维小波变换；而当n≥2时，对矩阵的变换不是传统的二维小波变换，而是一种被称为类小波的变换。该变换介于小波变换和小波包变换之间，其高频分辨率高于小波变换但低于小波包变换，计算时间慢于小波变换但快于小波包变换。

(2)根据矩阵方程求解的步骤可知，小波变换耗时主要取决于式(10-4)的矩阵变换的耗时，该步计算的复杂度为O(4pkN 2)。因为当矩阵的维数很大时，k是一个很大的量(k同矩阵的维数有正增长关系)，而小波变换的稀疏度(非零元素的比例)p为很小的量(1%左右)，因此有，这表明矩阵求解时间是大于小波变换时间的。对于很多物理问题，经小波变换预处理后，总体可以节约的时间大于90%以上。

(3)以上分析针对的是矩阵W的正交小波的情况，对于双正交小波来说，同样可以构造对偶的小波矩阵，满足双正交关系：

(10-11)这种正交或双正交关系使得变换后的矩阵条件数和变换前矩阵Z的条件数一致，不会造成矩阵的病态。

(4)上述算法能节省时间，是由于原始矩阵Z经小波变换和阈值化后得到的矩阵，可以达到很低的稀疏率。对于某些实际的物理问题，原始矩阵Z的元素分布类似于一种无规律振荡状态(如白噪声矩阵)，因此，经过小波变换后的矩阵并不能达到满意的稀疏率。如果强行增大阈值来改善稀疏化程度，会使最终的求解与原方程解之间的误差增大。10.2快速求解矩阵方程实例分析下面结合具体实例，给出矩阵求解操作的过程与MATLAB程序。

1.构造矩阵W的程序下面利用MATLAB编程，给出如何利用小波变换构造矩阵W。

MATLAB程序如下：

%构造小波矩阵

clear;

clc;[h,g]=wfilters('db7','d'); %分解低通和高通滤波器N=512; %矩阵维数(2的整数次幂)L=length(h); %滤波器长度rank_max=log2(N);

%最大层数rank_min=double(int8(log2(L)))+1;%最小层数ww=1; %预处理矩阵%构造矩阵Wforjj=rank_min:rank_maxnn=2^jj;%构造向量

p1_0=sparse([h,zeros(1,nn-L)]);p2_0=sparse([g,zeros(1,nn-L)]);%向量圆周移位

forii=1:nn/2p1(ii,:)=circshift(p1_0',2*(ii-1))';p2(ii,:)=circshift(p2_0',2*(ii-1))';end%构造正交矩阵

w1=[p1;p2];mm=2^rank_max-length(w1);w=sparse([w1,zeros(length(w1),mm);zeros(mm,length(w1)),eye(mm,mm)]);ww=ww*w;clearp1;图10.1小波矩阵非零元素分布

clearp2;endsaveww;在命令行输入命令：norm(ww*ww.'-eye(N,N))，显示结果为ans=5.9958e-012在命令行输入命令：norm(ww.'*ww-eye(N,N))，显示结果为ans=5.9965e-012由此可知，得到的结果都是10-12量级的很小的数，说明构造的矩阵是正交矩阵。用spy(ww)命令可以看出ww的非零元素分布，如图10.1所示。这说明小波矩阵是很稀疏的，因此小波变换的时间是很快的。图10.1小波矩阵非零元素分布

2.矩阵求解的程序下面利用MATLAB编程，给出如何利用小波变换进行矩阵求解。

MATLAB程序如下：

%利用小波变换求解矩阵方程

clear;

clc;

loadww;

N=512;%构造测试矩阵fori=1:Nforj=1:Nz(i,j)=1/(abs(i-j)+1);endendv=eye(N,1);result=z\v;%小波变换(稀疏矩阵乘法)z_t=ww*sparse(z)*ww.';v_t=ww*sparse(v);%稀疏化矩阵threshold=3*10^(-3);z_t=z_t.*sparse(abs(z_t)>(threshold*norm(abs(z_t),1)));%求解小波域矩阵方程j_t=z_t\v_t;图10.2经小波变换并压缩后的矩阵的非零点元素分布%反小波变换result_app=ww.'*j_t;%计算相对误差e=norm(result_app-result)/norm(result)%绘制稀疏度图spy(z_t)程序执行结果为e=0.0591非零元素分布如图10.2所示。图10.2经小波变换并压缩后的矩阵的非零点元素分布从上图可以看出，经过小波变换得到的矩阵方程的解和原方程解的相对误差小于5%，且稀疏度为2770/(512*512)=1.06%。将上面程序中的z改成随机矩阵randn(N