Part1-Introduction To The Wavelet Transform（簡介）

1、Origin of the wavelet transform:
The theories of Wavelet originate from diffierent areas of study:

• Engineering

• Time-frequency analysis and Multiresolution Analysis

• Computer Vision

• Pyramidal algorithm

• Physics

• Pure Mathematics

2、Multiresolution Analysis（MRA-多解析度分析與處理）

• Multiresolution analysis is about analyzing a signal based on the information appeared in different scales of such signal – to mimic human beings in analyzing signals.

• A signal of a certain “scale” refers to its best approximation at a certain resolution

• By “traveling” from the coarse scales toward the fine scales, one zooms in and arrives at a more exact representation of the given signal
  
  3、The Discrete Wavelet Transform(離散小波變換)
  A simple way to implement MAR is by using the Discrete Wavelet Transform(DWT):
  
  The definition of wavelet transform shows that the wavelet analysis is a measure of similarity the basis functions(wavelets) and the original function .The coefficients ，named H0 for Low Pass Filter

  , and H1 for High Pass Filter, caculated indicate how close the function is to the danghter wavelet at that particular scale.
  

Part2-Decomposition(DWT) and Reconstruction(Inverse DWT)–離散小波分解與重構

對於離散小波變換，由於很多小波函式不是正交函式，因此需要一個尺度係數（Scaling Coefficients）和一個小波係數（Wavelet Coefficients）.因此，原訊號函式可以分解成尺度函式（係數）和小波函式（係數）的線性組合，在這個函式中，尺度函式產生低頻部分，小波函式產生高頻部分。

1、一維Haar小波變換-分解

離散的訊號經過Haar小波變換時，首先會將這個訊號所攜帶的資訊進行壓縮，得到N/2個數據點進行儲存，那麼這些點資訊就是由這些Scaling Coefficients-尺度係數來表徵的。對於離散訊號而言，其可能具有高頻和低頻成分，而小波係數或細節係數代表它的高頻部分。

小波中的下采樣就是對訊號進行隔點取樣，目的就是為了將資訊進行壓縮儲存。
小波中的上取樣就是隔點插零，目的是為了重構訊號。

而對訊號的濾波過程，在數學上等效為訊號與濾波器衝激響應的卷積。

分解LP 濾波器（離散卷積運算元）：

分解HP濾波器（離散卷積運算元）：

離散時間卷積定理：

“離散卷積”是兩個離散序列x(n) 和h(n) 之間按照一定的規則將它們的有關序列值分別兩兩相乘再相加的一種特殊的運算。具體可用公式表示為:

其中 y（n）就是經過卷積運算以後所得到的一個新的序列。

在工程上離散卷積有著廣泛的應用，例如，在數字影象處理領域，為了將數字訊號進行濾波，可以將表示成離散序列的影象訊號C(n) 與數字濾波器的衝激響應h(n) 進行離散卷積執行。

多分辨分析的小波函式Ψ(t)和尺度函式φ(t)滿足雙尺度的差分方程：

多分辨分析每一層分解使訊號f(t)通過一低通濾波器和帶通濾波器，把訊號分解為低頻部分和高頻部分。低通濾波器的特性由小波函式Ψ(x)確定，帶通濾波器的特性由尺度函式φ(x)確定。分解後的係數由兩部分組成：低頻係數向量c1和高頻係數向量d1。低頻係數向量c1 是由訊號與低通濾波器（小波函式確定）的脈衝響應經過卷積運算得到，高頻係數向量d1是由訊號與帶通濾波器（尺度函式確定）經過卷積運算得到。

2、一維Haar小波變換–重構

逆變換過程：

3、二維Haar小波變換

二維影象訊號

對於二維影象訊號，可以用分別在水平和垂直方向進行濾波的方法實現二維小波多解析度分解。圖2.5為經過二維離散小波變換的分解後子影象的劃分。其中:
(l)LL子帶是由兩個方向利用低通小波濾波器卷積後產生的小波係數，它是影象的近似表示。
(2)HL子帶是在行方向利用低通小波濾波器卷積後，再用高通小波濾波器在列方向卷積而產生的小波係數，它表示影象的水平方向奇異特性。（水平子帶）
(3)LH子帶是在行方向利用高通小波濾波器卷積後，再用低通小波濾波器在列方向卷積而產生的小波係數，它表示影象的垂直方向奇異特性。（垂直子帶）
(4)HH子帶是由兩個方向利用高通小波濾波器卷積後產生的小波係數，它表示影象的對角邊緣特性。（對角子帶）
第一個字母表示列方向的處理，第二個字母表示行方向的處理，影象的奇異特性通過低通時保留，通過高通時被濾除。

小波去噪方法也就是尋找從實際訊號空間到小波函式空間的最佳映像，以便得到原訊號的最佳恢復。

目前，小波去噪的方法大概可以分為三大類:

第一類方法–小波變換模極大值去噪法

利用小波變換模極大值原理去噪，即根據訊號和噪聲在小波變換各尺度上的不同傳播特性，剔除由噪聲產生的模極大值點，保留訊號所對應的模極大值點，然後利用所餘模極大值點重構小波係數，進而恢復訊號;

第二類方法–小波係數相關性去噪法

對含噪訊號作小波變換之後，計算相鄰尺度間小波係數的相關性，根據相關性的大小區別小波係數的型別，從而進行取捨，然後直接重構訊號;

第三類方法–小波變換閾值去造法

小波閾值去噪方法，該方法認為訊號對應的小波係數包含有訊號的重要資訊，其幅值較大，但數目較少，而噪聲對應的小波係數是一致分佈的，個數較多，但幅值小。

4、小波閥值收縮去噪法:

1、 小波閥值去噪的基本思想:

Donoho提出的小波閥值去噪的基本思想是將訊號通過小波變換（採用Mallat演算法）後，訊號產生的小波係數含有訊號的重要資訊，將訊號經小波分解後小波係數較大，噪聲的小波係數較小，並且噪聲的小波係數要小於訊號的小波係數，通過選取一個合適的閥值，大於閥值的小波係數被認為是有訊號產生的，應予以保留，小於閥值的則認為是噪聲產生的，置為零從而達到去噪的目的。其基本步驟為：
（1）分解：選定一種層數為N的小波對訊號進行小波分解；
（2）閥值處理過程：分解後通過選取一合適的閥值，用閥值函式對各層係數進行量化；
（3）重構：用處理後的係數重構訊號。

2、小波閥值去噪的基本問題

小波閥值去噪的基本問題包括三個方面：小波基的選擇，閥值的選擇，閥值函式的選擇。
（1）小波基的選擇：通常我們希望所選取的小波滿足以下條件：正交性、高消失矩、緊支性、對稱性或反對稱性。但事實上具有上述性質的小波是不可能存在的，因為小波是對稱或反對稱的只有Haar小波，並且高消失矩與緊支性是一對矛盾，所以在應用的時候一般選取具有緊支的小波以及根據訊號的特徵來選取較為合適的小波。
（2）閥值的選擇：直接影響去噪效果的一個重要因素就是閥值的選取，不同的閥值選取將有不同的去噪效果。目前主要有通用閥值（VisuShrink）、SureShrink閥值、Minimax閥值、BayesShrink閥值等。
（3）閥值函式的選擇：閥值函式是修正小波係數的規則，不同的反之函式體現了不同的處理小波係數的策略。最常用的閥值函式有兩種：一種是硬閥值函式，另一種是軟閥值函式。還有一種介於軟、硬閥值函式之間的Garrote函式。
另外，對於去噪效果好壞的評價，常用訊號的信噪比（SNR）與估計訊號同原始訊號的均方根誤差（MSE）來判斷。

例項：基於通用閾值和SURE閾值的小波去噪方法的MATLAB實現

%
% 2d wavelet shrinkage image denoising with universal threshold
%
close all

% Load image
imagefilename = 'Lena512Gray.bmp';
thewavelet = 'db4';
noise_var = 300;

x = double(imread(imagefilename));
n = sqrt(noise_var)*randn(size(x));
xn = x + n;

% Make sure the additive noise will not make the signal component to have
% magnitude larger than 255.
for i=1:size(x,1)
    for j=1:size(x,2)
        if (xn(i,j) > 255) 
            xn(i,j) = 255;
        end
    end
end

std_noise = sqrt(sum(sum((xn-x).^2))/length(x)^2);

figure
imagesc([x,xn]);
axis image
colormap(gray);
title('original, noisy');


% 2d wavelet transform
% 1 levels decomposition with periodic extension.
dwtmode('per');
[ca1,ch1,cv1,cd1] = dwt2(x,  thewavelet);
[na1,nh1,nv1,nd1] = dwt2(xn, thewavelet);
[na2,nh2,nv2,nd2] = dwt2(na1, thewavelet);

% Display
figure
imagesc([ca1,ch1;cv1,cd1]);
axis image
colormap(gray);
title('1st level decomposition with periodic extension.');
figure
imagesc([na1,nh1;nv1,nd1]);
axis image
colormap(gray);
title('1st level decomposition with periodic extension of noisy image.');


% First normalize the wavelet coefficients such that the noise component
% has a variance of 1 
nor_nh1 = nh1/std_noise;
nor_nv1 = nv1/std_noise;
nor_nd1 = nd1/std_noise;

% Find the threshold!
disp('Universal threshold used for level 1')
thr_uni = zeros(1,3);
thr_uni(:) = sqrt(2*log(size(x,1)*size(x,2)));

% Shrinkage (SURE)

thr_sure = zeros(1,3);
thr_sure(1) = thselect(nor_nh1,'rigrsure');
thr_sure(2) = thselect(nor_nv1,'rigrsure');
thr_sure(3) = thselect(nor_nd1,'rigrsure');

disp('SURE threshold used for level 1');

% Shrinkage (universal threshold)
dnh1_uni = wthresh(nor_nh1,'s',thr_uni(1));
dnv1_uni = wthresh(nor_nv1,'s',thr_uni(2));
dnd1_uni = wthresh(nor_nd1,'s',thr_uni(3));
% Shrinkage (SURE)
dnh1_sure = wthresh(nor_nh1,'s',thr_sure(1));
dnv1_sure = wthresh(nor_nv1,'s',thr_sure(2));
dnd1_sure = wthresh(nor_nd1,'s',thr_sure(3));

% Denormalize the denoised wavelet coefficients
dnh1_uni = dnh1_uni*std_noise;
dnv1_uni = dnv1_uni*std_noise;
dnd1_uni = dnd1_uni*std_noise;

dnh1_sure = dnh1_sure*std_noise;
dnv1_sure = dnv1_sure*std_noise;
dnd1_sure = dnd1_sure*std_noise;



% Reconstruction (universal threshold)
rx = idwt2(ca1,ch1,cv1,cd1,thewavelet,size(x));
rdnx_uni = idwt2(na1,dnh1_uni,dnv1_uni,dnd1_uni,thewavelet,size(x));

% Reconstruction (SURE)
rdnx_sure = idwt2(na1,dnh1_sure,dnv1_sure,dnd1_sure,thewavelet,size(x));


figure
imagesc([x,xn,rdnx_uni,rdnx_sure]);
axis image
colormap(gray);
title('original, noisy, denoised(universal), denoised(SURE)');

res = sprintf('SNR of noisy image: %f dB', snr(x, xn));
res1 = sprintf('SNR of denoised image (universal threshold): %f dB', snr(x, rdnx_uni));
res2 = sprintf('SNR of denoised image (SURE): %f dB', snr(x, rdnx_sure));
disp(res);
disp(res1);
disp(res2);

Fig.1 Lena 原圖

Fig .2 加入噪聲之後的影象

計算各影象的SNR（信噪比）得到如下資料：
SNR of noisy image: 38.468755 dB
SNR of denoised image (universal threshold): 40.259733 dB
SNR of denoised image (SURE): 41.010018 dB

結論：小波變換在影象去噪等方面有著廣泛應用，實現小波閾值收縮去噪的方法主要有通用閾值、軟硬閾值、SURE閾值等。1階SURE演算法去噪較通用閾值演算法去噪，SNR提高了大約1dB，去噪效果顯著。