小波變換

  傳統的訊號理論，是建立在Fourier分析基礎上的，而Fourier變換作為一種全域性性的變化，其有一定的侷限性，如不具備區域性化分析能力、不能分析非平穩訊號等。在實際應用中人們開始對Fourier變換進行各種改進，以改善這種侷限性，如STFT（短時傅立葉變換）。由於STFT採用的的滑動窗函式一經選定就固定不變，故決定了其時頻解析度固定不變，不具備自適應能力。

  小波分析很好的解決了這個問題。小波直接把傅立葉變換的基給換了——將無限長的三角函式基換成了有限長的會衰減的小波基。這樣不僅能夠獲取頻率，還可以定位到時間了。

  常見小波函式：Haar、Daubechies、Biorthogonal、Coiflets、Symlets、Morlet、Mexican Hat、Meyer、Gaus、Dmeyer、ReverseBior、Cgau、Cmor、Fbsp、Shan.

程式碼示例

function y=wtfusion(x1,x2,N,wname)

%函式功能：
%     函式x=wtfusion(x1,x2,N,wname)將兩幅原影象x1,x2進行基於小波變換的影象融合，得到融合後的影象y
%     近似分量採用加權平均的融合規則，各細節分量採用基於區域特性量測的融合規則
%輸入引數：
%     x1----輸入原影象1
%     x2----輸入原影象2
%     N----小波分解的層數
%     wname----小波基函式
%輸出引數：
%     y----原影象融合後得到的影象
%-----------------------------------------------------------------%
 


x1=double(x1);                   %將uint8的影象資料型別轉換成double型進行資料處理
x2=double(x2);

 %將原影象x1,x2分別進行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層高頻係數,分別是水平,垂直,對角高頻
 %S(N+2行2列)的結構是儲存各層分解係數長度的,即第一行是A(N)的長度（其實是A(N)的原矩陣的行數和列數）,
 

 %第二行是H(N)|V(N)|D(N)|的長度,第三行是H(N-1)|V(N-1)|D(N-1)的長度,
 %倒數第二行是H(1)|V(1)|D(1)長度,最後一行是X的長度(大小)

[C1,S1]=wavedec2(x1,N,wname); 
[C2,S2]=wavedec2(x2,N,wname);  

A1=appcoef2(C1,S1,wname,N);            %提取出小波分解的近似分量
A2=appcoef2(C2,S2,wname,N);
A=0.5*A1+0.5*A2;                       %近似分量的融合規則採用加權平均的方法

%仿照matlab中近似分量和細節分量的儲存方式，把融合後的近似分量和細節分量轉成行向量，然後存入向量C中
%這樣做是為了方便重構原影象

a=reshape(A,1,S1(1,1)*S1(1,2));        %將A轉換成行向量
C=a;

for i=N:-1:1                           %迴圈從第N層到第1層    
    [H1,V1,D1]=detcoef2('all',C1,S1,i);    %提取出小波分解的各層細節分量
    [H2,V2,D2]=detcoef2('all',C2,S2,i);
    H=f(H1,H2);                            %各層細節分量的融合規則採用基於區域特性量測的融合規則
    V=f(V1,V2);
    D=f(D1,D2);
    h=reshape(H,1,S1(N+2-i,1)*S1(N+2-i,2));%分別將融合後的細節分量轉成行向量，並存入行向量C中
    v=reshape(V,1,S1(N+2-i,1)*S1(N+2-i,2));
    d=reshape(D,1,S1(N+2-i,1)*S1(N+2-i,2));
    C=[C,h,v,d];
end

S=S1;
y=waverec2(C,S,wname);      %重構原影象
figure(1);imshow(uint8(y));title('基於小波變換的融合影象')
end

function y=f(x1,x2)

%函式功能：
%       y=f(x1,x2)將兩幅原影象x1和x2基於區域特性量測的融合規則進行融合,得到融合後的影象y
%       首先計算兩幅影象的匹配度，若匹配度大於閾值，說明兩幅影象對應區域性能量較接近，
%       因此採用加權平均的融合方法；若匹配度小於閾值，說明兩幅影象對應區域性能量相差較大，
%       因此選取區域性區域能量較大的小波係數作為融合影象的小波係數
%輸入引數：
%      x1----輸入原影象1
%      x2----輸入原影象2
%輸出引數：
%      y----融合後的影象
%------------------------------------------------------------%

w=1/16*[1 2 1;2 4 2;1 2 1];   %權係數
E1=conv2(x1.^2,w,'same');     %分別計算兩幅影象相應分解層上對應區域性區域的“能量”
E2=conv2(x2.^2,w,'same');
M=2*conv2(x1.*x2,w,'same')./(E1+E2);%計算兩幅影象對應區域性區域的匹配度
T=0.7;                              %定義匹配閾值
Wmin=1/2-1/2*((1-M)/(1-T));
Wmax=1-Wmin;
[m,n]=size(M);

for i=1:m
    for j=1:n
        if M(i,j)<T                %如果匹配度小於匹配閾值，說明兩幅影象對應區域性區域能量距離較遠；
            if E1(i,j)>=E2(i,j)    %那麼就直接選取區域能量較大的小波係數
                y(i,j)=x1(i,j);
            else
                y(i,j)=x2(i,j);
            end
        else                       %如果匹配度大於匹配閾值，說明兩幅影象對應區域性區域能量比較接近；
            if E1(i,j)>=E2(i,j)    %那麼就採用加權的融合演算法
                y(i,j)=Wmax(i,j)*x1(i,j)+Wmin(i,j)*x2(i,j);
            else
                y(i,j)=Wmin(i,j)*x1(i,j)+Wmax(i,j)*x2(i,j);
            end
        end
    end
end

end

函式呼叫

x1 = imread('lena1.jpg');
x2 = imread('lena2.jpg');
out = wtfusion(x1,x2,2,'db1');

subplot(1,3,1);imshow(x1);title('左模糊');
subplot(1,3,2);imshow(x2);title('右模糊');
subplot(1,3,3);imshow(uint8(out));title('基於小波變換的影象融合');

執行結果