【影像處理筆記】小波變換

灣仔碼農發表於2022-12-17

0 引言

曾經有人問我有關haar的東西，我沒答上來，耿耿於懷，所以我從傅立葉變換學到小波變換再到haar小波，驀然回首，才發現他當時問的是haar特徵。但是，學都學了，不寫筆記很快就會忘記，忘記就等於沒學，只能趕緊記下來。關於小波變換，數字影像處理書上寫的實在難以理解，所以主要參考網上的文章（已在最後註明）。

Haar-like特徵是計算機視覺領域一種常用的特徵描述運算元，也稱為Haar特徵，這是因為Haar-like是受到一維haar小波的啟示而發明的，所以稱為類Haar特徵。

1 從傅立葉變換到小波變換

1.1 傅立葉變換

對於大多數訊號而言，傅立葉分析絕對是非常有用的，因為頻率分析在大多數情況下都非常重要。那麼為什麼我們還需要研究短時傅立葉變換呢（STFT）？原因是因為傅立葉分析有一個非常嚴重的缺點， 在將訊號從時間域變換到頻率域去的時候，把時間資訊丟失了。當我們在用傅立葉變化去分析一個具體訊號的時候，我們不知道哪個頻率是對應在哪個時間點出現的，在哪個時間點消失的。假設有兩個訊號，這兩個訊號都是由sin(t)和sin(5t)組成的，y1是先出現了sin(5t)，再出現了sin(t)，y2是先出現了sin(t)，再出現了sin(5t)。

進行傅立葉變換後的頻譜為

可以看出，變換後的結果是一模一樣的，都在w=1rad/s和w=5rad/s出現了峰值，這就可以說明FT的缺點了——FT只能提供頻域資訊，而完全丟失了時域資訊。不管某一頻率的訊號出現的時間是早還是晚，FT都是將它一視同仁地乘上sin和cos(FT的變換基函式)，然後在整個時間區間加和。因此，它不能提供某一頻率訊號出現的時間。如果一個訊號的頻率並不隨著時間變化，那麼我們稱它為平穩訊號。那麼知道哪一個頻率的訊號在哪一個時間點出現的就不那麼重要了。可是如果現實生活中我們研究的大多數訊號都是非平穩訊號，他們都許多非常短暫變化的特性，這些特點對於我們訊號分析的特點，傅立葉分析並不適合去做這種分析，而短時傅立葉變換則可以。

1.2 短時傅立葉變換

我們將訊號從中間截斷，左邊進行一次FT，右邊進行一次FT，便可以得到，在y1中(0, 25)的訊號是5rad/s的頻率，(25, 50)的訊號是5rad/s的頻率，y2恰好相反。這就是短時傅立葉變換的基本原理。截斷的方法在專業上將叫分窗——有一個窗子在訊號上從左向右滑動，每次只能看到訊號的一部分，相當於把整個時域過程分解成無數個等長的小過程，再傅立葉變換，就知道在哪個時間點上出現了什麼頻率了。

短時傅立葉變換輸出為三維圖形，分別為時間、頻率、強度三個軸（顏色即為對應時間、頻率下的訊號強度），時間頻率軸上可明顯觀察到該訊號的頻率成分，隨時間逐漸由20Hz線性增加到100Hz。

定義方窗函式為ywindow，窗長為width，將方窗函式向右平移了短時傅立葉變換可以表示為

這裡之所以

線上性代數中，基是描述、刻畫向量空間的基本工具。向量空間的基是它的一個特殊的子集，基的元素稱為基向量。向量空間中任意一個元素，都可以唯一地表示成基向量的線性組合。比如，向量空間的一個向量可以分解在x，y方向，同時在各個方向定義單位向量e1、e2，這樣任意一個向量都可以表示為a=xe1+ye2，這是二維空間的基。

在變換中，我們將原始訊號乘上的訊號稱為基函式。在傅立葉變換中，e-jwt是這個變換的基函式。STFT將基函式乘上一個方窗函式，形成了一個新的基函式：

假設用正弦函式sin(5t)表示原來的基函式

FT的基函式是在時域無限延伸的，因此，無論怎麼平移，都是任分佈在整個時域的，起不到分窗的作用。而STFT的基函式只在時域一段不為0，在剩下的時域都是0，因此，STFT的基函式的平移，就相當於自動加了窗子。STFT的本質是基函式的改變。

對於低頻訊號，為了更好地確定頻率，我們希望，時域區間寬一些，即時間不確定度
對於高頻訊號，為了更好地在時域定位，我們希望，時域區間窄一些，即時間不確定度

所以，對於時變非穩態訊號，高頻部分適合用小窗(短週期)，低頻部分適合用大窗(長週期)。然而，在一次STFT中，視窗的寬度是固定，即時域解析度是固定的，根據海森堡測不準原理，其頻域解析度也是固定的。也就是說，不論高頻低頻，其時域和頻域解析度都不可調，這是STFT的缺點。

1.3 連續小波變換

小波分析在時間域和頻率域同時具有良好的區域性化性質，即在低頻部分具有較高的頻率解析度和較低的時間解析度，在高頻部分具有較高的時間解析度和較低的頻率解析度，因此有“顯微鏡”之稱。小波變換具有對訊號的自適應能力，因而比傅立葉分析更適合處理非平穩訊號。小波做的改變就在於，將無限長的三角函式基換成了有限長的會衰減的小波基。

從公式可以看出，不同於傅立葉變換，變數只有頻率ω，小波變換有兩個變數：尺度a(scale)和平移量τ(translation)。尺度a控制小波函式的伸縮，平移量τ控制小波函式的平移。尺度就對應於頻率（反比），平移量τ就對應於時間。比如下圖，中間的小波s較小，相當於擠壓；右邊的s較大，相當於拉伸。

根據上面講的，滑動相當於分窗，窗的長度是基函式不為零的長度。中間的圖，

Step2：把小波向右移，距離為τ，得到小波ψ(t-τ)，重複步驟1。再把小波右移，得到小波ψ(t-2τ)，重複步驟1。按上述步驟一直重複下去，知道訊號f(t)結束。

Step3：擴充小波ψ(t)，例如擴充套件一倍，得到的小波函式ψ(t/s)，重複步驟1-2。

Step4：重複步驟1-3。

接下來我們對一個訊號就行一次連續小波變換（CWT）。下圖中藍色部分為小波函式，黃色部分為訊號。

如上圖，選擇較小的此時，窗子較窄（小波函式不為0的部分窄），時間解析度好，頻率解析度差。

如上圖，將此時，窗子變寬了（小波函式不為0的部分變寬），時間解析度變差，頻率解析度變好。

如上圖，將s 進一步增大，對小波母函式進行縮放，此時小波函式頻率再次降低，窗子更寬（小波函式不為0的部分更寬），用來篩選低頻部分。小波函式在時間軸上平移，每一次平移就先相乘，再積分，篩選出訊號中與自己頻率相近的部分。此時，窗子很寬（小波函式不為0的部分很寬），時間解析度差，頻率解析度很好。

1.4 總結

FT的基函式，是分佈在(−∞,+∞)
STFT的基函式，是用窗函式截斷的sin,cos （圖中是被高斯窗截斷的），具有了緊支撐性，時域平移等同於分窗，使得STFT既能篩選頻率，也能篩選時間。但是STFT基函式是：先確定頻率，再與窗函式相乘構成的。因此不同的頻率，具有同樣的時間和頻率解析度。另外，窗函式的長短也比較難以確定。
CWT的基函式，是小波函式，具有緊支撐性，時域平移等同於分窗，使得CWT既能篩選頻率，也能篩選時間。小波函式在改變頻率的時候，是透過“縮放”實現的，這使得小波函式在改變頻率的同時，改變了窗長。因此不同的頻率，具有不同的時間和頻率解析度，實現了解析度動態可調。

2 從小波變換到haar小波

2.1 離散小波變換

離散小波變換是對基本小波的尺度和平移進行離散化。為了解決計算量問題，縮放因子和平移引數都選擇2j的倍數。執行離散小波變換的有效方法是使用濾波器。該方法時Mallat再1988年開發的，叫做Mallat演算法。這種方法實際上是一種訊號的分解方法，在數字訊號處理中稱為雙通道子帶編碼。

首先將原始訊號作為輸入訊號，透過一組正交的小波基分解成高頻部分和低頻部分，A表示訊號的近似值（approximations），是大的縮放因子產生的係數，表示訊號的低頻分量。D表示訊號的細節值(detail)，是小的縮放因子產生的係數，表示訊號的高頻分量。然後將得到的低頻部分作為輸入訊號，又進行小波分解，得到下一級的高頻部分和低頻部分，以此類推。隨著小波分解的級數增加，其在頻域上的解析度就越高。這就是多解析度分析（MRA，MultiResolution Analysis）。

離散小波變換可以被表示成由低通濾波器和高通濾波器組成的一棵樹。原始訊號透過這樣的一對濾波器進行的分解叫做一級分解。訊號的分解過程可以進行多級分解，分解級數的多少取決於要被分析的資料和使用者的需要。小波分解樹只對訊號的低頻分量進行連續分解。把分解的係數還原成原始訊號的過程叫做小波重構，數學上叫做逆離散小波變換。

任何小波變換的基函式，其實就是對母小波和父小波縮放和平移的集合。縮放倍數都是2的級數，平移的大小和當前其縮放的程度有關。小波系統有很多種，不同的母小波，衍生的小波基就完全不同。小波展開的近似形式是這樣的：

其中的ψ_j,k(t)就是小波級數，這些級數的組合就形成了小波變換中的基。

2.2 haar小波函式

最簡單的基函式是haar基函式，是由一組分段常值函式組成的函式集合，其父小波，也叫尺度函式（scaling function）

其中，i=0,1,...,(2^j-1)。j是尺度因子，改變j使函式圖形縮小或者放大；i為平移引數，改變i使函式沿x軸方向平移。

尺度函式可以代表一個空間，可以看出，假如有一個函式他屬於一個某空間，那你將其在時域上平移，它還是屬於這個空間。但如果你對它頻域的放大或縮小，它就會相應移到下一個或者上一個空間了。在向量空間Vj中的每一個向量被包含在向量空間Vj+1中。我個人的理解，Φ₀⁰可以由Φ₀¹和Φ₁¹表示，Φ₀¹可以由Φ₀²和Φ₁²表示等等，即任何尺度平移的尺度函式，都可以用更加精細的尺度層面上的尺度函式構建出來。j越大，尺度函式所做的時域上的平移幅度會越小，相應的在j子空間裡面得到的f(t)表示粒度會很細，細節展現很多。不同的子空間有不同的解析度，這就是用不同的解析度去看目標訊號。

母小波

其中，i=0,1,...,(2^j-1)。j是尺度因子，改變j使函式圖形縮小或者放大；i為平移引數，改變i使函式沿x軸方向平移。

向量空間Wj中的每一個向量也被包含在向量空間Wj+1中。下面我們對一個原始訊號進行多解析度分析：

可以看出，在不同的子空間，對於同一個訊號就有不同的詮釋。詮釋最好的是這就是多解析度的意義。我們可以有巢狀的，由尺度函式演變的基函式集合，每一個集合都提供對原始訊號的某種近似，解析度越高，近似越精確。

從上面也可以看出，直接用尺度函式就可以，或者母小波就可以還原訊號。為什麼要用兩個呢？

回顧上面的向量空間圖，向量空間W0中的小波函式，和V0中的尺度函式，可以構建V1中的尺度函式。向量空間W1中的小波函式，和V1中的尺度函式，可以構建V2中的尺度函式。繼續推導下去，可以發現，尺度j下的小波函式，可以用來將Vj的基擴充到Vj+1中。這代表著，對任何一個子空間Vj，我們現在有兩種方法去得到它的正交基：（1）用它原本的基Φ_j,k（2）用上一級子空間的Φ_j-1,k和上一級子空間的小波函式ψ_j-1,k。第二種選擇能給我們帶來額外的好處，那就是我們可以迴圈不斷地用上一級子空間的尺度函式以及小波函式的組合來作為當前子空間的基。換句話說，如果針對V3這個子空間，它實際上就有四種不同的，但是等價的正交基。

下面是V3子空間的第二種可選擇的正交基，V2的尺度函式和W2的小波函式：（上面沒歸一化，這邊歸一化了）

左邊這四個基和原始訊號作內積，表徵的是訊號的平均。而右邊的這四個基和原始訊號作內積則表徵了在平均中丟失的訊號細節。得益於此，多解析度分析能夠對訊號在越來越寬的區域上取平均，等同於做低通濾波，而且，它還能保留因為平均而損失的訊號細節，等同於做高通濾波。也就是說，尺度函式和小波函式背後的物理意義是：小波函式等同於對訊號做高通濾波保留變化細節，而尺度函式等同於對訊號做低通濾波保留平滑的形狀。

所有訊號空間中的訊號都可以寫成組成這個基的函式的線性組合：

對應的係數的計算和平常一樣：

小波變換的基礎流程：

1. 選取合適的小波函式和尺度函式，從已有的訊號中，反算出係數c和d。

2. 對係數做對應處理

3. 從處理後的係數中重新構建訊號。

這裡的係數處理是區別你的應用的重點。比如影像或者影片壓縮，就希望選取能將能量聚集到很小一部分系數中的小波，然後拋棄那些能量很小的小波係數，只保留少數的這些大頭係數，再反變換回去。這樣的話，影像訊號的能量並沒有怎麼丟失，影像體積卻大大減小了。

3 應用

3.1 一維harr小波

示例：求只有4個畫素[9 7 3 5]的影像的哈爾小波變換系數

步驟1：求均值(averaging)，也叫Approximation。計算相鄰畫素對的平均值，得到一幅解析度比較低的新影像，新的影像的解析度是原來的1/2，相應的畫素值為：[8 4]

步驟2：求差值(differencing)。很明顯，用2個畫素表示這幅影像時，影像的資訊已經部分丟失。為了能夠從由2個畫素組成的影像重構出由4個畫素組成的原始影像，就需要儲存一些影像的細節係數(detail coefficient)，以便在重構時找回丟失的資訊。方法是使用這個畫素對的差值除以2，(9-7)/2=1，(3-5)/2=-1，結果為[1 -1]。到此為止，小波變換的一級分解結束，結果為[8 4 1 -1]。下面是小波變換的二級分解，即對低頻部分[8,4]進行分解。

步驟3：重複第1，2步，把由第一步分解得到的影像進一步分解成解析度更低的影像和細節係數，結果為[6 2 1 -1]。

在變換過程沒有丟失資訊，因為能夠從所記錄的資料中重構出原始影像。首先在解析度為1的影像上重構出解析度為2的影像，在解析度為2的影像上重構出解析度為4的影像。透過變換之後產生的細節係數的幅度值比較小，這就為影像壓縮提供了一種途徑，例如去掉一些微不足道的細節係數並不影響對重構影像的理解。上面的分解其實就是之前說的，可以在不同的向量空間用不同的基表示該影像，用不同的解析度去看影像

3.2 二維harr小波分解與重構

以離散餘弦變換為基礎的JPEG標準演算法，把圖片分塊處理會產生塊效應，在小波變換中，由於小波變換中適用的基函式的長度是可變的，因此無需把輸入影像進行分塊，避免了塊效應。下圖為進行多級小波變換後，影像各部分存放的是不同級數的小波係數。左上角的元素表示整個影像塊的畫素值的平均值，其餘是該影像塊的細節係數。如果從矩陣中去掉影像的某些細節係數，而且重構的影像質量仍然可以接受，則可實現壓縮。具體做法是設定一個閾值，當細節係數小於這個閾值，就當做0來看待。

示例小波影像分解→閾值處理→重構

# include<opencv.hpp>
# include<iostream>
using namespace std;
using namespace cv;
Mat waveletTransform2D(Mat& img, int level) {
    int width = img.cols;
    int height = img.rows;
    int depth = level;
    int depthcount = 1;
    Mat wavelet = Mat::zeros(img.size(), CV_32FC1);        
    Mat tmp = Mat::zeros(img.size(), CV_32FC1);
    Mat imgtmp = img.clone();
    imshow("src", imgtmp);
    imgtmp.convertTo(imgtmp, CV_32FC1, 1.0);

    while (depthcount <= depth)
    {
        height = img.rows / depthcount;
        width = img.cols / depthcount;

        // 水平
        for (int i = 0; i < height; i++)
        {
            for (int j = 0; j < width / 2; j++)
            {
                tmp.at<float>(i, j) = (imgtmp.at<float>(i, 2 * j) + imgtmp.at<float>(i, 2 * j + 1)) / 2; //水平方向求均值
                tmp.at<float>(i, j + width / 2) = (imgtmp.at<float>(i, 2 * j) - imgtmp.at<float>(i, 2 * j + 1)) / 2; //水平方向求差值
            }
        }

        // 豎直
        for (int i = 0; i < height / 2; i++)
        {
            for (int j = 0; j < width; j++)
            {
                wavelet.at<float>(i, j) = (tmp.at<float>(2 * i, j) + tmp.at<float>(2 * i + 1, j)) / 2; //得到總均值和水平方向差值
                wavelet.at<float>(i + height / 2, j) = (tmp.at<float>(2 * i, j) - tmp.at<float>(2 * i + 1, j)) / 2;//得到豎直方向均值和對角差值
            }
        }
        imgtmp = wavelet.clone();
        depthcount++;
    }
    return wavelet;
}
Mat invWaveletTransform2D(Mat& img, int level) {
    int width = img.cols;
    int height = img.rows;
    int depth = level;
    Mat tmp = Mat::ones(img.size(), CV_32FC1);
    Mat wavelet = img.clone();
    Mat imgtmp = img.clone();
  
    while (depth > 0)
    {
        height = img.rows / depth;
        width = img.cols / depth;

        // 豎直
        for (int i = 0; i < height; i+=2)
        {
            for (int j = 0; j < width; j++)
            {
                tmp.at<float>(i, j) = (imgtmp.at<float>(i/2, j) + imgtmp.at<float>(i/ 2+height / 2, j)); 
                tmp.at<float>(i + 1, j) = (imgtmp.at<float>(i/2, j) - imgtmp.at<float>(i/ 2 + height / 2, j)) ;
            }
        }

        // 水平
        for (int i = 0; i < height; i++)
        {
            for (int j = 0; j < width; j+=2)
            {
                wavelet.at<float>(i, j) = (tmp.at<float>(i, j/2) + tmp.at<float>(i, j/2+width/2)); 
                wavelet.at<float>(i, j + 1) = (tmp.at<float>(i, j/2) - tmp.at<float>(i, j/2 + width / 2)); 
            }
        }
        imgtmp = wavelet;
        depth--;
    }
    return wavelet;
}
int main()
{
    Mat img = imread("./2.tif", cv::IMREAD_GRAYSCALE);
    //Mat img = (Mat_<uchar>(4, 4)<<1, 2, 3, 4,5,6,7,8,9,10,11,12,13,14,15,16);
    
    // 1.小波變換
    Mat wavelet = waveletTransform2D(img, 2);
    // 2.閾值處理
    for (int i = 0; i < wavelet.rows; i++)
    {
        for (int j = 0; j < wavelet.cols; j ++)
        {
            if (abs(wavelet.at<float>(i, j)) < 15)
                wavelet.at<float>(i, j) = 0;
        }
    }
    // 3.逆向小波變換
    Mat dst = invWaveletTransform2D(wavelet, 2);
    dst.convertTo(dst, CV_8UC1);
    imshow("dst", dst);
    waitKey(0);
    return 0;
}

一般，有用訊號表現為低頻訊號或者是一些比較平穩的訊號，而噪聲訊號則通常表現為高頻訊號，因此根據噪聲與訊號在不同尺度（即不同頻率）上的小波譜具有不同表現的特點，將噪聲小波譜占主導地位的那些尺度上的小波分量去掉，而實現去噪。根據部落格基於小波變換實現影像增強：

小波分解後刪除高頻，即低通濾波，可以讓影像變得平滑,濾除影像中的噪聲。
小波分解後刪除低頻，即高通濾波，可以提取影像邊緣，進而銳化影像
小波分解後弱化高頻，增強低頻，可以增強對比度

該部落格用matlab中的wavedec2函式 [c,s]=wavedec2(X,N,’wname’)，c為各層分解係數，s為各層分解係數長度。c是一個行向量，c=[A(N)|H(N)|V(N)|D(N)|H(N-1)|V(N-1)|D(N-1)|H(N-2)|V(N-2)|D(N-2)|…|H(1)|V(1)|D(1)]，A(N)代表第N層低頻係數，H(N)|V(N)|D(N)代表第N層高頻係數，分別是水平，垂直，對角高頻，以此類推,到H(1)|V(1)|D(1)。比如，設一級高頻係數為0，也就是c的最後部分H(1)|V(1)|D(1)置為0，對應c++這邊，就是將分解得到的影像右半邊和下半邊都置為0。

參考：

1. 岡薩雷斯《數字影像處理（第四版）》Chapter7（所有圖片可在連結中下載）

2. 短時傅立葉變換（Short Time Fourier Transform）

3. 訊號處理基礎——傅立葉變換與短時傅立葉變換

4. 小波變換（深入淺出）

5. 從傅立葉變換進階到小波變換（三）

6. 影像處理-小波變換

7. 哈爾(Haar)小波變換的原理及opencv原始碼

8. 小波變換完美通俗講解系列之（二）

9. 小波變換（wavelet transform）的通俗解釋

10. 基於小波變換實現影像增強