mallat演算法

❶ Mallat演算法

6.9.1.1 尺度空間的有限分解

MRA框架表明f（t）∈L²（R）可分解為無窮個小波分量的直和，但在實際應用中，僅知道f（t）的近似函數。為不失一般性，可假設原信號是在a=1或j=0的解析度下測得的，用f⁰（t）表示，它屬於子空間V₀。而子空間V₀又可分解成兩個子空間。因此，在MRA框架下理解為f⁰（t）∈V₀，這樣就有如下的尺度空間的有限分解表現

V₀=V₁⊕W₁=（V₁⊕W₁）⊕W₀=（V₂⊕W₂）⊕W₁

=（V_j⊕W_j）⊕W_j-1⊕W_j-2⊕…⊕W₁

︙

=V_J+W_J⊕W_J-1⊕W_J-2⊕…⊕W₁

其中子空間及分量分別為

地球物理信息處理基礎

V₀的有限分解關系僅是MRA無窮分解中的一部分。因此，就子空間而言、就函數分量而言以及就頻率范圍而言，有限分解的含義都與MRA相同。例如，V₀中的元素f⁰（t）是有限頻率范圍的，f^J（t）∈V_J是f⁰（t）的最低頻表現，W_j中元素δ^j（t）是具有特定帶寬的，它們互不重疊，這些頻帶的總和就是f⁰（t）的頻率范圍。

為了數字計算和分析處理的目的，需要將f^j（t）和δ^j（t）用離散數據來表示。顯然，｛c_j，k｝∈l²是合適的，因為

表明f^j（t）∈V_j和｛c_j，k｝∈l²是一一對應的，或者說利用｛c_j，k｝就一定可恢復f^j（t）。根據同樣的理由，可用數據｛d_j，k｝∈l²來表示δ^j（t）∈W_j，這類數據｛c_j，k｝和｛d_j，k｝對應著特定的離散表現，j表示尺度方面的2進制規則離散，k表示整節點標記的時移離散。

6.9.1.2 分解演算法

分解演算法要實現的目標是：在｛φ（t-k）｝是標准正交基條件下，已知｛c_j-1，k｝、｛h_k｝和｛g_k｝，求出｛c_j，k｝和｛d_j，k｝。

ψ（t）關於φ（t）的兩尺度關系式（6-95）提供了一條由尺度函數φ（t）構造母小波ψ（t）的途徑。根據尺度函數

地球物理信息處理基礎

和小波函數

地球物理信息處理基礎

得到

地球物理信息處理基礎

用l取代2k+l，上式成為

地球物理信息處理基礎

計算f（t）∈L²（R）與上式兩端的內積，便得到如下計算小波級數系數的一個公式

地球物理信息處理基礎

式中c_j-1，l=＜f^j，φ_j-1，l＞=＜f，φ_j-1，l＞是信號f（t）與φ_j-1，l（t）的內積，即c_j，l=＜f，φ_j，l＞。信號在V_j的正交投影

，將兩尺度關系式的φ_j，k（t）代入定義式φ_j，k（t）=2^-j/2φ（2^-jt-k），有

地球物理信息處理基礎

由此可以計算f（t）與φ_j，k（t）的內積，即

地球物理信息處理基礎

在數學上，f^j（t）是平方可積連續函數，它是L²（R）的元素；c_j，l（l∈Z）是平方可和序列，它是平方可和序列矢量空間l²（Z）中的元素。

說明f^j（t）與c_j，k（k∈Z）有等效對應關系，在數學上稱l²（Z）與L²（R）之間存在著同構關系。顯然，c_j，k比f^j（t）更適合於在計算機上運算。c_j，k稱為在分辨2^j下的f（t）離散逼近，因此，可以把f⁰（t_k）=f（t_k）的采樣值作為｛c_j，k｝的初始數據，即｛c_0，k｝=｛f（t_k）｝。

綜合式（6-109）和（6-110），得到一般的分解公式

地球物理信息處理基礎

式（6-111）的實現過程如圖6-26所示，每個尺度j所存儲的數據｛c_j，k｝都是按整數編號的。以j尺度層為基礎來觀察j-1尺度層，j尺度層的采樣節點編號k對應著j-1尺度層上編號為2k的采樣節點，或者說。j-1尺度的采樣節點是在j尺度采樣節點基礎上均勻加密的結果；若以j-1尺度層為基礎來觀察j尺度層，則j-1尺度層上隔2取樣（「隔1取1」）的節點正好對應著j尺度層上的采樣節點。圖6-26表明了由j-1尺度向j尺度的變換過程，｛h_k｝可看做濾波器的單位沖激響應（權系數）。假設｛h_k｝僅有6個元素，則式（6-111）所表明的變換過程相當於把｛h_-2，h_-1，h₀，h₁，h₂，h₃，｝作為權值，其中心點h₀對准｛c_j，k｝後再作加權平均，即

地球物理信息處理基礎

這一實現過程是非常快捷的，相當於先計算序列c_j-1，k（k∈Z）與

的卷積，然後抽取偶數下標（隔1抽1）的卷積結果，即得到逼近級數的系數c_j，k（k∈Z）。

圖6-26 分解演算法計算｛c_j，k｝示意圖

圖6-26雖然僅表明由｛c_j-1，k｝計算｛c_j，k｝的變換過程，它同樣也表明了由｛c_j-1，k｝計算｛d_j，k｝的變換過程，只不過要將｛h_k｝換成｛g_k｝而已。對正交小波而言，g_k=（-1）^kh_-k+1是由｛h_k｝決定的。

由c_0，k開始，利用式（6-111）進行迭代運算，陸續計算出c_1，k、c_2，k等等，與此同時，利用c_0，k、c_1，k、c_2，k等值，同樣不斷計算出d_1，k、d_2，k等小波級數系數值。

式（6-111）所表明的計算過程由運算元表示會更簡單些。記A^j=｛c_j，k｝，D^j=｛d_j，k｝。記運算元H：l²➝l²，其運算意義如式（6-111）的第1式所示，即

地球物理信息處理基礎

同樣，記運算元G：l²➝l²，其運算意義如式（6-111）的第2式所示，即

地球物理信息處理基礎

採用運算元表示後，式（6-111）所表明的分解演算法結構如圖6-27所示。它表明

A^j=H^jA⁰；D^j=GH^j-1A⁰

即A⁰經H運算元j次作用後即可獲得A^j，A⁰經H運算元j-1次作用後再經G運算元作用1次即可獲得D^j。

圖6-27 分解演算法結構示意圖

圖6-27還表明，只要在細密采樣間隔的尺度層次上給定A⁰，就可利用分解演算法快速地獲得較粗采樣間隔尺度層上的有關數據A^j和D^j（0＜j≤J）。假設實際問題中A⁰有N個數據，則2尺度層上的A¹和D¹各有N/2個數據，依此類推；還假設｛h_k｝和｛g_k｝分別有M個數據；那麼，用A⁰計算A¹和D¹共需2MN/2次運算，從2到3尺度層共需2MN/4次運算，依此類推可知，要得到A^j和D^j（0＜j≤J）共需2MN（2^-1+2^-2+…2^-J）次運算。由此可見，分解演算法是快速的。

原信號f⁰（t）具有下列唯一分解

地球物理信息處理基礎

δ_j是信號在W_j（j=1，2，…，J）各子空間上的正交投影，它們是從一個較精細的逼近變成較粗略的逼近（兩個逼近的解析度相鄰近）時所丟失的信息，即

地球物理信息處理基礎

可將這一不斷降低逼近解析度的過程看成是「一層又一層地把信號進行剝皮」的過程。當J選得足夠大時，「剝」下來的信息總和

地球物理信息處理基礎

足夠多，將足以精確表示原信號f⁰（t），而最終的逼近信號f^J（t）的解析度已經非常低，這樣反而可以把式（6-113）當成原信號的估計，而把f^J（t）看成是估計誤差。這就是說，用小波級數在所有解析度下的全部系數（j=1，2，…，J）來代替原信號，其誤差f^J（t）可以任意小。按照這種解釋，式（6-111）演算法就是將f⁰（t）的信息（c_0，k是它的離散表示）表示成c_1，k～c_J，k等信息和一個估計誤差（實際上它是在解析度最低即2^J下的逼近）c_J，k。這一過程實際上是在一次又一次地改變著正交基（或子空間）。

6.9.1.3 信號重構演算法

重構演算法是分解演算法的逆過程。此時已知數據｛c_j，k｝和｛d_j，k｝（0≤j≤J），希望利用這些數據快速准確地重構出原始數據。

一般而言，相鄰兩解析度下的逼近信號存在著下列關系：

地球物理信息處理基礎

計算上式左右兩端與φ_j-l，k（t）的內積，得

地球物理信息處理基礎

為了便於討論運算過程，在上式中令m=p-2k，所以信號重構公式為

地球物理信息處理基礎

根據式（6-115）計算｛

｝的過程見圖6-28所示。為了理解式（6-115）和圖6-28所示的計算過程，應注意下面兩點：一是j-1尺度層偶數編號的采樣點對應著j尺度層上的采樣點，二是示意圖中假設｛h_k｝的有限數據為｛h_-2，h_-1，h₀，h₁，h₂，h₃，｝。

可以看出，計算｛

｝可分為兩步進行：

地球物理信息處理基礎

第一步計算j-1尺度層上的偶數編號采樣點處的｛

｝，此時可對j尺度層采樣節點循序推進，所用數據為｛h_-2，h₀，h₂，｝，按公式

地球物理信息處理基礎

計算出j-l尺度層上偶數編號節點處的相應值。

第二步計算j-1尺度層上的奇數編號的采樣點處的｛

｝，其循序推進過程仍同於第一步，但需另用｛h₃，h₁，h_-1｝，且按公式

地球物理信息處理基礎

圖6-28也表明了式（6-115）中｛d_j，k｝計算｛

｝的過程，此時要採用｛g_k｝，計算式應為

地球物理信息處理基礎

為了更清楚地表明｛c_j，p｝和｛d_j，p｝，0≤ j≤J-1重構｛c_0，p｝的關系，與前述的分解情況一樣，記A^j=｛c_j，p｝，D^j=｛d_j，p｝，採用運算元的記號H^*：l²➝l²，其運算意義為

地球物理信息處理基礎

記G^*：l²➝l²，其運算意義為

地球物理信息處理基礎

在這些重構運算元意義下，式（6-115）可表示為

地球物理信息處理基礎

於是從尺度層j=J到尺度層j=0的重構演算法過程可用圖6-29表示。重構演算法也是快速的，現估計其計算量。設｛h_m｝和｛g_m｝分別有M個數據，0尺度層｛c_j，0｝有N個數據，則尺度層1上｛c_j，1｝有N/2個數據，用｛c_1，p｝計算｛

｝的偶節點值和奇節點值分別需要（M/2）N/2次運算；由｛c_1，p｝計算｛

｝共需MN/2次運算；｛d_1，p｝計算｛

｝也需MN/2次運算；由1尺度層數據重構0尺度層數據共需MN次運算。可依此類推出其它相鄰尺度層之間採用重構演算法的運算量，最後，得總運算量為

2MN（2^-1+2^-2+…2^-J）

圖6-29 信號重構演算法示意圖

圖6-30所示的是信號分解與重構演算法的計算流程示意圖，圖中的「

」表示「按因子2抽取」或「隔1抽1」運算。符號「

」表示，「按因子2內插」或「兩相鄰樣本間內插一個零樣本」。輸入是c_0，k，相當於原信號的離散表示；輸出是不同尺度下的小波級數系數c_j，k（j∈Z）以及信號在子空間V_J的正交投影c_J，k。這是在做信號分析，其分解公式見（6-111）。與此相反，當輸入是d_j，k（j∈Z）和c_J，k，輸出是c_0，k時，為信號重建或合成，此時有

地球物理信息處理基礎

將合成與分析演算法合起來畫成圖6-30所示的形式，稱之為Mallat演算法。

圖6-30 信號分解與重構演算法示意圖

6.9.1.4 Mallat演算法實現中的一些問題

Mallat演算法是一種純數字的快速遞推演算法，在使用Mallat演算法時，有一些具體問題需引起注意。

（1）對正交尺度函數φ（t）而言，Mallat演算法中僅需數據｛c_0，k｝和｛h_k｝可進行快速的分解和重構遞推運算。要存儲的數據為｛c_J，k｝和｛d_j，k｝，0＜j≤J，這些有用的數據的存儲量等於｛c_0，k｝的數據存儲量。特別值得強調的是，Mallat演算法中隱含著兩類關系，一類是關於多分辨分析方面的，例如對0＜j≤J，有

地球物理信息處理基礎

f^j-1（t）=f^j（t）+δ^j（t）

地球物理信息處理基礎

另一類關系是由於尺度函數φ（t）平移正交性產生的，例如

c_j，k=＜f，φ_j，k＞；d_j，k=＜δ^j，ψ_j，k＞

g_k=（-1）^kh_-k+1，0＜j≤J

Mallat演算法正是利用了這些關系，在演算法實施過程中不需尺度函數φ（t）和小波函數ψ（t）的具體形式，只要求它們存在並找出｛h_k｝，就可以順利地進行分解和重構處理了。因此，只要查得正交尺度函數雙尺度方程的傳遞系數｛h_k｝，就可以應用Mallat演算法了。

順便指出，如果尺度函數φ（t）（例如樣條函數）不是平移正交的，它雖然可以生成MRA，但由此構造出的樣條小波ψ（t）僅關於尺度正交，沒有平移正交性。此時

地球物理信息處理基礎

所以Mallat演算法不再適用，必須另行推導相應的分解和重構演算法。

（2）初始數據｛c_0，k｝的選用，在正交小波分解中，φ（t）是正交尺度函數，0尺度層上的展180開系數｛c_0，k｝=＜f（t），φ_0，k（t）＞，用復雜的計算來確定初始數據｛c_0，k｝是不合算的，應採用變通處理辦法，即簡單採用｛c₀，k｝=｛f（t_k）｝，用細尺度層上的采樣值作為初始數據｛c_0，k｝。這種做法似乎有些不嚴密，但可以證明，雖然｛f（t_k）｝作為｛c_0，k｝的近似值時有微小誤差，但數據｛f（t_k）｝同樣有效地表現了f⁰（t）的變化波動狀況和有效的頻率范圍，這種替代不會影響對f⁰（t）的時頻分析；同時還應看到，用｛f（t_k）｝作初始數值，不僅可使問題簡單化，而且也可使Mallat演算法准確地分解和重構初始數據。總之，用｛f（t_k）｝替代｛c_0，k｝是實用方便的。

（3）分解層數和采樣間隔的關系，這個問題主要從以下幾方面考慮便可得出結論。

第一，因為最細的0尺度層的采樣間隔T決定了f⁰（t）的頻率范圍，由取樣效應可知，最大的頻率范圍為｜ω（f⁰）｜≤1/（2T）；同樣，最粗的J尺度層的采樣間隔為J=2^JT，最低的頻率范圍為｜ω（f^J）｜≤1/（2T）=2^-J｜ω（f⁰）｜。於是可從需要分辨的最高頻率和需要分辨的最低頻率這兩個指標來決定最細尺度層的采樣間隔和數據分解的層數。

第二，在最細的0尺度層上，應取用多少個數據才能滿足J個層次的數據分解呢？在Mallat演算法中，｛c_1，k｝的數據量僅為｛c_0，k｝數據量的一半，依此類推。同樣，在J尺度層至少要取用N_J個數據才能表現低頻量，於是推知，在0尺度層，至少取用N=2^JN_J個數據，才能滿足J個分解層次的需要。

（4）在Mallat演算法的運算中，需用到所存儲的數據外面的數據（見圖6-31），圖中實線框內數據是要存儲的。在分解演算法中，若｛h_k｝僅有6個數據，由圖6-26可知，要用到實線框外的數據c_0，-2和c_0，-1才能計算出c_1，-1，要用到c_0，N0+1、c_0，N0+2和c_0，N0+3才能計算出c_1，N1，其它層次的情形類似。由於細密采樣層（圖6-31中對應著j=0的尺度層次）中的近似函數f⁰（t）∈L²（R），當t➝∞時，f⁰（t）➝0，實線框內足夠多的采樣數據｛c_0，k｝，k=0，1，…，N₀，已反應了f⁰（t）的基本特性，f⁰（t）在圖6-31實線框外的數據幾乎消失，因此在實際分解過程中可簡單地令實線框外的數據為零。同樣，重構過程也會用到實線框以外的數據（見圖6-28），也可以簡單地令實線框外的數據為零。

圖6-31 Mallat演算法涉及沒有存儲的數據示意圖

另一種辦法也是常用的。在最細尺度層上較多地取用數據，在計算過程中適當地多存儲些數據，如圖6-31中虛線框所示。此時應以實線框存儲數據作為分解重構演算法和進行數字分析的依據。

6.9.1.5 Mallat 演算法所表現的頻域分解特點

有限尺度空間的正交小波子空間的直和分解關系，例如，

V₀=V₁⊕W₁=V₂⊕W₂⊕W₁=V₃⊕W₃⊕W₂⊕W₁

在Mallat演算法中是通過運算元H和G來表現的。數據A^j表徵f^j（t）∈V_j，數據D^j表徵δ^j（t）∈W_j，那麼在Mallat演算法中，A^j=HA^j-1；D^j=GA^j-1，從而實現了子空間的分解。這種子空間分解和運算元H與G之間的關系示意於圖6-32中（以V₀分解成3層為例）。

圖6-32 Mallat演算法所確定的數據分解和子空間分解的對應關系

因為f^j-1（t）和f^j（t）都是有限頻率范圍的，f^j（t）的頻率范圍僅是f^j-1（t）的相對低頻部分，δ^j（t）的頻率范圍是f^j-1（t）的相對高頻部分，所以δ^j（t）是f^j-1（t）的關於高頻帶的分量。換成子空間來描述，V_j-1所表現的頻率范圍被分解為兩部分，一部分是由V_j表現的低頻部分，另一部分是由W_j表現的頻帶部分。因此，V₀尺度空間的直和分解，表明把f⁰（t）的頻率范圍分解為由W₁、W₂、W₃和V₃所表現的頻帶，且這些頻帶是互不重疊的。這里所說的W₂所表現的頻帶寬度，也就是其基函數ψ_2，k（t）（窗函數）的頻窗寬度。由aD_ω可知，W₂所表現的頻帶寬是W₁的一半，尺度指標增加，則小波子空間所表現的頻帶寬分半減小。由於對任意尺度而言，小波子空間所表現的時頻窗面積是恆定不變的常數，所以對時頻窗的時窗寬度而言，尺度指標增加，則相應的時窗寬度成倍增加。圖6-33表示了正交小波分解用於分析t^*處局部時域信號時，在各個頻帶的時頻窗表現。

圖6-33 正交小波分解在各頻帶的時頻窗示意圖

❷ 什麼是Mallat演算法

小波變換的多解析度分析(或多尺度分析)是建立在函數概念上的理論，多解析度分析（Multiresolution Analysis，MRA）概念是由 S.Mallat和 Y.Meyer 在前人大量工作的基礎上於1986年提出的，從空間的概念上形象的說明了小波的多解析度特性，隨著尺度由大到小變化，在各尺度上可以由粗到細的觀察圖像的不同特徵。在大尺度時，觀察到圖像的輪廓，在小尺度的空間里，則可以觀察圖像的細節。
1989年，Mallat在小波變換多解析度分析理論與圖像處理的應用研究中受到塔式演算法的啟發，提出了信號的塔式多解析度分析與重構的快速演算法稱為馬拉特（Mallat）演算法。