2維小波變換經典程序
% FWT_DB.M;
% 此示意程序用DWT實現二維小波變換
% 編程時間2004-4-10，編程人沙威
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;clc;
T=256; % 圖像維數
SUB_T=T/2; % 子圖維數
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1.調原始圖像矩陣
load wbarb; % 下載圖像
f=X; % 原始圖像
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2.進行二維小波分解
l=wfilters('db10','l'); % db10（消失矩為10)低通分解濾波器沖擊響應（長度為20）
L=T-length(l);
l_zeros=[l,zeros(1,L)]; % 矩陣行數與輸入圖像一致，為2的整數冪
h=wfilters('db10','h'); % db10（消失矩為10)高通分解濾波器沖擊響應（長度為20）
h_zeros=[h,zeros(1,L)]; % 矩陣行數與輸入圖像一致，為2的整數冪
for i=1:T; % 列變換
row(1:SUB_T,i)=dyaddown( ifft( fft(l_zeros).*fft(f(:,i)') ) ).'; % 圓周卷積<->FFT
row(SUB_T+1:T,i)=dyaddown( ifft( fft(h_zeros).*fft(f(:,i)') ) ).'; % 圓周卷積<->FFT
end;
for j=1:T; % 行變換
line(j,1:SUB_T)=dyaddown( ifft( fft(l_zeros).*fft(row(j,:)) ) ); % 圓周卷積<->FFT
line(j,SUB_T+1:T)=dyaddown( ifft( fft(h_zeros).*fft(row(j,:)) ) ); % 圓周卷積<->FFT
end;
decompose_pic=line; % 分解矩陣
% 圖像分為四塊
lt_pic=decompose_pic(1:SUB_T,1:SUB_T); % 在矩陣左上方為低頻分量--fi(x)*fi(y)
rt_pic=decompose_pic(1:SUB_T,SUB_T+1:T); % 矩陣右上為--fi(x)*psi(y)
lb_pic=decompose_pic(SUB_T+1:T,1:SUB_T); % 矩陣左下為--psi(x)*fi(y)
rb_pic=decompose_pic(SUB_T+1:T,SUB_T+1:T); % 右下方為高頻分量--psi(x)*psi(y)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3.分解結果顯示
figure(1);
colormap(map);
subplot(2,1,1);
image(f); % 原始圖像
title('original pic');
subplot(2,1,2);
image(abs(decompose_pic)); % 分解后圖像
title('decomposed pic');
figure(2);
colormap(map);
subplot(2,2,1);
image(abs(lt_pic)); % 左上方為低頻分量--fi(x)*fi(y)
title('\Phi(x)*\Phi(y)');
subplot(2,2,2);
image(abs(rt_pic)); % 矩陣右上為--fi(x)*psi(y)
title('\Phi(x)*\Psi(y)');
subplot(2,2,3);
image(abs(lb_pic)); % 矩陣左下為--psi(x)*fi(y)
title('\Psi(x)*\Phi(y)');
subplot(2,2,4);
image(abs(rb_pic)); % 右下方為高頻分量--psi(x)*psi(y)
title('\Psi(x)*\Psi(y)');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 5.重構源圖像及結果顯示
% construct_pic=decompose_matrix'*decompose_pic*decompose_matrix;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
l_re=l_zeros(end:-1:1); % 重構低通濾波
l_r=circshift(l_re',1)'; % 位置調整
h_re=h_zeros(end:-1:1); % 重構高通濾波
h_r=circshift(h_re',1)'; % 位置調整

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
top_pic=[lt_pic,rt_pic]; % 圖像上半部分
t=0;
for i=1:T; % 行插值低頻

if (mod(i,2)==0)
topll(i,:)=top_pic(t,:); % 偶數行保持
else
t=t+1;
topll(i,:)=zeros(1,T); % 奇數行為零
end
end;
for i=1:T; % 列變換
topcl_re(:,i)=ifft( fft(l_r).*fft(topll(:,i)') )'; % 圓周卷積<->FFT
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
bottom_pic=[lb_pic,rb_pic]; % 圖像下半部分
t=0;
for i=1:T; % 行插值高頻
if (mod(i,2)==0)
bottomlh(i,:)=bottom_pic(t,:); % 偶數行保持
else
bottomlh(i,:)=zeros(1,T); % 奇數行為零
t=t+1;
end
end;
for i=1:T; % 列變換
bottomch_re(:,i)=ifft( fft(h_r).*fft(bottomlh(:,i)') )'; % 圓周卷積<->FFT
end;

construct1=bottomch_re+topcl_re; % 列變換重構完畢

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
left_pic=construct1(:,1:SUB_T); % 圖像左半部分
t=0;
for i=1:T; % 列插值低頻

if (mod(i,2)==0)
leftll(:,i)=left_pic(:,t); % 偶數列保持
else
t=t+1;
leftll(:,i)=zeros(T,1); % 奇數列為零
end
end;
for i=1:T; % 行變換
leftcl_re(i,:)=ifft( fft(l_r).*fft(leftll(i,:)) ); % 圓周卷積<->FFT
end;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
right_pic=construct1(:,SUB_T+1:T); % 圖像右半部分

t=0;
for i=1:T; % 列插值高頻
if (mod(i,2)==0)
rightlh(:,i)=right_pic(:,t); % 偶數列保持
else
rightlh(:,i)=zeros(T,1); % 奇數列為零
t=t+1;
end
end;
for i=1:T; % 行變換
rightch_re(i,:)=ifft( fft(h_r).*fft(rightlh(i,:)) ); % 圓周卷積<->FFT
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
construct_pic=rightch_re+leftcl_re; % 重建全部圖像

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 結果顯示
figure(3);
colormap(map);
subplot(2,1,1);
image(f); % 源圖像顯示
title('original pic');

subplot(2,1,2);
image(abs(construct_pic)); % 重構源圖像顯示
title('reconstructed pic');
error=abs(construct_pic-f); % 重構圖形與原始圖像誤值
figure(4);
mesh(error); % 誤差三維圖像
title('absolute error display');