语音信号隐马尔可夫模型仿真 第10页

附录二

Baum-Welch的完整的实验程序

function [lik,likv]=hmm_cl(X,T,K,Mu,Cov,P,Pi);

p=length(X(1,:));

N=length(X(:,1));

tiny=exp(-700);

if (rem(N,T)~=0)

  disp('Error: Data matrix length must be multiple of sequence length T');

  return;

end;

N=N/T;

alpha=zeros(T,K);

B=zeros(T,K);      % P( output | s_i)

k1=(2*pi)^(-p/2);

Scale=zeros(T,1);

likv=zeros(1,N);

for n=1:N

    B=zeros(T,K);

  iCov=inv(Cov);     

  k2=k1/sqrt(det(Cov));

  for i=1:T

    for l=1:K

      d= Mu(l,:)-X((n-1)*T+i,:);

      B(i,l)=k2*exp(-0.5*d*iCov*d');

    end;

  end;

  scale=zeros(T,1);

  alpha(1,:)=Pi(:)'.*B(1,:);

  scale(1)=sum(alpha(1,:));

  alpha(1,:)=alpha(1,:)/(scale(1)+tiny);

  for i=2:T

    alpha(i,:)=(alpha(i-1,:)*P).*B(i,:);

    scale(i)=sum(alpha(i,:));

    alpha(i,:)=alpha(i,:)/(scale(i)+tiny);

  end;

  likv(n)=sum(log(scale+(scale==0)*tiny));

  Scale=Scale+log(scale+(scale==0)*tiny);

end;

lik=sum(Scale);

**************************************************************

% function Z=rdiv(X,Y)

% row division: Z = X / Y row-wise

% Y must have one column

 function Z=rdiv(X,Y)

 if(length(X(:,1)) ~= length(Y(:,1)) | length(Y(1,:)) ~=1)

  disp('Error in RDIV');

  return;

end

 Z=zeros(size(X));

 for i=1:length(X(1,:))

  Z(:,i)=X(:,i)./Y;

end

**************************************************************

% row product

% function Z=rprod(X,Y)

 function Z=rprod(X,Y)

 [n m]=size(X);

 if(length(Y(:,1)) ~=n  | length(Y(1,:)) ~=1)

  disp('Error in RPROD');

  return;

end

 Z=X.*(Y*ones(1,m));

**************************************************************

% row sum

% function Z=rsum(X)

function Z=rsum(X)

Z=zeros(size(X(:,1)));

for i=1:length(X(1,:))

  Z=Z+X(:,i);

end

**************************************************************

function [Mu,Cov,P,Pi,LL]=hmm(X,T,K,cyc,tol)

p=length(X(1,:));

N=length(X(:,1));

if nargin<5   tol=0.0001; end;

if nargin<4   cyc=100; end;

if nargin<3   K=2; end;

if nargin<2   T=N; end;

 if (rem(N,T)~=0)

  disp('Error: Data matrix length must be multiple of sequence length T');

  return;

end;

N=N/T;

Cov=diag(diag(cov(X)));

Mu=randn(K,p)*sqrtm(Cov)+ones(K,1)*mean(X);

Pi=rand(1,K);

Pi=Pi/sum(Pi);

P=rand(K);

P=rdiv(P,rsum(P));

LL=[];

lik=0;

alpha=zeros(T,K);

beta=zeros(T,K);

gamma=zeros(T,K);

B=zeros(T,K);

k1=(2*pi)^(-p/2);

for cycle=1:cyc

  %%%% FORWARD-BACKWARD

  Gamma=[];

上一页  [1] [2] [3] [4] [5] [6] [7] [8] [9] [10]  ... 下一页  >>