hmm.m
上传用户:ay_070428
上传日期:2014-12-04
资源大小:11427k
文件大小:3k
源码类别:

语音合成与识别

开发平台:

Matlab

  1. % function [Mu,Cov,P,Pi,LL]=hmm(X,T,K,cyc,tol);
  2. % Gaussian Observation Hidden Markov Model
  3. %
  4. % X - N x p data matrix
  5. % T - length of each sequence (N must evenly divide by T, default T=N)
  6. % K - number of states (default 2)
  7. % cyc - maximum number of cycles of Baum-Welch (default 100)
  8. % tol - termination tolerance (prop change in likelihood) (default 0.0001)
  9. %
  10. % Mu - mean vectors
  11. % Cov - output covariance matrix (full, tied across states)
  12. % P - state transition matrix
  13. % Pi - priors
  14. % LL - log likelihood curve
  15. %
  16. % Iterates until a proportional change < tol in the log likelihood 
  17. % or cyc steps of Baum-Welch
  18. function [Mu,Cov,P,Pi,LL]=hmm(X,T,K,cyc,tol)
  19. p=length(X(1,:));
  20. N=length(X(:,1));
  21. if nargin<5   tol=0.0001; end;
  22. if nargin<4   cyc=100; end;
  23. if nargin<3   K=2; end;
  24. if nargin<2   T=N; end;
  25. if (rem(N,T)~=0)
  26.   disp('Error: Data matrix length must be multiple of sequence length T');
  27.   return;
  28. end;
  29. N=N/T;
  30. Cov=diag(diag(cov(X)));
  31. Mu=randn(K,p)*sqrtm(Cov)+ones(K,1)*mean(X);
  32. Pi=rand(1,K);
  33. Pi=Pi/sum(Pi);
  34. P=rand(K);
  35. P=rdiv(P,rsum(P));
  36. LL=[];
  37. lik=0;
  38. alpha=zeros(T,K);
  39. beta=zeros(T,K);
  40. gamma=zeros(T,K);
  41. B=zeros(T,K);
  42. k1=(2*pi)^(-p/2);
  43. for cycle=1:cyc
  44.   
  45.   %%%% FORWARD-BACKWARD 
  46.   
  47.   Gamma=[];
  48.   Gammasum=zeros(1,K);
  49.   Scale=zeros(T,1);
  50.   Xi=zeros(T-1,K*K);
  51.   
  52.   for n=1:N
  53.     
  54.     iCov=inv(Cov);
  55.     k2=k1/sqrt(det(Cov));
  56.     for i=1:T
  57.       for l=1:K
  58. d=Mu(l,:)-X((n-1)*T+i,:);
  59. B(i,l)=k2*exp(-0.5*d*iCov*d');
  60.       end;
  61.     end;
  62.     
  63.     scale=zeros(T,1);
  64.     alpha(1,:)=Pi.*B(1,:);
  65.     scale(1)=sum(alpha(1,:));
  66.     alpha(1,:)=alpha(1,:)/scale(1);
  67.     for i=2:T
  68.       alpha(i,:)=(alpha(i-1,:)*P).*B(i,:);
  69.       scale(i)=sum(alpha(i,:));
  70.       alpha(i,:)=alpha(i,:)/scale(i);
  71.     end;
  72.     
  73.     beta(T,:)=ones(1,K)/scale(T);
  74.     for i=T-1:-1:1
  75.       beta(i,:)=(beta(i+1,:).*B(i+1,:))*(P')/scale(i); 
  76.     end;
  77.     
  78.     gamma=(alpha.*beta); 
  79.     gamma=rdiv(gamma,rsum(gamma));
  80.     gammasum=sum(gamma);
  81.     
  82.     xi=zeros(T-1,K*K);
  83.     for i=1:T-1
  84.       t=P.*( alpha(i,:)' * (beta(i+1,:).*B(i+1,:)));
  85.       xi(i,:)=t(:)'/sum(t(:));
  86.     end;
  87.     
  88.     Scale=Scale+log(scale);
  89.     Gamma=[Gamma; gamma];
  90.     Gammasum=Gammasum+gammasum;
  91.     Xi=Xi+xi;
  92.   end;
  93.   
  94.   %%%% M STEP 
  95.   
  96.   % outputs
  97.   Mu=zeros(K,p);
  98.   Mu=Gamma'*X;
  99.   Mu=rdiv(Mu,Gammasum');
  100.   
  101.   % transition matrix 
  102.   sxi=rsum(Xi')';
  103.   sxi=reshape(sxi,K,K);
  104.   P=rdiv(sxi,rsum(sxi));
  105.   
  106.   % priors
  107.   Pi=zeros(1,K);
  108.   for i=1:N
  109.     Pi=Pi+Gamma((i-1)*T+1,:);
  110.   end
  111.   Pi=Pi/N;
  112.   
  113.   % covariance
  114.   Cov=zeros(p,p);
  115.   for l=1:K
  116.     d=(X-ones(T*N,1)*Mu(l,:));
  117.     Cov=Cov+rprod(d,Gamma(:,l))'*d;
  118.   end;
  119.   Cov=Cov/(sum(Gammasum));
  120.   
  121.   oldlik=lik;
  122.   lik=sum(Scale);
  123.   LL=[LL lik];
  124.   fprintf('cycle %i log likelihood = %f ',cycle,lik);  
  125.   
  126.   if (cycle<=2)
  127.     likbase=lik;
  128.   elseif (lik<oldlik) 
  129.     fprintf('violation');
  130.   elseif ((lik-likbase)<(1 + tol)*(oldlik-likbase)|~finite(lik)) 
  131.     fprintf('n');
  132.     break;
  133.   end;
  134.   fprintf('n');
  135. end