MATLAB 2维小波变换程序源代码

0

%  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');

3