% Directorio
myFolder = 'C:\Users\jc_ec\Documents\ProDiVoz\Matlab\digitos_11025Hz';
% Leer digitos de los locutores 01 al 10
[y01,Fs]=leerdig('?01*.wav',myFolder);
[y02,Fs]=leerdig('?02*.wav',myFolder);
[y03,Fs]=leerdig('?03*.wav',myFolder);
[y04,Fs]=leerdig('?04*.wav',myFolder);
[y05,Fs]=leerdig('?05*.wav',myFolder);
[y06,Fs]=leerdig('?06*.wav',myFolder);
[y07,Fs]=leerdig('?07*.wav',myFolder);
[y08,Fs]=leerdig('?08*.wav',myFolder);
[y09,Fs]=leerdig('?09*.wav',myFolder);
[y10,Fs]=leerdig('?10*.wav',myFolder);
N=220; % nro. de muestras por frame
D=N/Fs;  % duracion de un frame
%
% Calculo de los coeficientes cepstral
%
Q=15;    % numero de coeficientes cepstral
yframes01=cell(size(y01));
yframes01w=cell(size(y01));
ycepstrum01=cell(size(y01));
yframes02=cell(size(y02));
yframes02w=cell(size(y02));
ycepstrum02=cell(size(y02));
yframes03=cell(size(y03));
yframes03w=cell(size(y03));
ycepstrum03=cell(size(y03));
yframes04=cell(size(y04));
yframes04w=cell(size(y04));
ycepstrum04=cell(size(y04));
yframes05=cell(size(y05));
yframes05w=cell(size(y05));
ycepstrum05=cell(size(y05));
yframes06=cell(size(y06));
yframes06w=cell(size(y06));
ycepstrum06=cell(size(y06));
yframes07=cell(size(y07));
yframes07w=cell(size(y07));
ycepstrum07=cell(size(y07));
yframes08=cell(size(y08));
yframes08w=cell(size(y08));
ycepstrum08=cell(size(y08));
yframes09=cell(size(y09));
yframes09w=cell(size(y09));
ycepstrum09=cell(size(y09));
yframes10=cell(size(y10));
yframes10w=cell(size(y10));
ycepstrum10=cell(size(y10));
for i=1:30
    y01{i}=y01{i}-mean(y01{i});
    NF01=floor(length(y01{i})/N);
    y01{i}=y01{i}(1:N*NF01);
    yframes01{i}=reshape(y01{i},N,NF01);
    yframes01w{i}=yframes01{i}.*kron(ones(1,NF01),hamming(N));
    ycepstrum01{i}=rceps(yframes01w{i});
    ycepstrum01{i}=ycepstrum01{i}(2:Q,:);
end
for i=1:30
    y02{i}=y02{i}-mean(y02{i});
    NF02=floor(length(y02{i})/N);
    y02{i}=y02{i}(1:N*NF02);
    yframes02{i}=reshape(y02{i},N,NF02);
    yframes02w{i}=yframes02{i}.*kron(ones(1,NF02),hamming(N));
    ycepstrum02{i}=rceps(yframes02w{i});
    ycepstrum02{i}=ycepstrum02{i}(2:Q,:);
end
for i=1:30
    y03{i}=y03{i}-mean(y03{i});
    NF03=floor(length(y03{i})/N);
    y03{i}=y03{i}(1:N*NF03);
    yframes03{i}=reshape(y03{i},N,NF03);
    yframes03w{i}=yframes03{i}.*kron(ones(1,NF03),hamming(N));
    ycepstrum03{i}=rceps(yframes03w{i});
    ycepstrum03{i}=ycepstrum03{i}(2:Q,:);
end
for i=1:30
    y04{i}=y04{i}-mean(y04{i});
    NF04=floor(length(y04{i})/N);
    y04{i}=y04{i}(1:N*NF04);
    yframes04{i}=reshape(y04{i},N,NF04);
    yframes04w{i}=yframes04{i}.*kron(ones(1,NF04),hamming(N));
    ycepstrum04{i}=rceps(yframes04w{i});
    ycepstrum04{i}=ycepstrum04{i}(2:Q,:);
end
for i=1:30
    y05{i}=y05{i}-mean(y05{i});
    NF05=floor(length(y05{i})/N);
    y05{i}=y05{i}(1:N*NF05);
    yframes05{i}=reshape(y05{i},N,NF05);
    yframes05w{i}=yframes05{i}.*kron(ones(1,NF05),hamming(N));
    ycepstrum05{i}=rceps(yframes05w{i});
    ycepstrum05{i}=ycepstrum05{i}(2:Q,:);
end
for i=1:30
    y06{i}=y06{i}-mean(y06{i});
    NF06=floor(length(y06{i})/N);
    y06{i}=y06{i}(1:N*NF06);
    yframes06{i}=reshape(y06{i},N,NF06);
    yframes06w{i}=yframes06{i}.*kron(ones(1,NF06),hamming(N));
    ycepstrum06{i}=rceps(yframes06w{i});
    ycepstrum06{i}=ycepstrum06{i}(2:Q,:);
end
for i=1:30
    y07{i}=y07{i}-mean(y07{i});
    NF07=floor(length(y07{i})/N);
    y07{i}=y07{i}(1:N*NF07);
    yframes07{i}=reshape(y07{i},N,NF07);
    yframes07w{i}=yframes07{i}.*kron(ones(1,NF07),hamming(N));
    ycepstrum07{i}=rceps(yframes07w{i});
    ycepstrum07{i}=ycepstrum07{i}(2:Q,:);
end
for i=1:30
    y08{i}=y08{i}-mean(y08{i});
    NF08=floor(length(y08{i})/N);
    y08{i}=y08{i}(1:N*NF08);
    yframes08{i}=reshape(y08{i},N,NF08);
    yframes08w{i}=yframes08{i}.*kron(ones(1,NF08),hamming(N));
    ycepstrum08{i}=rceps(yframes08w{i});
    ycepstrum08{i}=ycepstrum08{i}(2:Q,:);
end
for i=1:30
    y09{i}=y09{i}-mean(y09{i});
    NF09=floor(length(y09{i})/N);
    y09{i}=y09{i}(1:N*NF09);
    yframes09{i}=reshape(y09{i},N,NF09);
    yframes09w{i}=yframes09{i}.*kron(ones(1,NF09),hamming(N));
    ycepstrum09{i}=rceps(yframes09w{i});
    ycepstrum09{i}=ycepstrum09{i}(2:Q,:);
end
for i=1:30
    y10{i}=y10{i}-mean(y10{i});
    NF10=floor(length(y10{i})/N);
    y10{i}=y10{i}(1:N*NF10);
    yframes10{i}=reshape(y10{i},N,NF10);
    yframes10w{i}=yframes10{i}.*kron(ones(1,NF10),hamming(N));
    ycepstrum10{i}=rceps(yframes10w{i});
    ycepstrum10{i}=ycepstrum10{i}(2:Q,:);
end
%
% ENTRENAMIENTO - Elijo los primeros 30
%
q=1;
M=50;
[pi00, A00, B00, mu00, Sigma00] = init_mhmm(ycepstrum01, q, M,'diag');
[ll_trace01, pi01, A01, mu01, sigma01, B01] = mhmm_em(ycepstrum01(1:30), pi00, A00, mu00, Sigma00, B00,'max_iter', 30,'verbose',0);
[pi10, A10, B10, mu10, Sigma10] = init_mhmm(ycepstrum02, q, M,'diag');
[ll_trace02, pi02, A02, mu02, sigma02, B02] = mhmm_em(ycepstrum02(1:30), pi10, A10, mu10, Sigma10, B10,'max_iter', 30,'verbose',0);
[pi20, A20, B20, mu20, Sigma20] = init_mhmm(ycepstrum03, q, M,'diag');
[ll_trace03, pi03, A03, mu03, sigma03, B03] = mhmm_em(ycepstrum03(1:30), pi20, A20, mu20, Sigma20, B20,'max_iter', 30,'verbose',0);
[pi30, A30, B30, mu30, Sigma30] = init_mhmm(ycepstrum04, q, M,'diag');
[ll_trace04, pi04, A04, mu04, sigma04, B04] = mhmm_em(ycepstrum04(1:30), pi30, A30, mu30, Sigma30, B30,'max_iter', 30,'verbose',0);
[pi40, A40, B40, mu40, Sigma40] = init_mhmm(ycepstrum05, q, M,'diag');
[ll_trace05, pi05, A05, mu05, sigma05, B05] = mhmm_em(ycepstrum05(1:30), pi40, A40, mu40, Sigma40, B40,'max_iter', 30,'verbose',0);
[pi50, A50, B50, mu50, Sigma50] = init_mhmm(ycepstrum06, q, M,'diag');
[ll_trace06, pi06, A06, mu06, sigma06, B06] = mhmm_em(ycepstrum06(1:30), pi50, A50, mu50, Sigma50, B50,'max_iter', 30,'verbose',0);
[pi60, A60, B60, mu60, Sigma60] = init_mhmm(ycepstrum07, q, M,'diag');
[ll_trace07, pi07, A07, mu07, sigma07, B07] = mhmm_em(ycepstrum07(1:30), pi60, A60, mu60, Sigma60, B60,'max_iter', 30,'verbose',0);
[pi70, A70, B70, mu70, Sigma70] = init_mhmm(ycepstrum08, q, M,'diag');
[ll_trace08, pi08, A08, mu08, sigma08, B08] = mhmm_em(ycepstrum08(1:30), pi70, A70, mu70, Sigma70, B70,'max_iter', 30,'verbose',0);
[pi80, A80, B80, mu80, Sigma80] = init_mhmm(ycepstrum09, q, M,'diag');
[ll_trace09, pi09, A09, mu09, sigma09, B09] = mhmm_em(ycepstrum09(1:30), pi80, A80, mu80, Sigma80, B80,'max_iter', 30,'verbose',0);
[pi90, A90, B90, mu90, Sigma90] = init_mhmm(ycepstrum10, q, M,'diag');
[ll_trace10, pi10, A10, mu10, sigma10, B10] = mhmm_em(ycepstrum10(1:30), pi90, A90, mu90, Sigma90, B90,'max_iter', 30,'verbose',0);
%=================================
% Universal Bakground Model (UBM)
%=================================
ycepstrumUBM=[ycepstrum01(:);ycepstrum02(:);ycepstrum03(:);ycepstrum04(:);ycepstrum05(:);...
    ycepstrum06(:);ycepstrum07(:);ycepstrum08(:);ycepstrum09(:);ycepstrum10(:)];
[pi0UBM, A0UBM, B0UBM, mu0UBM, Sigma0UBM] = init_mhmm(ycepstrumUBM, q, M,'diag');
[ll_traceUBM, piUBM, AUBM, muUBM, sigmaUBM, BUBM] = mhmm_em(ycepstrumUBM(1:300),...
    pi0UBM, A0UBM, mu0UBM, Sigma0UBM, B0UBM,'max_iter', 30,'verbose',0);
%========================
% RECONOCIMIENTO
%========================
yincognita=ycepstrum08;
[loglik01, errors01] = mhmm_logprob(yincognita, pi01, A01, mu01, sigma01, B01);
[loglik02, errors02] = mhmm_logprob(yincognita, pi02, A02, mu02, sigma02, B02);
[loglik03, errors03] = mhmm_logprob(yincognita, pi03, A03, mu03, sigma03, B03);
[loglik04, errors04] = mhmm_logprob(yincognita, pi04, A04, mu04, sigma04, B04);
[loglik05, errors05] = mhmm_logprob(yincognita, pi05, A05, mu05, sigma05, B05);
[loglik06, errors06] = mhmm_logprob(yincognita, pi06, A06, mu06, sigma06, B06);
[loglik07, errors07] = mhmm_logprob(yincognita, pi07, A07, mu07, sigma07, B07);
[loglik08, errors08] = mhmm_logprob(yincognita, pi08, A08, mu08, sigma08, B08);
[loglik09, errors09] = mhmm_logprob(yincognita, pi09, A09, mu09, sigma09, B09);
[loglik10, errors10] = mhmm_logprob(yincognita, pi10, A10, mu10, sigma10, B10);
%
[loglikUBM, errorsUBM] = mhmm_logprob(yincognita, piUBM, AUBM, muUBM, sigmaUBM, BUBM);
loglikelihoodratio=[loglik01,loglik02,loglik03,loglik04,loglik05,loglik06,loglik07,...
    loglik08,loglik09,loglik10]/loglikUBM;
[maxi,locutor]=max(loglikelihoodratio);