function [ESN,PSN,MUN,SIN]=emlht(w,ES,PS,MU,SI,zm)
%function [ESN,PSN,MUN,SIN]=emlht(w,ES,PS,MU,SI,zm)
% updates LH subband of HMT model once
% by tying within scale for each subband.
%
% Author: H. Choi
% Last update : 12/14/1998
%
% the data structure of 2-D DWT follows the format used
% in the Rice wavelelet matlab toolbox (see README)
%
% internal variables:
% M : no. of mixture densities
% P : size of image (PxP pixels)
% level : no. of levels of HMT
%
% Input :
% w : data PxP matrix (wavelet transform of an image)
% ES : state transition matrix (MxMxPxP)
% PS : state probability matrix (MxPxP)
% MU : mean matrix (MxPxP)
% SI : variance matrix (MxPxP)
% zm : type of density functions
% zm = 1 : zero mean (do not update MU)
% zm = 0 : nonzero mean (update MU)
%
% Output:
% Updated ES PS MU SI in ESN PSN MUN SIN
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);
%UP step
wtmp = repmat(w,[1 1 M]);
wtmp = shiftdim(wtmp,2);
si=1; ei=2^(level-1); 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);
%clear MUtmp SItmp;
for k=level:-1:2
J=2^(k-1);J2=J*J; si = 1; ei = J; sj = J+1; ej = 2*J;
EStmp = reshape(ES(:,:,si:ei,sj:ej),M,M*J2);
if M == 2
%%%%%% For M=2 the following is faster
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
% For general M (not equal to 2) use the following
%
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) = squeeze(sum(BEtmp,1));
sni = 1; eni = J/2; snj = J/2+1; enj = J;
%construct betachild matrix here
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;
%construct BE(:,pai(i),paj(j),dindex) matrix
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 BCtmp Btmp;
%DOWN step
%initialize
AL(:,1,2) = PS(:,1,2);
for k=2:level
J = 2^(k-1); J2=J*J;
si=1; ei=J; sj=J+1; ej=2*J;
sni = 1; eni = J/2; snj = J/2+1; enj = J;
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;
%compute probabilities
for k=2:level
J=2^(k-1); J2=J*J;
si=1; ei=J; sj=J+1; ej=2*J;
sni = 1; eni = J/2; snj = J/2+1; enj = J;
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;
%compute P2
if M == 2
% For M=2 the following may be faster
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
% For general M (not equal to 2) use the following
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(:,1,2)=AL(:,1,2).*BE(:,1,2)./repmat(sum(AL(:,1,2).*BE(:,1,2),1),[M 1 1]);
clear temp BEtmp EStmp Atmp BERtmp;
%PSP=PS; ESP=ES; MUP=MU; SIP=SI;
%M step
PS(:,1,2)=P1(:,1,2);
for k=2:level
J=2^(k-1); J2=J*J;
si=1; ei=J; sj=J+1; ej=2*J;
sni = 1; eni = J/2; snj = J/2+1; enj = J;
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
% do not update MU if zero mean densities
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; %k
ESN=ES; PSN=PS; MUN=MU; SIN=SI;
没有合适的资源?快使用搜索试试~ 我知道了~
资源推荐
资源详情
资源评论
收起资源包目录
MATLAB实现隐马尔可夫预测【数学建模、科学计算算法】.zip (36个子文件)
MATLAB实现隐马尔可夫预测【数学建模、科学计算算法】
whmt1
psnr.m 88B
vec2mat.m 1KB
example2.m 541B
hdenoise.m 4KB
lenamodel.mat 16.5MB
ex512.m 511B
peppers.png 157KB
up.m 976B
lena512 256KB
README 3KB
barbara.png 174KB
dfilters.m 3KB
emhht.m 5KB
posthl.m 3KB
gauss.m 81B
hmtdeno.m 1KB
Idwt.m 802B
emhlt.m 5KB
example3.m 1KB
postlh.m 3KB
Idwt2.m 3KB
emlht.m 5KB
lena256 64KB
example1.m 508B
hmtdeno0.m 1KB
daubcqf.m 2KB
Iidwt2.m 2KB
posthh.m 3KB
testcon.m 187B
ex256.m 511B
makeh1.m 84B
example3.m~ 954B
makenoise.m 1KB
hmtmodel.m 2KB
hmttrain.m 2KB
Iidwt.m 936B
共 36 条
- 1
资源评论
不脱发的程序猿
- 粉丝: 24w+
- 资源: 5804
上传资源 快速赚钱
- 我的内容管理 展开
- 我的资源 快来上传第一个资源
- 我的收益 登录查看自己的收益
- 我的积分 登录查看自己的积分
- 我的C币 登录后查看C币余额
- 我的收藏
- 我的下载
- 下载帮助
安全验证
文档复制为VIP权益,开通VIP直接复制
信息提交成功