wavelet_lab2.m

主要是学习小波的基础

clear
clc
[X,map]=imread('LENA.bmp');   % 载入图像LENA.bmp，需要更改路径到指定的文件
X=double(X);
% nbcol=255;
%cod_X=wcodemat(X,255);
%figure
colormap(pink(255));
% subplot(2,2,1);
image(X);
%image(cod_X);       % 显示原始图像

% 生成小波滤波器组
[Lo_D,Hi_D,Lo_R,Hi_R]=wfilters('sym3');

sx=size(X);         % 计算原始图像信号的系数矩阵大小
%coefs=X;
c=X; 
len=length(X);      % 计算信号的长度      

% 第1级的边缘延拓
ext_c=[];
ext_c=[ext_c c];
lf=16;              % 设置延拓的扩展长度
type=0;             % 选定延拓方式
switch type
  case 0            % 连续延拓
       for t=1:lf
           ext_c=[ext_c(:,1) ext_c ext_c(:,len)];       % 对信号的列延拓
       end
       ext_c=ext_c';
       for k=1:lf
           ext_c=[ext_c(:,1) ext_c ext_c(:,len)];       % 对信号的行延拓
       end
       ext_c=ext_c';

  case 1            % 周期延拓
       for t=1:lf
           ext_c=[ext_c(:,len-t+1) ext_c ext_c(:,t)];   % 对信号的列延拓
       end
       ext_c=ext_c';
       for t=1:lf
           ext_c=[ext_c(:,len-t+1) ext_c ext_c(:,t)];   % 对信号的行延拓
       end
       ext_c=ext_c';  

  case 2            % 对称延拓
       for t=1:lf
           ext_c=[ext_c(:,t) ext_c ext_c(:,len-t+1)];   % 对信号的列延拓
       end
       ext_c=ext_c';
       for t=1:lf
           ext_c=[ext_c(:,t) ext_c ext_c(:,len-t+1)];   % 对信号的行延拓
       end
       ext_c=ext_c';  

  case 3            % 补零延拓
       temp=zeros(len,lf);
       ext_c=[temp ext_c temp];                         % 对信号的列延拓
       ext_c=ext_c';
       temp=zeros(len+2*lf,lf);
       ext_c=[temp ext_c temp];                         % 对信号的行延拓
       ext_c=ext_c';               
end

coefs=ext_c;                                 % 第1级延拓后的信号系数矩阵
    
% 分解
for i=1:3
    len=(sx(1)+lf*2)/(2^(i-1));              % 计算当前级数上，进行分解的信号长度
    dwt_c=coefs(1:len,1:len);                % 在coefs中选定进行变换的信号
        
    if i>1 
    % 第2，3级的边缘延拓
    ext_c=[];
    ext_c=[ext_c dwt_c];    
    %type=0;
    switch type
      case 0            % 连续延拓
           for t=1:lf
               ext_c=[dwt_c(:,1) ext_c dwt_c(:,len)];
           end
           ext_c=ext_c';
           for k=1:lf
               ext_c=[ext_c(:,1) ext_c ext_c(:,len)];
           end
           ext_c=ext_c';

      case 1            % 周期延拓
           for t=1:lf
               ext_c=[ext_c(:,len-t+1) ext_c ext_c(:,t)];
           end
           ext_c=ext_c';
           for t=1:lf
               ext_c=[ext_c(:,len-t+1) ext_c ext_c(:,t)];
           end
           ext_c=ext_c';  

      case 2            % 对称延拓
           for t=1:lf
               ext_c=[ext_c(:,t) ext_c ext_c(:,len-t+1)];
           end
           ext_c=ext_c';
           for t=1:lf
               ext_c=[ext_c(:,t) ext_c ext_c(:,len-t+1)];
           end
           ext_c=ext_c';  
    
      case 3            % 补零延拓
           temp=zeros(len,lf);
           ext_c=[temp ext_c temp];
           ext_c=ext_c';
           temp=zeros(len+2*lf,lf);
           ext_c=[temp ext_c temp];
           ext_c=ext_c';               
    end
    
    dwt_c=ext_c;
end
    
    sizeKEPT=[len len];                             % 设置分解信号的保留长度

    tem1=wkeep(conv2(dwt_c,Lo_D(:)'),sizeKEPT);     % 对信号的每一行与分解低通滤波器作卷积，并取中间sizeKEPT个值
    tem2=wkeep(conv2(dwt_c,Hi_D(:)'),sizeKEPT);     % 对信号的每一行与分解高通滤波器作卷积，并取中间sizeKEPT个值
    dwt_c=[tem1(:,1:2:len) tem2(:,1:2:len)];        % 对两次卷积后的信号分别作1/2采样，合并得到第一步分解的信号
    
    tem=conv2(dwt_c',Lo_D(:)');                     % 对信号的每一列与分解低通滤波器作卷积
    tem1=wkeep(tem',sizeKEPT);                      % 取中间sizeKEPT个值
    tem=conv2(dwt_c',Hi_D(:)');                     % 对信号的每一行与分解低通滤波器作卷积
    tem2=wkeep(tem',sizeKEPT);                      % 取中间sizeKEPT个值
    dwt_c=[tem1(1:2:len,:);tem2(1:2:len,:)];        % 对两次卷积后的信号分别作1/2采样，合并得到第二步分解的信号

    coefs(1:len,1:len)=dwt_c;                       % 分解结果替换到coefs中；当i=3时，得到的coefs即为分解完成的图像信号系数矩阵
end    

% 显示分解后的各级图像分量
cod_coefs=wcodemat(coefs,255);
figure
colormap(pink(255));
image(cod_coefs);
%subplot(2,2,4);
%image(cod_coefs);
%axis('square')

% 重构
icoefs=coefs;                                 
for j=1:3
    t=(sx(1)+lf*2)/(2^(3-j));                      % 计算当前级数上，进行重构的信号长度
    sizeKEPT=[t t];                                % 设置重构信号的保留长度

    idwt_c=zeros(t,t);                             
    idwt_c(1:2:t,:)=icoefs(1:t/2,1:t);             % 对信号前一半的列作上采样，即隔列插0
    tem=conv2(idwt_c',Lo_R(:)');                   % 与重构低通滤波器做卷积
    tem1=wkeep(tem',sizeKEPT,'c',[0 0]);           % 取中间sizeKEPT个值
 
    idwt_c=zeros(t,t);
    idwt_c(1:2:t,:)=icoefs((1+t/2):t,1:t);         % 对信号后一半的列作上采样，即隔列插0
    tem=conv2(idwt_c',Hi_R(:)');                   % 与重构高通滤波器做卷积 
    tem2=wkeep(tem',sizeKEPT,'c',[0 0]);           % 取中间sizeKEPT个值
    idwt_c=tem1+tem2;                              % 合并两次结果，得到第一步重构的信号
    icoefs(1:t,1:t)=idwt_c;
    
    idwt_c=zeros(t,t);
    idwt_c(:,1:2:t)=icoefs(1:t,1:t/2);                                   % 对信号前一半的行作上采样，即隔行插0
    tem1=wkeep(conv2(idwt_c,Lo_R(:)'),sizeKEPT,'c',[0 0]);               % 与重构低通滤波器做卷积，并取中间sizeKEPT个值 
    idwt_c=zeros(t,t);
    idwt_c(:,1:2:t)=icoefs(1:t,(1+t/2):t);                               % 对信号后一半的行作上采样，即隔行插0
    tem2=wkeep(conv2(idwt_c,Hi_R(:)'),sizeKEPT,'c',[0 0]);               % 与重构高通滤波器做卷积，并取中间sizeKEPT个值 
    idwt_c=tem1+tem2;                                                    % 合并两次结果，得到第二步重构的信号

    icoefs(1:t,1:t)=idwt_c;                        % 重构结果替换到icoefs中；当i=3时，得到的icoefs即为重构完成的图像信号系数矩阵
end
icoefs=icoefs(17:272,17:272);                      % 去掉延拓时扩展的长度，得到与原始图像大小相同的重构图像

% 显示重构图像
cod_icoefs=wcodemat(icoefs,255);
figure
colormap(pink(255));
%subplot(2,2,3);
image(cod_icoefs);
%axis('square')


% 计算峰值信噪比PSNR，结果将直接显示在Matlab的Command Window窗口中
I=X-icoefs;
I=I.*I;
[m,n]=size(I);
I=cumsum(I);
sum=cumsum(I(m,:));
mse=sum(n)/(m*n);
psnr=10*log10(255*255/mse)