更新hh子带的hmt.m

用MATLAB编写的HMT参数训练的HH子带更新子函数

function [ESN,PSN,MUN,SIN]=EMHHT(w,ES,PS,MU,SI,zm)
M=size(ES,1);
P=size(w,1);
level=log2(P);
BE=zeros(M,P,P);
BEP=zeros(M,P,P);
BER=zeros(M,P,P);

AL=zeros(M,P,P);
P1=zeros(M,P,P);
P2=zeros(M,M,P,P);
wtmp=repmat(w,[1 1 M]);
wtmp=shiftdim(wtmp,2);
si=2^(level-1)+1;ei=P;sj=2^(level-1)+1;ej=P;
gtmp=gauss(wtmp,MU,SI);
scale=repmat(mean(gtmp,1),[M 1 1]);
BE(:,si:ei,sj:ej)=gtmp(:,si:ei,sj:ej)./scale(:,si:ei,sj:ej);

for k=level:-1:2
    J=2^(k-1);J2=J*J;si=J+1;ei=2*J;sj=J+1;ej=2*J;
    EStmp=reshape(ES(:,:,si:ei,sj:ej),M,M*J2);
    if M==2
        BEtmp=zeros(M,M*J2);
        BEtmp(:,1:M:(M*J2))=reshape(BE(:,si:ei,sj:ej),M,J2);
        BEtmp(:,2:M:(M*J2))=BEtmp(:,1:M:(M*J2));
    else
        BEtmp=zeros(M,M*4^(k-1));
        for m=1:M
            BEtmp(:,m:M:(M*4^(k-1)))=reshape(BE(:,si:ei,sj:ej,:),M,4^(k-1));
        end;
    end
     BEtmp=reshape(EStmp.*BEtmp,[M M J J]);
     BEP(:,si:ei,sj:ej)=aqueeze(sum(BEtmp, 1));
     
     sni=J/2+1;eni=si-1;snj=J/2+1;enj=sj-1;
     BCtmp=BEP(:,si:2:ei,sj:2:ej);
     BCtmp=BCtmp.*BEP(:,si+1:2:ei,sj:2:ej);
     BCtmp=BCtmp.*BEP(:,si:2:ei,sj+1:2:ej);
     BCtmp=BCtmp.*BEP(:,si+1:2:ei,sj+1:2:ej);
     scaletmp=repmat(mean(BCtmp,1),[M 1 1]);
     scale(:,sni:eni,snj:enj)=scale(:,sni:eni,snj:enj).*scaletmp;
     BE(:,sni:eni,snj:enj)=gtmp(:,sni:eni,snj:enj)./scale(:,sni:eni,snj:enj).*BCtmp;
     
     Btmp=zeros(M,J,J);
     Btmp(:,1:2:J,1:2:J)=BE(:,sni:eni,snj:enj);
     Btmp(:,2:2:J,1:2:J)=BE(:,sni:eni,snj:enj);
     Btmp(:,1:2:J,2:2:J)=BE(:,sni:eni,snj:enj);
     Btmp(:,2:2:J,2:2:J)=BE(:,sni:eni,snj:enj);
     BER(:,si:ei,sj:ej)=Btmp.BEP(:,si:ei,sj:ej);
end

clear EStmp BEtmp Btmp
AL(:,2,2)=PS(:,2,2);
for k=2:level
    J=2^(k-1);J2=J*J;
    si=J+1;ei=2*J;sj=J+1;ej=2*J;
    sni=J/2+1;eni=si-1;snj=J/2+1;enj=sj-1;
    Atmp=zeros(M,J,J);
    Atmp(:,1:2:J,1:2:J)=AL(:,sni:eni,snj:enj);
    Atmp(:,2:2:J,1:2:J)=AL(:,sni:eni,snj:enj);
    Atmp(:,1:2:J,2:2:J)=AL(:,sni:eni,snj:enj);
    Atmp(:,2:2:J,2:2:J)=AL(:,sni:eni,snj:enj);
    
    Atmp=repmat(reshape(Atmp.*BER(:,si:ei,sj:ej),1,M*J2),[M 1]);
    EStmp=reshape(ES(:,:,si:ei,sj:ej),M,M*J2);
    ALtmp=reshape(EStmp.*Atmp,[M M J J]);
    AL(:,si:ei,sj:ej)=squeeze(sum(ALtmp,2));
end;

clear Atmp EStmp ALtmp;
for k=2:level
    J=2^(k-1);J2=J*J;
    si=J+1;ei=2*J;sj=J+1;ej=2*J;
    sni=J/2+1;eni=si-1;snj=J/2+1;enj=sj-1;
    
    temp=repmat(sum(AL(:,si:ei,sj:ej).*BE(:,si:ei,sj:ej),1),[M 1]);
    P1(:,si:ei,sj:ej)=AL(:,si:ei,sj:ej).*BE(:,si:ei,sj:ej)./temp;
   if m==2
    BEtmp=zeros(M,M*J2);
    BEtmp(:,1:M:(M*J2))=reshape(BE(:,si:ei,sj:ej),M,J2);
    BEtmp(:,2:M:(M*J2))=BEtmp(:,1:M:(M*J2));
   else
       BEtmp=zeros(M,M*J2);
       for m=1:M
           BEtmp(:,m:M:(M*J2))=reshape(BE(:,si:ei,sj:ej,:),M,J2);
       end
   end
   
   BEtmp=reshape(BEtmp,[M M J J]);
   
   EStmp=ES(:,:,si:ei,sj:ej);
   Atmp=zeros(M,J,J);
   Atmp(:,1:2:J,1:2:J)=AL(:,sni:eni,snj:enj);
   Atmp(:,2:2:J,1:2:J)=AL(:,sni:eni,snj:enj);
   Atmp(:,1:2:J,2:2:J)=AL(:,sni:eni,snj:enj);
   Atmp(:,2:2:J,2:2:J)=AL(:,sni:eni,snj:enj);
   Atmp=repmat(reshape(Atmp,1,M*J2),[M 1]);
   Atmp=reshape(Atmp,[M M J J]);
   
   BERtmp=repmat(reshape(BER(:,si:ei,sj:ej),1,M*J2),[M 1]);
   BERtmp=reshape(BERtmp,[M M J J]);
   temp=repmat(reshape(temp,1,M*J2),[M 1]);
   temp=reshape(temp,[M M J J]);
   P2(:,:,si:ei,sj:ej)=BEtmp.*EStmp.*Atmp.*BERtmp./temp;
end
P1(:,2,2)=AL(:,2,2).*BE(:,2,2)./REPMAT(SUM(AL(:,2,2).BE(:,2,2),1),[M 1 1]);
 clear temp BEtmp EStmp Atmp BERtmp;
 PS(:,2,2)=P1(:,2,2);
 
 for k=2:level
     J=2^(k-1);J2=J*J;
     si=J+1;ei=2*J;sj=J+1;ej=2*J;
     sni=J/2+1;eni=si-1;snj=J/2+1;enj=sj-1;
     
     pstmp=sum(sum(P1(:,si:ei,sj:ej),3),2)/J2;
     pstmp=pstmp.*(pstmp>1e-4)+1e-4*(pstmp<=1e-4);
     PS(:,si:ei,sj:ej)=repmat(pstmp,[1 J J]);
     
     if zm==0
         mutmp=sum(sum(wtmp(:,si:ei,sj:ej).*P1(:,si:ei,sj:ej),3),2)/J2;
         MU(:,si:ei,sj:ej)=repmat(mutmp,[1 J J])./PS(:,si:ei,sj:ej);
     end
     
     sitmp=sum(sum((wtmp(:,si:ei,sj:ej)-MU(:,si:ei,sj:ej)).^2.*P1(:,si:ei,sj:ej),3),2)/J2;
     SI(:,si:ei,sj:ej)=repmat(sitmp,[1 J J])./PS(:,si:ei,sj:ej);
     
     estmp=sum(sum(P2(:,:,si:ei,sj:ej),4),3)/J2;
     ptmp=[PS(:,sni,snj)':PS(:,sni,snj)'];
     ES(:,:,si:ei,sj:ej)=repmat(estmp,[1 1 J J])./repmat(ptmp,[1 1 J J]);
 end
 ESN=ES;PSN=PS;MUN=MU;SIN=SI;