第 7 章・小波與其他影像轉換

影像處理中階26 分鐘2026年10月4日

先備知識: 第 6 章・色彩影像處理

你將學到

  • 為什麼所有線性影像轉換(傅立葉、餘弦、沃爾什–哈達瑪、哈爾、小波)其實是同一件事:用另一組積木來描述影像,這組積木叫做基底(basis)。
  • 如何把轉換寫成兩次矩陣乘法,以及為什麼可分離(separable)的二維轉換很省計算。
  • DCT、沃爾什–哈達瑪、斜轉換(slant)與哈爾轉換的基底影像長什麼樣子,以及 DCT 為何是 JPEG 的主力。
  • 如何用時頻平面(time–frequency plane)理解傅立葉分析與小波分析的差別。
  • 多解析度分析、尺度函數與小波函數如何導出快速小波轉換,也就是一個由兩個濾波器組成的濾波器組。
  • 如何用 PyWavelets 做二維離散小波轉換,並應用在去雜訊、邊緣偵測與壓縮。

先看全貌

第 4 章把影像拆成許多正弦波的總和。這只是其中一種積木,並不是唯一的選擇。本章退一步問一個更一般的問題:給定任意一組積木,我要怎麼算出影像裡每塊積木各有多少?又該挑哪一組積木最適合手上的工作?

「挑積木」的嚴謹說法是線性代數:內積、基底與基底轉換。這套工具一旦建立,本章的每一種轉換都只是一個不同的矩陣。小波會占最多篇幅,因為它是第一種在空間與頻率兩邊同時具有局部性的基底。這個特性讓小波成為 JPEG 2000 標準的核心 [15],也是現代稀疏表示,以及近年部分神經網路層背後的重要想法。

預備知識

白話版:一個有 NN 個取樣的訊號,就是 NN 維空間中的一支箭頭(向量)。轉換就是量這支箭頭在 NN 支參考箭頭上各投下多長的影子。如果參考箭頭兩兩垂直且長度為 1,影子就包含全部資訊,把影子加回去就能還原原本的箭頭。

內積空間

內積(inner product)⟨f,g⟩\langle \mathbf{f}, \mathbf{g} \rangle 接收兩個向量、回傳一個數,用來衡量兩者指向同一方向的程度。對長度為 NN 的實數或複數向量:

⟨f,g⟩=∑x=0N−1f(x) g∗(x),∥f∥=⟨f,f⟩.\langle \mathbf{f}, \mathbf{g} \rangle = \sum_{x=0}^{N-1} f(x)\, g^*(x), \qquad \|\mathbf{f}\| = \sqrt{\langle \mathbf{f}, \mathbf{f} \rangle}.

其中 f(x)f(x)、g(x)g(x) 是 f\mathbf{f}、g\mathbf{g} 的第 xx 個元素,g∗g^* 是複數共軛(對實數資料沒有作用),∥f∥\|\mathbf{f}\| 是 f\mathbf{f} 的長度,也就是範數(norm)。對連續函數,總和改成積分:⟨f,g⟩=∫f(x)g∗(x) dx\langle f, g \rangle = \int f(x) g^*(x)\,dx。當兩個向量的內積為零時稱為正交(orthogonal):它們彼此毫無共同成分。

正交歸一基底

若一組向量 {s0,…,sN−1}\{\mathbf{s}_0, \dots, \mathbf{s}_{N-1}\} 每一個長度都是 1,而且兩兩正交,即 ⟨su,sv⟩=δ(u−v)\langle \mathbf{s}_u, \mathbf{s}_v \rangle = \delta(u - v)(u=vu = v 時 δ\delta 為 1,否則為 0),就稱為正交歸一基底(orthonormal basis)。此時任何向量都能乾淨地拆開:

t(u)=⟨f,su⟩,f=∑u=0N−1t(u) su.t(u) = \langle \mathbf{f}, \mathbf{s}_u \rangle, \qquad \mathbf{f} = \sum_{u=0}^{N-1} t(u)\, \mathbf{s}_u .

t(u)t(u) 稱為轉換係數。第一式是分析(正向)步驟,第二式是合成(反向)步驟。正交歸一基底還會保持能量:∑x∣f(x)∣2=∑u∣t(u)∣2\sum_x |f(x)|^2 = \sum_u |t(u)|^2(帕塞瓦定理,Parseval’s theorem),所以我們才能談「每個係數裝了多少能量」。

雙正交基底

有時最好用的積木並不互相垂直。這時需要兩組向量:一組 {s~u}\{\tilde{\mathbf{s}}_u\} 負責量測(分析),一組 {su}\{\mathbf{s}_u\} 負責重建(合成)。若滿足

⟨su,s~v⟩=δ(u−v),t(u)=⟨f,s~u⟩,f=∑ut(u) su,\langle \mathbf{s}_u, \tilde{\mathbf{s}}_v \rangle = \delta(u - v), \qquad t(u) = \langle \mathbf{f}, \tilde{\mathbf{s}}_u \rangle, \qquad \mathbf{f} = \sum_u t(u)\, \mathbf{s}_u ,

就稱為雙正交(biorthogonal)基底對,量測用的 s~u\tilde{\mathbf{s}}_u 稱為對偶(dual)基底。雙正交換來的是設計自由度:例如 JPEG 2000 採用雙正交小波濾波器,可以同時做到長度短且左右對稱,這個組合用正交小波很難達成 [15]。

框架

框架(frame)是一組向量數量多於維度的生成集合,所以帶有冗餘。它的定義性質是:量到的能量被固定的上下界夾住,

A ∥f∥2  ≤  ∑k∣⟨f,sk⟩∣2  ≤  B ∥f∥2,0<A≤B<∞.A\,\|\mathbf{f}\|^2 \;\le\; \sum_k |\langle \mathbf{f}, \mathbf{s}_k \rangle|^2 \;\le\; B\,\|\mathbf{f}\|^2, \qquad 0 < A \le B < \infty .

AA、BB 稱為框架界(frame bounds)。若 A=BA = B,則稱為緊框架(tight frame),重建公式就是 f=1A∑k⟨f,sk⟩sk\mathbf{f} = \frac{1}{A}\sum_k \langle \mathbf{f}, \mathbf{s}_k \rangle \mathbf{s}_k。一個經典的例子:平面上三支夾角 120° 的單位箭頭。它們能張成整個平面,任兩支都不正交,卻構成 A=B=3/2A = B = 3/2 的緊框架。當你需要穩健性或平移不變性時,多花一些係數換冗餘是值得的,例如不抽樣小波轉換(undecimated wavelet transform),以及「現代觀點」會提到的散射轉換。

import numpy as np

v = np.array([3.0, 1.0, 2.0])
e = np.eye(3)                         # standard orthonormal basis
c = e @ v                             # coefficients = inner products <v, e_k>
print(c, np.allclose(c @ e, v))       # synthesis: sum_k c_k e_k

# A biorthogonal pair in 2-D: analysis rows A, synthesis columns B = A^{-1}
A = np.array([[1.0, 0.0], [1.0, 1.0]])
B = np.linalg.inv(A)
x = np.array([2.0, 5.0])
print(A @ x, B @ (A @ x))             # [2. 7.] -> back to [2. 5.]

以矩陣表示的轉換

白話版:把量測用的向量疊成矩陣的每一列。訊號乘上這個矩陣,就一次量完所有係數;再乘上反矩陣,就把訊號拼回來。

對長度為 NN 的一維訊號 f\mathbf{f},把基底向量 su∗T\mathbf{s}_u^{*T} 依序放進 N×NN \times N 轉換矩陣 A\mathbf{A} 的各列,則

t=A f,f=A−1t.\mathbf{t} = \mathbf{A}\,\mathbf{f}, \qquad \mathbf{f} = \mathbf{A}^{-1}\mathbf{t} .

若基底是正交歸一的,A\mathbf{A} 就是么正矩陣(unitary,A−1=A∗T\mathbf{A}^{-1} = \mathbf{A}^{*T};實數情形稱為正交矩陣,A−1=AT\mathbf{A}^{-1} = \mathbf{A}^{T}),反轉換幾乎不用成本。

對 M×NM \times N 的影像 F\mathbf{F},一般的二維線性轉換需要 MN×MNMN \times MN 的巨大矩陣。實務上幾乎所有影像轉換都是可分離的:先對每一行做一維轉換,再對每一列做一維轉換。寫成矩陣:

T=A F BT,F=A−1 T (BT)−1,\mathbf{T} = \mathbf{A}\,\mathbf{F}\,\mathbf{B}^{T}, \qquad \mathbf{F} = \mathbf{A}^{-1}\,\mathbf{T}\,(\mathbf{B}^{T})^{-1} ,

其中 A\mathbf{A}(M×MM \times M)作用在行上,B\mathbf{B}(N×NN \times N)作用在列上;方形影像且兩方向使用同一種轉換時 B=A\mathbf{B} = \mathbf{A}。可分離性把直接計算 N×NN \times N 轉換的成本從 O(N4)O(N^4) 降到 O(N3)O(N^3),而快速演算法(FFT、快速 DCT、快速小波轉換)還能再進一步降低。

import numpy as np
from scipy.fft import dct, dctn
from skimage import data, img_as_float

N = 8
A = dct(np.eye(N), norm="ortho", axis=0)  # row u = u-th cosine basis vector

print(np.allclose(A @ A.T, np.eye(N)))    # orthonormal -> True

F = img_as_float(data.camera())[200:208, 200:208]
T = A @ F @ A.T                           # forward 2-D transform
print(np.allclose(T, dctn(F, norm="ortho")))   # same as scipy -> True
print(np.allclose(A.T @ T @ A, F))             # inverse -> True

相關運算

白話版:每個係數都在回答一個問題:「我的影像有多像這個圖樣?」這個問題本身就是相關(correlation)運算。

把分析公式和第 3 章的互相關比一比。係數

t(u)=⟨f,su⟩=∑xf(x) su∗(x)t(u) = \langle \mathbf{f}, \mathbf{s}_u \rangle = \sum_{x} f(x)\, s_u^*(x)

正是 ff 與基底函數 sus_u 在位移為零時的相關值。∣t(u)∣|t(u)| 大代表 ff 和 sus_u 很像;t(u)=0t(u) = 0 代表 ff 完全不含 sus_u 的成分。由此得到兩個結論:

  • 轉換就是樣板比對。 每個基底函數都是一個樣板,轉換一次替所有樣板打分數。傅立葉轉換替正弦波打分數,哈爾轉換替階梯狀圖樣打分數。
  • 小波是在許多位移與尺度上做相關。 小波基底包含同一個圖樣的各種平移與伸縮版本。和某一尺度的所有平移版本做相關,就是一次濾波,所以本章後面的小波轉換會變成一組濾波器。

時頻平面上的基底函數

白話版:有些積木能告訴你事情何時發生,有些能告訴你它的音高,但沒有一塊能把兩者都說得完美。時頻平面就是描繪這種取捨的地圖。

對一個能量為 1 的基底函數 s(x)s(x),用 ∣s(x)∣2|s(x)|^2 與 ∣S(ω)∣2|S(\omega)|^2 的標準差定義它在時間(空間)與頻率上的展開程度:

σx2=∫(x−xˉ)2∣s(x)∣2 dx,σω2=12π∫(ω−ωˉ)2∣S(ω)∣2 dω,\sigma_x^2 = \int (x - \bar{x})^2 |s(x)|^2\,dx, \qquad \sigma_\omega^2 = \frac{1}{2\pi}\int (\omega - \bar{\omega})^2 |S(\omega)|^2\,d\omega ,

其中 S(ω)S(\omega) 是 ss 的傅立葉轉換,xˉ\bar{x}、ωˉ\bar{\omega} 分別是兩個能量分布的質心。訊號版的海森堡測不準原理(Heisenberg uncertainty principle)指出

σx σω  ≥  12,\sigma_x\, \sigma_\omega \;\ge\; \tfrac{1}{2},

只有高斯形狀的函數能取到等號。把每個基底函數畫成一個以 (xˉ,ωˉ)(\bar{x}, \bar{\omega}) 為中心、寬 σx\sigma_x、高 σω\sigma_\omega 的矩形,稱為海森堡盒(Heisenberg box)或時頻磚。磚的面積不能小於某個下限,但形狀可以自由選擇。

四種時頻鋪磚方式:取樣是直條、傅立葉是橫條、短時傅立葉是均勻方格、小波是高頻處窄而高、低頻處寬而矮的二進位鋪法。
圖 7.1 — 四種基底如何鋪滿時頻平面,每塊磚的面積都相同。取樣在時間上最精確;傅立葉/DCT 基底在頻率上最精確;短時傅立葉轉換對所有頻率採用同一種折衷;小波則依頻率調整磚的形狀。

圖 7.1 呈現四種標準情形:

  • 取樣(單位基底): 位置完全精確,沒有任何頻率資訊。
  • 傅立葉或 DCT: 頻率完全精確,沒有位置資訊。一條銳利的邊緣會擴散到所有係數上。
  • 短時傅立葉轉換(short-time Fourier transform, STFT):把訊號切成一個個窗,各自做傅立葉轉換。每塊磚形狀相同,一次決定就套用到所有頻率 [2]。
  • 小波: 高頻的磚在時間上窄、頻率上高;低頻的磚則寬而矮。這很符合自然影像:邊緣是短暫的高頻事件,平緩的明暗變化則是緩慢的低頻成分。

基底影像

白話版:二維轉換的每個係數都對應一張小圖。影像就是這些小圖的加權總和,係數就是權重。

對矩陣 A\mathbf{A}(第 uu 列為 auT\mathbf{a}_u^T)定義的可分離轉換,可以把反轉換改寫成外積的總和:

F=∑u=0N−1∑v=0N−1T(u,v) Su,v,Su,v=au∗av∗T.\mathbf{F} = \sum_{u=0}^{N-1} \sum_{v=0}^{N-1} T(u, v)\, \mathbf{S}_{u,v}, \qquad \mathbf{S}_{u,v} = \mathbf{a}_u^{*} \mathbf{a}_v^{*T} .

每個 N×NN \times N 矩陣 Su,v\mathbf{S}_{u,v} 就是一張基底影像(basis image)。T(u,v)T(u,v) 告訴你要加入多少第 (u,v)(u, v) 張基底影像;uu 對應垂直方向的圖樣,vv 對應水平方向。把 N2N^2 張基底影像排在一起,是理解一種轉換「看到什麼」最快的方法。

DCT、沃爾什–哈達瑪、斜轉換與哈爾轉換各自的 8×8 組 8×8 基底影像。DCT 是平滑的餘弦;沃爾什–哈達瑪是黑白棋盤;斜轉換含有線性斜坡;哈爾則是灰色零背景上的局部方塊。
圖 7.2 — 四種轉換的 8×8 基底影像,向下代表垂直頻率(序率)增加,向右代表水平頻率增加;每格縮放到 [−1, 1]。注意哈爾基底影像和其他三種不同,大部分區域為零(灰色):它們是局部的。

四種轉換左上角的基底影像都是常數,它量的是平均亮度,稱為 DC 係數。越往右、往下,圖樣振盪越快。DCT 平滑地振盪,沃爾什–哈達瑪在 ±1\pm 1 之間突然切換,斜轉換多了線性斜坡,哈爾則把每個圖樣限制在格子中的一小塊區域。

與傅立葉相關的轉換

白話版:傅立葉家族用波當積木。其中餘弦轉換在區塊邊界的表現最好,所以你手機裡的 JPEG 檔用的就是它。

以矩陣看離散傅立葉轉換

第 4 章的一維 DFT 是一個矩陣轉換:

Au,x=1N e−j2πux/N,u,x=0,…,N−1.A_{u,x} = \frac{1}{\sqrt{N}}\, e^{-j 2\pi u x / N}, \qquad u, x = 0, \dots, N-1 .

Au,xA_{u,x} 是第 uu 列第 xx 行的元素,j=−1j = \sqrt{-1}。乘上 1/N1/\sqrt{N} 後矩陣是么正的。對壓縮來說它有兩個缺點:係數是複數,而且 DFT 默默把訊號當成週期性的。若區塊左右邊界的值不同,週期延伸就會出現跳躍,而跳躍需要大量高頻係數來表示。

離散餘弦轉換

(第二型)離散餘弦轉換(discrete cosine transform, DCT)[3] 使用實數餘弦:

t(u)=α(u)∑x=0N−1f(x)cos⁡ ⁣[(2x+1) u π2N],α(0)=1N,    α(u)=2N  (u≥1).t(u) = \alpha(u) \sum_{x=0}^{N-1} f(x) \cos\!\left[\frac{(2x+1)\,u\,\pi}{2N}\right], \qquad \alpha(0) = \sqrt{\tfrac{1}{N}},\;\; \alpha(u) = \sqrt{\tfrac{2}{N}} \;(u \ge 1) .

uu 是頻率索引,α(u)\alpha(u) 讓基底正交歸一。DCT 等於(差一個縮放)把訊號鏡射成長度 2N2N 的偶對稱延伸後再做 DFT。鏡射消除了邊界跳躍,能量因此集中在少數低頻係數。這種能量集中(energy compaction)特性,正是 JPEG 標準對 8×88 \times 8 影像區塊做 DCT 的原因 [4]。Ahmed、Natarajan 與 Rao 在 1974 年提出 DCT,並指出它在維納濾波與率失真(rate–distortion)準則下的表現,與統計上最佳的 Karhunen–Loève 轉換非常接近 [3]。

離散正弦轉換

離散正弦轉換(discrete sine transform, DST)改用正弦。第一型 DST 的矩陣為

Au,x=2N+1 sin⁡ ⁣[(x+1)(u+1)πN+1],A_{u,x} = \sqrt{\frac{2}{N+1}}\, \sin\!\left[\frac{(x+1)(u+1)\pi}{N+1}\right] ,

對應的是訊號的奇對稱延伸。它是實數、正交歸一且對稱的矩陣(A=AT=A−1\mathbf{A} = \mathbf{A}^T = \mathbf{A}^{-1}),適合兩端接近零的訊號,例如預測殘差。在 SciPy 中可寫成 scipy.fft.dst(x, type=1, norm="ortho")。

沃爾什–哈達瑪轉換

白話版:每個基底向量都只由 +1+1 和 −1-1 組成。這樣轉換完全不需要乘法,只要加減。

N=2nN = 2^n 階的哈達瑪矩陣以遞迴方式建構:

H1=[ 1 ],H2N=[HNHNHN−HN],AWHT=1N HN.\mathbf{H}_1 = [\,1\,], \qquad \mathbf{H}_{2N} = \begin{bmatrix} \mathbf{H}_N & \mathbf{H}_N \\ \mathbf{H}_N & -\mathbf{H}_N \end{bmatrix}, \qquad \mathbf{A}_{\text{WHT}} = \frac{1}{\sqrt{N}}\,\mathbf{H}_N .

HN\mathbf{H}_N 的每個元素都是 ±1\pm 1,各列互相正交,因此 AWHT\mathbf{A}_{\text{WHT}} 是正交歸一(而且對稱)的。HN\mathbf{H}_N 的列天生是自然順序(哈達瑪順序)。做分析時,改依序率(sequency,一列中正負號改變的次數,扮演頻率的角色)排序會更好讀。依序率排列的列稱為沃爾什函數,圖 7.2 就是它們的二維基底影像。

沃爾什–哈達瑪轉換(Walsh–Hadamard transform, WHT)因為計算便宜,很早就被用於影像編碼 [5]。對自然影像而言,它的能量集中能力不如 DCT,因為方波基底不擅長逼近平滑的明暗變化(見後面的圖 7.7)。在講求速度或硬體簡單的場合,以及作為快速隨機化演算法的元件時,它仍然很有用。

import numpy as np
from scipy.linalg import hadamard

H = hadamard(8)                                    # natural (Hadamard) order
changes = (np.diff(np.sign(H), axis=1) != 0).sum(1)
W = H[np.argsort(changes)] / np.sqrt(8)            # sequency order
print((np.diff(np.sign(W), axis=1) != 0).sum(1))   # [0 1 2 3 4 5 6 7]

x = np.array([10, 12, 11, 13, 40, 42, 41, 43], float)
print(np.round(W @ x, 2))

這個測試訊號是一個帶有小漣漪的階梯。它的 WHT 幾乎全落在係數 0(平均值)與係數 1(恰好一次正負號變化,對應那個階梯);其餘只剩幾個小數值。

斜轉換

白話版:影像中很多區域不是平的,而是逐漸變亮或變暗。斜轉換準備了一塊「斜坡」積木,讓這種區域只需要一兩個係數。

Pratt、Chen 與 Welch 為影像編碼設計了斜轉換(slant transform)[6]。它和 WHT 一樣是正交歸一、遞迴建構的,從 S2=12[111−1]\mathbf{S}_2 = \frac{1}{\sqrt{2}}\begin{bmatrix}1 & 1 \\ 1 & -1\end{bmatrix} 開始:

SN=12 QN[SN/200SN/2],\mathbf{S}_{N} = \frac{1}{\sqrt{2}}\, \mathbf{Q}_N \begin{bmatrix} \mathbf{S}_{N/2} & \mathbf{0} \\ \mathbf{0} & \mathbf{S}_{N/2} \end{bmatrix} ,

其中 QN\mathbf{Q}_N 是一個稀疏的 N×NN \times N 混合矩陣,內含兩個常數

aN=3N24(N2−1),bN=N2−44(N2−1),a_N = \sqrt{\frac{3N^2}{4(N^2-1)}}, \qquad b_N = \sqrt{\frac{N^2-4}{4(N^2-1)}} ,

它們的取值讓 SN\mathbf{S}_N 的第二列成為等差遞減的階梯(「斜」向量),同時保持矩陣正交歸一。本章的繪圖腳本(scripts/figures/dip_ch07.py)就是這樣建出 S8\mathbf{S}_8,並以數值檢查這兩個性質。圖 7.2 中,斜轉換第一列與第一行的基底影像都是看得見的斜坡。斜轉換有 O(Nlog⁡2N)O(N \log_2 N) 的快速演算法,但實務上同樣能處理斜坡的 DCT 取代了它。

哈爾轉換

白話版:把相鄰的兩個值拿來,記下它們的平均和差。再對平均值重複同樣的步驟。這就是哈爾轉換,也是最簡單的小波。

哈爾轉換(Haar transform)[1] 的基底函數要不是常數,就是一個「先上後下」的階梯,具有不同的寬度與位置。對 N=2nN = 2^n 個取樣,用尺度 pp 與位置 qq 來編號,k≥1k \ge 1 時 k=2p+qk = 2^p + q。不計歸一化,第 kk 個基底向量在區間 [q2p,q+12p)\left[\frac{q}{2^p}, \frac{q+1}{2^p}\right) 的前半段為 +1+1、後半段為 −1-1,其餘為 0(位置以訊號長度的比例表示);k=0k = 0 的向量是常數。N=4N = 4 時:

AHaar=[121212121212−12−1212−12000012−12].\mathbf{A}_{\text{Haar}} = \begin{bmatrix} \tfrac{1}{2} & \tfrac{1}{2} & \tfrac{1}{2} & \tfrac{1}{2} \\[2pt] \tfrac{1}{2} & \tfrac{1}{2} & -\tfrac{1}{2} & -\tfrac{1}{2} \\[2pt] \tfrac{1}{\sqrt{2}} & -\tfrac{1}{\sqrt{2}} & 0 & 0 \\[2pt] 0 & 0 & \tfrac{1}{\sqrt{2}} & -\tfrac{1}{\sqrt{2}} \end{bmatrix} .

第 0、1 列看整段訊號;第 2、3 列各自只看一半。這是它和前面所有轉換最關鍵的不同:哈爾基底函數是局部的。邊緣這類局部事件只會改變覆蓋到它的少數幾個哈爾係數。在時頻圖上,哈爾正好對應圖 7.1 的二進位(dyadic)鋪法。它是通往小波轉換的橋梁:以哈爾小波做的離散小波轉換,就是哈爾轉換。

小波轉換

白話版:用許多種縮放倍率看同一張影像。每往下一層,保留一張更模糊的版本,並只記下模糊時被抹掉的細節。模糊版本疊成一座金字塔,而那些細節就是小波係數。

多解析度分析

多解析度分析(multiresolution analysis, MRA)由 Mallat 形式化 [7],把上面的想法變成一組基底。先取一個尺度函數(scaling function)φ(x)\varphi(x),做出平移與伸縮的版本:

φj,k(x)=2j/2 φ(2jx−k),\varphi_{j,k}(x) = 2^{j/2}\, \varphi(2^j x - k) ,

其中 jj 是尺度(jj 越大,函數越窄、細節越精細),kk 是整數平移量,2j/22^{j/2} 讓能量維持為 1。令 VjV_j 為 {φj,k}k\{\varphi_{j,k}\}_k 所張成的空間,也就是在解析度 jj 下能表示的所有訊號。MRA 要求四件事:

  1. φ0,k\varphi_{0,k} 彼此正交歸一(或至少是 V0V_0 的穩定基底)。
  2. 空間是巢狀的,⋯⊂V−1⊂V0⊂V1⊂⋯\cdots \subset V_{-1} \subset V_0 \subset V_1 \subset \cdots:粗尺度看得到的東西,細尺度也看得到。
  3. 所有 VjV_j 唯一的共同元素是 f(x)=0f(x) = 0。
  4. 當 j→∞j \to \infty,任何平方可積函數都能被逼近到任意精確。

因為 V0⊂V1V_0 \subset V_1,尺度函數本身必須能由更細的版本組合出來,於是得到精細化方程式(refinement equation,亦稱 dilation equation):

φ(x)=∑nhφ(n) 2 φ(2x−n),\varphi(x) = \sum_n h_\varphi(n)\, \sqrt{2}\, \varphi(2x - n) ,

係數 hφ(n)h_\varphi(n) 稱為尺度函數係數,稍後它會變成低通濾波器。

小波函數

從 Vj+1V_{j+1} 降到 VjV_j 時遺失的細節,落在互補空間 WjW_j 中,使得 Vj+1=Vj⊕WjV_{j+1} = V_j \oplus W_j(⊕\oplus 表示 Vj+1V_{j+1} 中每個元素都能唯一拆成 VjV_j 的一部分加上 WjW_j 的一部分;對正交小波,兩部分互相正交)。WjW_j 由小波函數(wavelet function)張成:

ψj,k(x)=2j/2 ψ(2jx−k),ψ(x)=∑nhψ(n) 2 φ(2x−n).\psi_{j,k}(x) = 2^{j/2}\, \psi(2^j x - k), \qquad \psi(x) = \sum_n h_\psi(n)\, \sqrt{2}\, \varphi(2x - n) .

對正交小波,小波函數係數可由尺度係數經調變與時間反轉得到:hψ(n)=(−1)n hφ(1−n)h_\psi(n) = (-1)^n\, h_\varphi(1 - n),它們構成高通濾波器。一再重複這個拆分,得到

L2(R)=Vj0⊕Wj0⊕Wj0+1⊕⋯ ,L^2(\mathbb{R}) = V_{j_0} \oplus W_{j_0} \oplus W_{j_0+1} \oplus \cdots ,

其中 L2(R)L^2(\mathbb{R}) 是所有有限能量函數構成的空間:一個粗略的近似,加上每個更細尺度的細節。

哈爾小波的 hφ=[12,12]h_\varphi = [\tfrac{1}{\sqrt{2}}, \tfrac{1}{\sqrt{2}}]、hψ=[12,−12]h_\psi = [\tfrac{1}{\sqrt{2}}, -\tfrac{1}{\sqrt{2}}]:就是平均與差。Daubechies 展示了如何建構具有緊支撐(compact support,即有限長度濾波器)且平滑度可調的正交小波,其正則性隨濾波器長度線性增加 [8]。PyWavelets [10] 中的 dbN 家族就是這些小波,N 是消失矩(vanishing moments)的個數:對 m=0,…,N−1m = 0, \dots, N-1 都有 ∫xmψ(x) dx=0\int x^m \psi(x)\,dx = 0。具有 NN 個消失矩的小波對次數低於 NN 的多項式完全「視而不見」,因此平滑區域的細節係數幾乎為零。

haar、db2、db4 的尺度函數與小波圖形。haar 是方塊與階梯;db2 鋸齒狀;db4 較平滑也較長。
圖 7.3 — 以 PyWavelets 的級聯演算法(cascade algorithm)算出的尺度函數(上)與小波(下):Haar、db2、db4。濾波器越長,函數越平滑、也越長。

小波級數展開

有了這兩族函數,連續訊號 f(x)f(x) 可以展開成

f(x)=∑kcj0(k) φj0,k(x)+∑j=j0∞∑kdj(k) ψj,k(x),f(x) = \sum_k c_{j_0}(k)\, \varphi_{j_0,k}(x) + \sum_{j=j_0}^{\infty} \sum_k d_j(k)\, \psi_{j,k}(x) , cj0(k)=⟨f,φj0,k⟩,dj(k)=⟨f,ψj,k⟩.c_{j_0}(k) = \langle f, \varphi_{j_0,k} \rangle, \qquad d_j(k) = \langle f, \psi_{j,k} \rangle .

j0j_0 是任選的起始(最粗)尺度。cj0(k)c_{j_0}(k) 稱為近似(或尺度)係數,dj(k)d_j(k) 稱為細節(或小波)係數。注意每個係數又是一個內積,與預備知識中完全相同。Unser 與 Blu 說明了消失矩、逼近階數與平滑度這些性質在此展開式中從何而來 [9]。

一維離散小波轉換

對取樣訊號 f(x)f(x),x=0,…,M−1x = 0, \dots, M - 1,M=2JM = 2^J,總和變成有限項:

Wφ(j0,k)=1M∑xf(x) φj0,k(x),Wψ(j,k)=1M∑xf(x) ψj,k(x),j≥j0,W_\varphi(j_0, k) = \frac{1}{\sqrt{M}} \sum_{x} f(x)\, \varphi_{j_0,k}(x), \qquad W_\psi(j, k) = \frac{1}{\sqrt{M}} \sum_{x} f(x)\, \psi_{j,k}(x), \quad j \ge j_0 , f(x)=1M[∑kWφ(j0,k) φj0,k(x)+∑j=j0J−1∑kWψ(j,k) ψj,k(x)].f(x) = \frac{1}{\sqrt{M}} \left[ \sum_k W_\varphi(j_0, k)\, \varphi_{j_0,k}(x) + \sum_{j=j_0}^{J-1} \sum_k W_\psi(j, k)\, \psi_{j,k}(x) \right] .

WφW_\varphi、WψW_\psi 是離散小波轉換(discrete wavelet transform, DWT)的近似係數與細節係數,此處的 φj,k(x)\varphi_{j,k}(x)、ψj,k(x)\psi_{j,k}(x) 是在 xx 取樣的連續函數。對正交小波,係數總數恰好等於 MM:DWT 是一次基底轉換,而不是擴張。

快速小波轉換

白話版:你完全不需要去算 φ\varphi 或 ψ\psi 的值。兩個短濾波器加上「每隔一個取一個」,一層一層做下去就夠了。

把精細化方程式代入 DWT,就得到 Mallat 的金字塔演算法,今天稱為快速小波轉換(fast wavelet transform, FWT)[7]:

Wψ(j,k)=∑nhψ(n−2k) Wφ(j+1,n),Wφ(j,k)=∑nhφ(n−2k) Wφ(j+1,n).W_\psi(j, k) = \sum_n h_\psi(n - 2k)\, W_\varphi(j+1, n), \qquad W_\varphi(j, k) = \sum_n h_\varphi(n - 2k)\, W_\varphi(j+1, n) .

讀法:把較細一層的近似 Wφ(j+1,⋅)W_\varphi(j+1, \cdot) 與 hψh_\psi(或 hφh_\varphi)做相關,再每隔一個取一個輸出(2 倍降取樣)。用訊號處理的語言,這就是雙通道分析濾波器組(analysis filter bank):低通分支產生下一層近似,高通分支產生細節。把低通輸出再送進同一組濾波器,如此反覆。每一層資料量減半,所以總成本是 O(M)O(M),比 FFT 的 O(Mlog⁡M)O(M \log M) 還快。

反向 FWT 是合成濾波器組:每個分支先 2 倍升取樣(插零)、經合成濾波器濾波,再相加。對正交小波,合成濾波器就是分析濾波器的時間反轉,整組濾波器能達成完美重建(perfect reconstruction):即使每個分支的降取樣都造成混疊(aliasing),兩個分支的混疊項會互相抵消,輸出與輸入完全相同。Vetterli 與 Kovačević 的教科書從子頻帶編碼(subband coding)的角度深入發展了這個觀點 [2]。

import numpy as np
import pywt

x = np.array([4.0, 6.0, 10.0, 12.0, 8.0, 6.0, 5.0, 5.0])
w = pywt.Wavelet("haar")
h0, h1 = np.array(w.dec_lo), np.array(w.dec_hi)   # [.707 .707], [-.707 .707]

def analysis(x):
    lo = np.convolve(x, h0)[1::2]   # filter, then keep every 2nd sample
    hi = np.convolve(x, h1)[1::2]
    return lo, hi

a1, d1 = analysis(x)
a2, d2 = analysis(a1)
print("a2", a2, "d2", d2, "d1", d1.round(3))

ref = pywt.wavedec(x, "haar", level=2)            # [a2, d2, d1]
print(all(np.allclose(p, q) for p, q in zip(ref, [a2, d2, d1])))
print(np.allclose(pywt.waverec(ref, "haar"), x))

這兩行濾波器組與 pywt.wavedec 的結果完全一致。注意 (5, 5) 這一對的細節係數為 0:沒有變化,就沒有細節。

二維小波轉換

白話版:先沿著列跑一維濾波器組,再沿著行跑一次。你會得到一張小而模糊的影像,以及三張分別對應不同方向的「邊緣圖」。

二維 DWT 是可分離的。一層分解使用一個二維尺度函數與三個二維小波,各由一維函數相乘而成:

φ(x,y)=φ(x)φ(y),ψH(x,y)=ψ(x)φ(y),ψV(x,y)=φ(x)ψ(y),ψD(x,y)=ψ(x)ψ(y).\varphi(x, y) = \varphi(x)\varphi(y), \quad \psi^H(x, y) = \psi(x)\varphi(y), \quad \psi^V(x, y) = \varphi(x)\psi(y), \quad \psi^D(x, y) = \psi(x)\psi(y) .

依照課本慣例,xx 是列索引(垂直位置),yy 是行索引。ψH\psi^H 沿垂直方向變化、水平方向平滑,所以對水平邊緣有反應;ψV\psi^V 對垂直邊緣有反應;ψD\psi^D 則對對角細節與角點有反應。一層分解把 M×NM \times N 影像變成四個 M2×N2\frac{M}{2} \times \frac{N}{2} 的子頻帶(subband):

子頻帶垂直方向濾波水平方向濾波常見稱呼PyWavelets 名稱
近似低通低通LLcA
水平細節高通低通LH 或 HL(慣例不一)cH
垂直細節低通高通HL 或 LHcV
對角細節高通高通HHcD

由於不同教科書與函式庫對 LH/HL 的標記剛好相反,本章一律使用不會混淆的 H、V、D。對近似子頻帶重複分解,就得到圖 7.4 那種熟悉的巢狀排列。

import numpy as np
import pywt
from skimage import data, img_as_float

img = img_as_float(data.camera())
LL1, (H1, V1, D1) = pywt.dwt2(img, "haar")   # one level
LL2, (H2, V2, D2) = pywt.dwt2(LL1, "haar")   # second level on LL1
print(img.shape, LL1.shape, LL2.shape)       # (512, 512) (256, 256) (128, 128)

e = lambda a: float((a ** 2).sum())
total = e(img)
print(f"energy in LL2: {e(LL2) / total:.4f}")

coeffs = pywt.wavedec2(img, "haar", level=2)  # same thing in one call
print(np.allclose(pywt.waverec2(coeffs, "haar"), img))

兩層之後,128×128128 \times 128 的近似子頻帶只占係數的十六分之一,卻裝了 camera 影像約 99% 的能量。其餘 15/16 的係數大多接近零,這正是壓縮與去雜訊所利用的稀疏性。

左:camera 測試影像。右:兩層哈爾分解,左上角是小張的近似影像,其餘細節子頻帶顯示水平、垂直與對角方向的邊緣。
圖 7.4 — camera 影像的兩層哈爾 DWT。細節子頻帶顯示的是絕對值,且各自以其 99.5 百分位數重新縮放,好讓小係數看得見;若用同一個比例尺,它們幾乎全黑。三腳架垂直的腳在 V1 中特別亮,建築物的水平輪廓則出現在 H1。

小波封包

標準 DWT 只會再拆分低通子頻帶。小波封包(wavelet packet)分解則連細節子頻帶也全部再拆,形成一棵完整的二元樹(二維時是四元樹)。LL 層之後,二維封包樹有 4L4^L 片葉子。任何一組能不重疊地覆蓋頻率軸的節點都是合法的正交歸一基底,所以這棵樹是一座基底圖書館。Coifman 與 Wickerhauser 的最佳基底(best basis)演算法能有效率地搜尋這座圖書館,挑出使某個可加成本(例如歸一化係數的熵)最小的基底 [11]。封包特別適合紋理,因為紋理在中高頻有大量能量,而一般 DWT 只把它們放在一個寬寬的子頻帶裡。

import pywt
from skimage import data, img_as_float

img = img_as_float(data.camera())
wp = pywt.WaveletPacket2D(img, "haar", maxlevel=2)
print([n.path for n in wp.get_level(2)][:6], len(wp.get_level(2)))  # 16 subbands

節點路徑記錄每一層走的分支(a 為近似,h、v、d 為細節),所以 "hd" 是水平細節子頻帶裡的對角細節。

應用

以門檻值去雜訊

白話版:在小波域裡,乾淨的影像只是少數幾個大數字,白雜訊則是散布各處的許多小數字。把小的歸零,再轉回來就好。

正交轉換會把標準差為 σ\sigma 的白高斯雜訊,變成每個子頻帶中標準差同樣為 σ\sigma 的白高斯雜訊,而影像的能量則集中在少數係數上。對細節係數 ww 與門檻值 TT,兩種門檻規則為:

ηhard(w)={w,∣w∣>T0,∣w∣≤Tηsoft(w)=sign⁡(w) max⁡(∣w∣−T, 0).\eta_{\text{hard}}(w) = \begin{cases} w, & |w| > T \\ 0, & |w| \le T \end{cases} \qquad \eta_{\text{soft}}(w) = \operatorname{sign}(w)\, \max(|w| - T,\, 0) .

硬門檻(hard thresholding)只決定留或砍;軟門檻(soft thresholding)還會把留下的係數往零縮 TT,得到連續的規則,假影較少,但對比會稍微損失。兩種經典的 TT:

  • 通用門檻值(universal threshold,VisuShrink):nn 個取樣時 T=σ2ln⁡nT = \sigma\sqrt{2 \ln n},出自 Donoho 與 Johnstone [12]。Donoho 證明用此 TT 做軟門檻,所得估計有很高機率至少與真實訊號一樣平滑 [13]。雜訊強度通常由最細的對角子頻帶穩健地估計:σ^=median⁡(∣w∣)/0.6745\hat{\sigma} = \operatorname{median}(|w|)/0.6745 [12]。
  • BayesShrink:Chang、Yu 與 Vetterli 為每個子頻帶設定各自的門檻 T=σ^2/σ^XT = \hat{\sigma}^2 / \hat{\sigma}_X,其中 σ^X\hat{\sigma}_X 是該子頻帶中乾淨係數的估計標準差 [14]。它會依各子頻帶含有多少訊號自動調整,通常優於通用門檻值。
import numpy as np
import pywt
from skimage import data, img_as_float
from skimage.metrics import peak_signal_noise_ratio as psnr

rng = np.random.default_rng(7)
clean = img_as_float(data.camera())
noisy = clean + rng.normal(0, 0.1, clean.shape)

def wavelet_denoise(y, wavelet="db4", level=3, mode="soft"):
    c = pywt.wavedec2(y, wavelet, mode="periodization", level=level)
    sigma = np.median(np.abs(c[-1][2])) / 0.6745       # noise level from finest HH
    T = sigma * np.sqrt(2 * np.log(y.size))            # universal threshold
    c = [c[0]] + [tuple(pywt.threshold(d, T, mode) for d in lvl) for lvl in c[1:]]
    return pywt.waverec2(c, wavelet, mode="periodization")

for mode in ["hard", "soft"]:
    out = np.clip(wavelet_denoise(noisy, mode=mode), 0, 1)
    print(mode, round(psnr(clean, out, data_range=1), 1), "dB")
print("noisy", round(psnr(clean, np.clip(noisy, 0, 1), data_range=1), 1), "dB")

在這個例子中,含雜訊影像的 PSNR 為 20.4 dB,硬門檻 25.7 dB,軟門檻 24.7 dB。通用門檻值本身就很保守(設計目標是幾乎清除所有雜訊),再加上軟門檻的收縮,就會過度平滑;圖 7.5 那種卡通般的外觀就是代價。BayesShrink 這類依子頻帶調整的門檻,或改用不抽樣(平移不變)的轉換,可以同時減輕模糊與方塊假影。scikit-image 在 skimage.restoration.denoise_wavelet 中包裝了這些變體。

左:硬門檻與軟門檻函數圖。右:camera 影像的含雜訊、硬門檻與軟門檻結果裁切,附 PSNR。
圖 7.5 — camera 影像的小波去雜訊(高斯雜訊,在 [0, 1] 範圍內 σ = 0.1),使用 db4、三層分解與通用門檻值。硬門檻保留較銳利的邊緣,但留下零星的亮點;軟門檻較平滑,但也較模糊。

拖動滑桿,比較含雜訊的輸入與軟門檻的結果:

camera 影像含雜訊與小波去雜訊後的裁切比較
含雜訊軟門檻

邊緣偵測

細節子頻帶是三個方向的高通響應,因此小波能當成快速的邊緣偵測器:把近似子頻帶設為零,再做反轉換。剩下的是影像減去它的粗略版本,主要集中在邊緣。只保留部分細節子頻帶就能挑選邊緣方向,如圖 7.6。越粗的層級得到的邊緣越粗,但對雜訊越穩健,這種多尺度行為類似於選擇高斯微分邊緣偵測器的 σ\sigma。

coins 影像、近似子頻帶歸零後的邊緣,以及只保留垂直方向的邊緣。
圖 7.6 — 以 coins 影像兩層哈爾 DWT 得到的邊緣。中:近似子頻帶設為零。右:只保留垂直細節子頻帶,凸顯每枚硬幣左右兩側的邊緣。

壓縮預告

壓縮建立在能量集中之上:轉換、保留少數大係數、再有效率地編碼。圖 7.7 對整張 512×512512 \times 512 影像做四種轉換,只保留最大的 5% 係數。兩種小波比 DCT 高出 2 dB 以上、比 WHT 高出 4 dB 以上,而且誤差的樣貌不同:傅立葉類的基底會在邊緣周圍產生振鈴,並到處損失紋理;小波則保持邊緣清晰,主要在平滑或紋理區域損失細節。

camera 影像,以及只用最大 5% 係數、分別以 DCT、沃爾什–哈達瑪、Haar、db4 重建的結果與 PSNR。
圖 7.7 — 只用最大 5% 係數的重建結果(整張影像轉換,顯示裁切)。PSNR:DCT 28.5 dB、沃爾什–哈達瑪 26.5 dB、Haar 31.0 dB、db4 31.1 dB。

真正的編解碼器比「留下前 5%」細緻得多。JPEG 使用區塊 DCT 加上量化與熵編碼 [4];JPEG 2000 則使用多層二維 DWT 與雙正交濾波器,再對每個子頻帶做位元平面編碼 [15]。第 8 章會詳細介紹這些流程。

現代觀點

小波在 1980 年代末以統一理論之姿進入影像處理,1990 與 2000 年代成為去雜訊與 JPEG 2000 的骨幹,如今則多半是一種成分:存在於稀疏模型、受訊號處理啟發的網路層,以及高效率的生成模型之中。以下是值得一讀的綜論與里程碑論文,以及各自的貢獻。

基礎:Mallat(1989)與 Daubechies(1988)。 Mallat 的論文 [7] 為影像引入多解析度分析,建立正交小波與雙通道濾波器組之間的連結,並提出讓轉換只需 O(N)O(N) 的金字塔演算法。它也提出了分成水平、垂直、對角三個子頻帶的二維可分離分解,至今每個函式庫仍沿用。Daubechies [8] 補上了最後一塊拼圖:具有緊支撐且平滑度可選的正交小波,讓濾波器又短又實用。兩篇合起來,就解釋了 pywt.wavedec2(img, "db4") 為何存在。

教科書等級的綜覽:Vetterli 與 Kovačević(1995)。 《Wavelets and Subband Coding》[2] 從訊號處理的角度建立整套理論:濾波器組、完美重建、時頻鋪磚與編碼。作者買回版權並免費公開全書,是本章濾波器組觀點最好的開放參考資料。

什麼讓小波好用:Unser 與 Blu(2003)。 〈Wavelet theory demystified〉[9] 把任何尺度函數寫成一個 B 樣條(B-spline)與一個殘餘「不規則」部分的卷積,藉此以自成一體的方式重建古典理論。它的主要結論是:實務上在乎的五個性質,也就是逼近階數、多項式重現、消失矩、多尺度微分性質與平滑度,全都由那個 B 樣條因子單獨決定。想知道為什麼 db4 在平滑影像上比 Haar 好,這篇論文給了答案。

從基底到字典:Rubinstein、Bruckstein 與 Elad(2010)。 〈Dictionaries for sparse representation modeling〉[16] 綜覽了從解析式轉換(傅立葉、DCT、小波及其方向性後繼者)到從資料學習的字典這段演變。它的重點是通往現代方法的橋梁:小波「少數係數就能解釋影像」的想法保留了下來,但基底變成可以訓練的東西。這條路線直接通往影像復原中的學習式先驗。

標準的實際樣貌:Skodras、Christopoulos 與 Ebrahimi(2001)。 他們的 JPEG 2000 概論 [15] 說明了 DWT、雙正交濾波器、子頻帶量化與嵌入式編碼如何在真實編解碼器中組合,以及多解析度帶來了哪些功能(解析度可調性、感興趣區域編碼、無失真模式)。

古典去雜訊理論:Donoho 與 Johnstone;Chang、Yu 與 Vetterli。 這幾篇收縮(shrinkage)論文 [12]、[13]、[14] 證明在小波域套一個簡單的非線性函數,對廣泛的訊號類別都接近最佳,而且依子頻帶調整門檻在實務上更有幫助。所有「對轉換係數做門檻來去雜訊」的方法,包括深度學習之前許多手工設計的先驗,都源自這些結果。

深度學習改變了什麼、沒改變什麼

散射網路。 Bruna 與 Mallat 的小波散射轉換(wavelet scattering transform)[17] 把小波卷積、取模(modulus)非線性與平均化層層串接。它是一個濾波器固定為小波、而非經由訓練的卷積網路,具有平移不變性且對形變穩定。第一層的行為類似 SIFT 類描述子,更高層則補充更多資訊,能分辨傅立葉功率頻譜相同的紋理。它在當時的手寫數字與紋理分類上取得最佳結果,並成為解釋 CNN 為何有效的參考點:一個小波框架(如預備知識所述,是冗餘的)加上簡單的非線性,就已經能得到不變性。

把小波當成池化與降取樣層。 跨步卷積與最大池化在降取樣前沒有做適當的低通濾波,因此會混疊。Li 等人在標準分類網路中以 DWT 層取代這些運算(WaveCNet),只保留低頻子頻帶,提升了對雜訊的穩健性 [19]。在影像復原方面,Liu 等人的 MWCNN 以二維 DWT 及其反轉換取代 U-Net 的池化與反池化 [18]。由於 DWT 可逆,網路縮小特徵圖時不會遺失資訊,又能便宜地得到很大的感受野。這種「可逆且能抑制混疊的降取樣」想法,後來出現在許多影像復原與生成模型中。

大感受野。 Finder 等人的 WTConv 層(ECCV 2024)在多層 DWT 的各子頻帶上做小卷積,因此有效感受野隨層數指數成長,而參數量只隨核大小對數成長 [20]。它可以直接放進 ConvNeXt、MobileNetV2 等架構中替換原本的層,作者報告它提升了對影像損壞的穩健性,並使模型更偏重形狀而非紋理。

沒有改變的部分。 DCT 仍然在 JPEG 裡,DCT 及其整數近似也仍是傳統影像與視訊編解碼器的核心工具。手工調校的小波門檻去雜訊,在品質上已被訓練出來的網路超越,但它仍是最快的強基準、不需要訓練資料,而且有理論保證。最重要的是,本章的詞彙(基底、框架、濾波器組、多解析度、混疊、稀疏性)正是用來分析卷積網路在做什麼的詞彙。

重點整理

  • 線性轉換就是換基底。每個係數都是一個內積,也就是影像與某個基底函數的相關值。
  • 正交歸一轉換就是矩陣乘法;可分離的二維情形為 T=AFAT\mathbf{T} = \mathbf{A}\mathbf{F}\mathbf{A}^T,反轉換用轉置即可,且能量守恆。
  • 因為隱含的偶對稱延伸沒有邊界跳躍,DCT 對平滑影像的能量集中優於 DFT 與 WHT,這也是 JPEG 採用它的原因。
  • WHT 與斜轉換分別以簡單性(只有 ±1\pm 1)或明確的斜坡向量換取部分能量集中能力;哈爾則多了局部性。
  • 時頻平面呈現了取捨:傅立葉基底頻率精確、取樣時間精確,小波則依頻率調整磚的形狀。
  • 多解析度分析把兩個短濾波器變成 O(N)O(N) 的快速小波轉換;在二維中,每一層產生一個近似子頻帶與水平、垂直、對角三個細節子頻帶。
  • 小波域的稀疏性帶動了去雜訊(硬/軟門檻)、快速邊緣圖與壓縮(JPEG 2000)。
  • 現代網路把小波當成固定特徵擷取器(散射)、可逆降取樣器(MWCNN、WaveCNet),以及便宜的大感受野層(WTConv)。

練習

  1. 親手做一個框架。 取平面上角度為 90°、210°、330° 的三個單位向量。分別對 f=(1,0)\mathbf{f} = (1, 0) 與 f=(0,1)\mathbf{f} = (0, 1) 計算 ∑k⟨f,sk⟩2\sum_k \langle \mathbf{f}, \mathbf{s}_k \rangle^2。框架界是多少?要如何由三個係數重建 f\mathbf{f}?
提示

兩個總和都等於 3/23/2;事實上對任何單位向量 f\mathbf{f} 總和都是 3/23/2,所以這是 A=B=3/2A = B = 3/2 的緊框架。重建公式為 f=23∑k⟨f,sk⟩sk\mathbf{f} = \frac{2}{3}\sum_k \langle \mathbf{f}, \mathbf{s}_k \rangle \mathbf{s}_k。

  1. 可分離的成本。 對一張 256×256256 \times 256 影像,計算下列情形所需的乘法次數:(a) 把一般二維線性轉換寫成一個 65536×6553665536 \times 65536 矩陣;(b) 以稠密 A\mathbf{A} 做可分離形式 AFAT\mathbf{A}\mathbf{F}\mathbf{A}^T;(c) 三層哈爾 FWT。
提示

(a) N4≈4.3×109N^4 \approx 4.3 \times 10^9。(b) 兩次各 N3N^3 的乘法,2N3≈3.4×1072N^3 \approx 3.4 \times 10^7。(c) 每一層以 2 個係數的濾波器、在半速率下處理目前近似的每一列與每一行;總量是 N2N^2 乘上一個小常數,約數十萬次。

  1. 為什麼要鏡射? 取 f=[1,2,3,4,5,6,7,8]f = [1, 2, 3, 4, 5, 6, 7, 8](一個斜坡),用 SciPy 計算它的 DFT 與 DCT,分別數一數要多少個係數才能涵蓋 99% 的能量,並用隱含的週期延伸與偶對稱延伸解釋差異。
提示

週期延伸會從 8 跳回 1,需要很多 DFT 諧波。偶對稱延伸 …2,1,1,2,…,8,8,7,…\ldots 2, 1, 1, 2, \ldots, 8, 8, 7, \ldots 是連續的,所以 DCT 只需要兩三個係數。

  1. 手算哈爾。 對 [2,2,6,6,3,5,9,9][2, 2, 6, 6, 3, 5, 9, 9] 做兩層哈爾 DWT(平均與差都乘上 1/21/\sqrt{2})。哪些細節係數是零?為什麼?用 pywt.wavedec 驗證你的答案。
提示

第一層的配對是 (2,2)、(6,6)、(3,5)、(9,9),只有 (3,5) 的差不為零。第二層比較相鄰兩對縮放後的和。注意 PyWavelets 的正負號慣例:用它的哈爾濾波器,每個細節係數等於(前 − 後)/2/\sqrt{2},所以 (3, 5) 這一對得到 −2-\sqrt{2}。

  1. 實際調整門檻。 修改去雜訊程式碼,改用 (a) 門檻值 3σ^3\hat{\sigma},以及 (b) 依子頻帶計算的 BayesShrink 門檻 σ^2/σ^X\hat{\sigma}^2 / \hat{\sigma}_X,其中每個子頻帶的 σ^X=max⁡(w2‾−σ^2,0)\hat{\sigma}_X = \sqrt{\max(\overline{w^2} - \hat{\sigma}^2, 0)}。在 σ=0.05\sigma = 0.05 與 σ=0.2\sigma = 0.2 下,和通用門檻值比較 PSNR。
提示

512×512512 \times 512 影像的通用門檻值約為 5σ^5\hat{\sigma},所以 (a) 會保留更多係數。軟門檻時,(b) 應該是三者中最好的。若 σ^X=0\hat{\sigma}_X = 0,就把整個子頻帶設為零。

  1. 對平移的敏感度。 把 camera 影像水平平移一個像素,對兩個版本各做一層哈爾 DWT,比較 V 子頻帶的能量。為什麼平移一個像素會讓係數變化這麼大?不抽樣轉換(pywt.swt2)又會如何表現?
提示

2 倍降取樣讓 DWT 對平移敏感:落在同一個哈爾配對內部的邊緣會產生細節係數,落在兩對之間的同一條邊緣卻不會。平穩(不抽樣)轉換省略降取樣,因此是一個冗餘框架,係數會跟著影像一起平移。

參考文獻

  1. R. C. Gonzalez and R. E. Woods, Digital Image Processing, 4th ed., Pearson, 2018, Ch. 7. publisher page
  2. M. Vetterli and J. Kovačević, Wavelets and Subband Coding, Prentice Hall, 1995(作者提供免費版本). book site
  3. N. Ahmed, T. Natarajan and K. R. Rao, “Discrete Cosine Transform,” IEEE Transactions on Computers, 1974. doi
  4. G. K. Wallace, “The JPEG still picture compression standard,” IEEE Transactions on Consumer Electronics, 1992. doi
  5. W. K. Pratt, J. Kane and H. C. Andrews, “Hadamard transform image coding,” Proceedings of the IEEE, 1969. doi
  6. W. K. Pratt, W.-H. Chen and L. R. Welch, “Slant transform image coding,” IEEE Transactions on Communications, 1974. doi
  7. S. G. Mallat, “A theory for multiresolution signal decomposition: the wavelet representation,” IEEE Transactions on Pattern Analysis and Machine Intelligence, 1989. doi
  8. I. Daubechies, “Orthonormal bases of compactly supported wavelets,” Communications on Pure and Applied Mathematics, 1988. doi
  9. M. Unser and T. Blu, “Wavelet theory demystified,” IEEE Transactions on Signal Processing, 2003. doi
  10. G. R. Lee, R. Gommers, F. Waselewski, K. Wohlfahrt and A. O’Leary, “PyWavelets: A Python package for wavelet analysis,” Journal of Open Source Software, 2019. doi
  11. R. R. Coifman and M. V. Wickerhauser, “Entropy-based algorithms for best basis selection,” IEEE Transactions on Information Theory, 1992. doi
  12. D. L. Donoho and I. M. Johnstone, “Ideal spatial adaptation by wavelet shrinkage,” Biometrika, 1994. doi
  13. D. L. Donoho, “De-noising by soft-thresholding,” IEEE Transactions on Information Theory, 1995. doi
  14. S. G. Chang, B. Yu and M. Vetterli, “Adaptive wavelet thresholding for image denoising and compression,” IEEE Transactions on Image Processing, 2000. doi
  15. A. Skodras, C. Christopoulos and T. Ebrahimi, “The JPEG 2000 still image compression standard,” IEEE Signal Processing Magazine, 2001. doi
  16. R. Rubinstein, A. M. Bruckstein and M. Elad, “Dictionaries for sparse representation modeling,” Proceedings of the IEEE, 2010. doi
  17. J. Bruna and S. Mallat, “Invariant scattering convolution networks,” arXiv:1203.1513, 2012. arXiv
  18. P. Liu, H. Zhang, K. Zhang, L. Lin and W. Zuo, “Multi-level wavelet-CNN for image restoration,” CVPR Workshops (NTIRE), 2018. arXiv
  19. Q. Li, L. Shen, S. Guo and Z. Lai, “Wavelet integrated CNNs for noise-robust image classification,” CVPR, 2020. arXiv
  20. S. E. Finder, R. Amoyal, E. Treister and O. Freifeld, “Wavelet convolutions for large receptive fields,” ECCV, 2024. arXiv