function [mdl] = train(X, Y)
  %% X1 is a D*N matrix
  classes = unique(Y);
  for i = 1:length(classes),
    clsidx = find(Y==classes(i));
    mdl(i).class = classes(i);
    [mdl(i).mu, mdl(i).Sigma]= getgauss(X(:,clsidx));
  end

function [mu,sigma]=getgauss(X1)
  mu = mean(X1');
  sigma = cov(X1');