《基于MATLAB的小波分析應(yīng)用》課件第10章

第10章小波變換與矩陣

方程求解10.1小波變換快速求解矩陣方程

10.2快速求解矩陣方程實(shí)例分析10.1小波變換快速求解矩陣方程采用積分方程法后，偏微分方程可以轉(zhuǎn)化為下列矩陣方程：ZJ=V

對(duì)于一般迭代的解法，例如共軛梯度法，需要k·O(N?2)計(jì)算復(fù)雜度，其中，k是迭代需要的步數(shù)。(10-1)用小波變換求解(10-1)矩陣方程的步驟如下：

(1)構(gòu)造正交小波變換矩陣W(10-2)(10-3)式中，H和G分別是由消失矩為P的分解低通濾波器h和分解高通濾波器g構(gòu)成的矩陣。它們的每行是由長(zhǎng)度為N/2n-1的向量[h(0)，h(1)，…h(huán)(2P-2)，h(2P-1)，0，0，…，0]和向量[g(0)，g(1)，…g(2P-2)，g(2P-1)，0，0，…0]分別圓周2移位得到。根據(jù)濾波器的正交性質(zhì)，可知WWT=WTW=I。

(2)由于U=WZ，對(duì)矩陣Z和向量V分別作二維和一維小波變換(10-4)(10-5)經(jīng)變換后，矩陣方程(10-1)等價(jià)于(10-6)式中，。稱(chēng)式(10-6)為小波域矩陣方程。

(3)通過(guò)硬閾值方法，稀疏化矩陣。(10-7)式中，閾值的選擇基于L1范數(shù)標(biāo)準(zhǔn)(10-8)小波域矩陣經(jīng)稀疏化后轉(zhuǎn)化為(10-9)

(4)采用共軛梯度(CG)或者廣義最小余量法(GMRES)求解小波域的解。

(5)通過(guò)逆小波變換，得到矩陣方程的近似解。(10-10)對(duì)于利用小波變化加速矩陣方程的求解，需注意以下幾點(diǎn)：

(1)若式(10-2)中，時(shí)，該變換等價(jià)于一維小波變換；而當(dāng)n≥2時(shí)，對(duì)矩陣的變換不是傳統(tǒng)的二維小波變換，而是一種被稱(chēng)為類(lèi)小波的變換。該變換介于小波變換和小波包變換之間，其高頻分辨率高于小波變換但低于小波包變換，計(jì)算時(shí)間慢于小波變換但快于小波包變換。

(2)根據(jù)矩陣方程求解的步驟可知，小波變換耗時(shí)主要取決于式(10-4)的矩陣變換的耗時(shí)，該步計(jì)算的復(fù)雜度為O(4pkN?2)。因?yàn)楫?dāng)矩陣的維數(shù)很大時(shí)，k是一個(gè)很大的量(k同矩陣的維數(shù)有正增長(zhǎng)關(guān)系)，而小波變換的稀疏度(非零元素的比例)p為很小的量(1%左右)，因此有，這表明矩陣求解時(shí)間是大于小波變換時(shí)間的。對(duì)于很多物理問(wèn)題，經(jīng)小波變換預(yù)處理后，總體可以節(jié)約的時(shí)間大于90%以上。

(3)以上分析針對(duì)的是矩陣W的正交小波的情況，對(duì)于雙正交小波來(lái)說(shuō)，同樣可以構(gòu)造對(duì)偶的小波矩陣，滿(mǎn)足雙正交關(guān)系：

(10-11)這種正交或雙正交關(guān)系使得變換后的矩陣條件數(shù)和變換前矩陣Z的條件數(shù)一致，不會(huì)造成矩陣的病態(tài)。

(4)上述算法能節(jié)省時(shí)間，是由于原始矩陣Z經(jīng)小波變換和閾值化后得到的矩陣，可以達(dá)到很低的稀疏率。對(duì)于某些實(shí)際的物理問(wèn)題，原始矩陣Z的元素分布類(lèi)似于一種無(wú)規(guī)律振蕩狀態(tài)(如白噪聲矩陣)，因此，經(jīng)過(guò)小波變換后的矩陣并不能達(dá)到滿(mǎn)意的稀疏率。如果強(qiáng)行增大閾值來(lái)改善稀疏化程度，會(huì)使最終的求解與原方程解之間的誤差增大。10.2快速求解矩陣方程實(shí)例分析下面結(jié)合具體實(shí)例，給出矩陣求解操作的過(guò)程與MATLAB程序。

1.構(gòu)造矩陣W的程序下面利用MATLAB編程，給出如何利用小波變換構(gòu)造矩陣W。

MATLAB程序如下：

%構(gòu)造小波矩陣

clear;

clc;[h,g]=wfilters('db7','d'); %分解低通和高通濾波器N=512; %矩陣維數(shù)(2的整數(shù)次冪)L=length(h); %濾波器長(zhǎng)度rank_max=log2(N);

%最大層數(shù)rank_min=double(int8(log2(L)))+1;%最小層數(shù)ww=1; %預(yù)處理矩陣%構(gòu)造矩陣Wforjj=rank_min:rank_maxnn=2^jj;%構(gòu)造向量

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%構(gòu)造正交矩陣

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))，顯示結(jié)果為ans=5.9958e-012在命令行輸入命令：norm(ww.'*ww-eye(N,N))，顯示結(jié)果為ans=5.9965e-012由此可知，得到的結(jié)果都是10-12量級(jí)的很小的數(shù)，說(shuō)明構(gòu)造的矩陣是正交矩陣。用spy(ww)命令可以看出ww的非零元素分布，如圖10.1所示。這說(shuō)明小波矩陣是很稀疏的，因此小波變換的時(shí)間是很快的。圖10.1小波矩陣非零元素分布

2.矩陣求解的程序下面利用MATLAB編程，給出如何利用小波變換進(jìn)行矩陣求解。

MATLAB程序如下：

%利用小波變換求解矩陣方程

clear;

clc;

loadww;

N=512;%構(gòu)造測(cè)試矩陣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經(jīng)小波變換并壓縮后的矩陣的非零點(diǎn)元素分布%反小波變換result_app=ww.'*j_t;%計(jì)算相對(duì)誤差e=norm(result_app-result)/norm(result)%繪制稀疏度圖spy(z_t)程序執(zhí)行結(jié)果為e=0.0591非零元素分布如圖10.2所示。圖10.2經(jīng)小波變換并壓縮后的矩陣的非零點(diǎn)元素分布從上圖可以看出，經(jīng)過(guò)小波變換得到的矩陣方程的解和原方程解的相對(duì)誤差小于5%，且稀疏度為2770/(512*512)=1.06%。將上面程序中的z改成隨機(jī)矩陣randn(N