% wavdate
% 用HMM模型進行孤立詞的語音識別

function featurecluster
load wavdate         %調入wavda
for T00=2:2:16
    T00
for fg=1:37
    fg
    alpha=(fg-1)/36;
% train
M=36;       %語音信號0-9，A-Z
N=10;       %說話人1-10
for i = 1:M   
    i
	sample1= [];
	for k = 1:N  
        a=wavdate((k-1)*72+2*i-1).endPoint(1);
        b=wavdate((k-1)*72+2*i-1).endPoint(2);
        x = filter([1 -0.9375], 1, double(wavdate((k-1)*72+2*i-1).y(a:b)));
		%sample1(k).data = melcepst(x,16000,'M',12,29,256,80); %提取mel頻率倒譜系數12階
        %sample1(k).data = lpcc(x);                            %提取lpcc12階及差分
        %sample1(k).data=gbfenmfccfeature(x,alpha);
        %sample1(k).data=gbfcosmelcepst(alpha,x,16000,'M',12,29,256,80,0,0.5);%效果好
        %sample1(k).data=gbfenlpccfeature(x,alpha);
        sample1(k).data=gbfcosmelcepstzhouqi(T00,alpha,x,16000,'M',12,29,256,80,0,0.5);%效果好
    end
end
D0(fg)=Dratio(sample1)
end
figure(T00+10);plot([0:1/36:1],D0,'k');hold on;
clear D0
end


%%%%%%%%%%%%%%5
function D=Dratio(x)
% input--- x is a struct array variable, cluding m persons feature matrix f(n*k), x(i).f
%thereinto, n indicate speech number, k indicate feature vector dimension
[n,m]=size(x);% n is equal to "1", and m indicate the speakers number
%calculate W
n
m
for i=1:m
    [n,k]=size(x(i).data);
    u(i,:)=mean(x(i).data);
    for j=1:n
        xu=x(i).data(j,:)-u(i,:);
        y(i).z(j).a=xu'*xu;
    end
    summ=0;
    for j=1:n
        summ=summ+y(i).z(j).a;
    end
    y(i).b=summ/n;
end
sumb=0;
for i=1:m
    sumb=sumb+y(i).b;
end
W=sumb/m;
%calculate B
Bv=cov(u);
uu=mean(u);
for i=1:m
    ud=u(i,:)-uu;
    y(i).c=ud'*ud;
end
sumc=0;
for i=1:m
    sumc=sumc+y(i).c;
end
B=sumc/m;
%calculate D
D=trace(W^(-1)*Bv);