Automatyczne rozpoznawanie mowy
Autorzy: Ryszard Makowski, Paweł Świętojański, Robert Wielgat
Pliki: rozdzial_14.zip
Przykład 14.1
Przykład (plik przyklad_14_01.m)przedstawia sposób dokonywania podziału sygnału mowy na ramki. Została tu wykorzystana funkcja "ramkowanie", której listing przedstawiono poniżej oraz w programie 14.1 na str. 509. Pominięcie argumentu wyjściowego w wywołaniu funkcji "ramkowanie" powoduje wyświetlenie 10 pierwszych ramek sygnału. Funkcję z argumentem wyjściowym zastosowano w procedurze "melcepstrum.m" (przedstawioną w przykładzie 14.3) służącej do obliczania współczynników MFCC.% Przykład 14.1
% Podział sygnału na ramki
[x,fpr,nbit]=wavread('AF1K1_3_.WAV'); % wczytanie sygnału dźwiękowego do wektora x
twind = 30; % długość okna obserwacji (ramki danych) w
% milisekundach
tstep = 15; % przesunięcie pomiędzy kolejnymi
% położeniami okna w milisekundach
n = (twind*0.001)*fpr; % liczba próbek okna (ramki danych)
m = (tstep*0.001)*fpr; % liczba próbek przesunięcia pomiędzy
% kolejnymi położeniami okna
ramkowanie(x,hamming(n),m); % wywołanie funkcji dzielącej sygnał w wektorze x na ramki.
% wywołanie funkcji bez argumentu wyjściowego
% powoduje wyświetlenie pierwszych dziesięciu
% ramek
ramkowanie.m
function f=ramkowanie(x,win,inc)
%RAMKOWANIE dzieli sygnał na pokrywające się ramki. Ramki są zapisywane w rzędach macierzy f.
%
% Argumenty wejściowe: x sygnał wejściowy
% win dyskretne okno czasowe np. okno Hamminga
% inc przesunięcie między początkami ramek (okien) w próbkach
%
% Argumenty wyjściowe: f macierz z sygnałem podzielonym na ramki w rzędach są zapisywane famki sygnału
% wymnożone przez funkcję okna, a indeks kolumny jest indeksem czasu
nx=length(x(:)); % długość sygnału w próbkach
lw =length(win); % długość okna (ramki) w próbkach
w = win(:)';
nli=nx-lw+inc;
nf = fix((nli)/inc); % liczba ramek w sygnale
na=nli-inc*nf;
fn=zeros(nf,lw);
indf= inc*(0:(nf-1)).'; % wektor z początkami ramek
inds = (1:lw); % wektor z indeksami próbek w ramce
fn(:) = x(indf(:,ones(1,lw))+inds(ones(nf,1),:)); % podział sygnału na ramki
f = fn .* w(ones(nf,1),:); % wymnażanie ramek sygnału przez funkcję okna
% obejrzyj pierwsze 10 ramek sygnału
if nargout<1
for i=1:10
subplot(3,1,1);
plot(fn(i,:)); xlabel('Indeks ramki'); title('Sygnał w ramce');
subplot(3,1,2)
plot(w); xlabel('Indeks ramki'); title('Funkcja okna');
subplot(3,1,3)
plot(f(i,:)); xlabel('Indeks ramki'); title('Zokienkowany sygnał w ramce');
pause
end
end
Przykład 14.2
Przykład (plik przyklad_14_02.m) zawiera wywołanie funkcji "melbankf" z programu 14.2 przedstawionego na stronie 513 oraz na listingu poniżej. Funkcja ma za zadanie obliczać wagi filtrów trójkątnych wykorzystywanych do obliczania współczynników MFCC. Funkcja melbankf bez argumentu wyjściowego wyświetla graficznie wagi filtrów, natomiast ta sama funkcja z podanym argumentem wyjściowym zwraca te wagi i jest używana w funkcji "melcepstrum".% Przykład 14.2
% Obliczanie współczynników filtrów trójk±tnych i ich graficzna
% interpretacja
J=10; % liczba filtrów (okien częstotliwościwoych) w banku filtrów
% pokrywających się na połowie szerokości w skali melowej
n=512; % długość FFT
fpr=16000; % fs częstotliwość próbkowania w Hz
fl=0; % fl Częstotliwość odpowiadająca początkowi pasma najniższego
% filtru, wyrażona jako ułamek fs,
fh=0.5 % fh Częstotliwość odpowiadająca końcowi pasma najwyższego
% filtru, wyrażona jako ułamek fs,
melbankf(J,n,fpr,fl,fh); % funkcja obliczająca współczynniki filtrów
% trójkątnych. Wywołanie funkcji bez
% argumentów wyjściowych powoduje graficzne
% pokazanie współczynników filtrów na
% wykresie.
melbankf.m
function [G,mn,mx]=melbankf(J,n,fs,fl,fh)
%MELBANKF oblicza wagi filtrów trójkątnych dla banku filtrów z melową skalą
% częstotliwości. Wagi są zapisywane w macierzy rzadkiej G
%
% Argumenty wejściowe:
% J liczba filtrów (okien częstotliwościwoych) w banku filtrów
% pokrywających się na połowie szerokości w skali melowej
% n długość fft
% fs częstotliwość próbkowania w Hz
% fl Częstotliwość odpowiadająca początkowi pasma najniższego
% filtru, wyrażona jako ułamek fs, najczęściej przyjmuje się fl=0
% fh Częstotliwość odpowiadająca końcowi pasma najwyższego
% filtru, wyrażona jako ułamek fs, z reguły przyjmuje się fh=0.5
%
% Argumenty wyjściowe:
% G macierz rzadka zawierająca wagi filtrów trójkątnych. Rozmiar
% macierzy G wynosi: size(G)=[J,mx-mn+1]
% mn indeks pierwszego niezerowego prążka FFT
% mx indeks ostatniego niezerowego prążka FFT
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%-------------------------------------------------------------------------
% OBLICZANIE WAG FILTRÓW
%-------------------------------------------------------------------------
mflh=[fl fh]*fs; % częstotliwości graniczne w Hz
mflh=frq2mel(mflh); % konwersja częstotliwości granicznych do skali
% melowej. Kod funkcji frq2mel znajduje się w pliku
% frq2mel.m znajdującym się na stronie internetowej:
% http://www.teledsp.kt.agh.edu.pl/rozdzial_14
melrng=mflh(2)-mflh(1);% zakres częstotliwości w melach
fn2=floor(n/2); % indeks prążka FFT o największej dodatniej
% częstotliwości - częstotliwość Nyquista jeżeli n
% jest parzyste
melinc=melrng/(J+1); % odstęp między środkowymi częstotliwościami filtrów
% w melach
% Obliczanie 4 prążków FFT odpowiadających następującym częstotliwościom
% blim = [początek_1-go_fitru środek_1-go_fitru środek_p-tego_fitru
% koniec_p-tego_filtru]. Kod funkcji mel2freq znajduje się w pliku
% mel2freq.m znajdującym się na stronie internetowej:
% http://www.teledsp.kt.agh.edu.pl/rozdzial_14
blim=mel2frq(mflh(1)+[0 1 J J+1]*melinc)*n/fs;
mc=mflh(1)+(1:J)*melinc; % środkowe częstotliwości filtrów w melach
b1=floor(blim(1))+1; % nr wymaganego do obliczeń prążka FFT o
% najmniejszym indeksie (może być ujemny)
b4=min(fn2,ceil(blim(4))-1); % nr wymaganego do obliczeń prążka FFT o
% największym indeksie
% rozmieszczanie wszystkich użytecznych prążków FFT względem środków
% filtrów
pf=(frq2mel((b1:b4)*fs/n)-mflh(1))/melinc;
% usuwanie wszystkich nieprawidłowych wartości z wektora pf, które mogą
% się pojawić na skutek błędów zaokrągleń
if pf(1)<0 pf(1)=[]; b1=b1+1; end
if pf(end)>=J+1 pf(end)=[]; b4=b4-1; end
fp=floor(pf); % Wektor na podstawie którego można przypisać
% prążki FFT do filtrów, prążek i należy do
% filtrów o numerach fp(1+i-b1)+[0 1]
pm=pf-fp; % mnożniki dla filtru o wyższym numerze dla danego
% prążka FFT
k2=find(fp>0,1); % k2+b1 jest indeksem pierwszego prążka FFT,
% należącego do dwóch filtrów
k3=find(fp<J,1,'last'); % k2+b1 jest indeksem ostatniego prążka FFT,
% należącego do dwóch filtrów
k4=numel(fp); % k4+b1 jest ostatnim prążkiem, który należy do
% jakiegokolwiek filtru
if isempty(k2)
k2=k4+1;
end
if isempty(k3)
k3=0;
end
r=[1+fp(1:k3) fp(k2:k4)];% indeksy rzędów w macierzy rzadkiej G. Indeks
% rzędu odpowiada numerowi filtra trójkątnego
c=[1:k3 k2:k4]; % indeksy kolumn w macierzy rzadkiej G.
% Indeks kolumny odpowiada numerowi prążka FFT
v=[pm(1:k3) 1-pm(k2:k4)];% Wartości wag filtrów trójkątnych zapisywanych
% w macierzy rzadkiej G. Do komórki macierzy G o
% indeksach [r(k),c(k)] jest zapisywana wartość
% v(k) -> G(r(k),c(k))=v(k).
mn=b1+1; % obliczenie indeksu pierwszego niezerowego prążka FFT
mx=b4+1; % obliczenie indeksu ostatniego niezerowego prążka FFT
if b1<0
c=abs(c+b1-1)-b1+1; % konwersja ujemnych indeksów prążków na dodatnie
end
G=sparse(r,c,v); % zapisywanie wag filtrów do macierzy rzadkiej
%-------------------------------------------------------------------------
% WIZUALIZACJA
%-------------------------------------------------------------------------
% rysowanie wykresów wag filtrów z banku filtrów w przypadku braku
% argumentow wyjściowych
%
if ~nargout % rysowanie idealizowanych ch-k częstotliwościowych filtrów
ng=201; % 201 punktów
me=mflh(1)+(0:J+1)'*melinc;
fe=mel2frq(me); % definiowanie częstotliwości
xg=mel2frq(repmat(linspace(0,1,ng),J,1).*repmat(me(3:end)-me(1:end-2),...
1,ng)+repmat(me(1:end-2),1,ng));
v=1-abs(linspace(-1,1,ng));
v=repmat(v,J,1);
plot(xg',v','b');
set(gca,'xlim',[fe(1) fe(end)]);
xlabel(['Częstotliwo¶ć [Hz]']);
pause
end
Przykład 14.3
Zawarta w przykładzie (plik przyklad_14_03.m) funkcja "melcepstrum", której kod przedstawiono poniżej oraz w programie 14.3 ze strony 517, dokonuje obliczania współczynników MFCC w ramkach sygnału wejściowego. Funkcja "melcepstrum" korzysta z przedstawionych wcześniej funkcji "ramkowanie" oraz "melbankf". Ponadto do obliczania współczynników została również użyta procedura "rdct" wykonująca dyskretną transformację kosinusową na zlogarytmowanych współczynnikach widmowych. Wywołanie funkcji melcepstrum bez argumentów wyjściowych powoduje wizualizację macierzy współczynników MFCC reprezentującej przetwarzane słowo. Natomiast wywołania funkcji "melcepstrum" z argumentami wyjściowymi zostały użyte w funkcjach "rozpoznawanie_DTW_melcepstrum" i "rozpoznawanie_HMM_melcepstrum". Dla większej przejrzystości programu nie uwzględniono w nim procedur do obliczania pierwszej i drugiej pochodnej.% Przykład 14.3
% Obliczanie współczynników MFCC i ich wizualizacja
[s,fs,nbit]=wavread('AF1K1_3_.WAV'); % wczytanie pliku WAV
% s - wektor z sygnałem
% fs - częstotliwość próbkowania
w='0'; % Jeżeli w='0', to został uwzględniony zerowy współczynnik MFCC
nc=12; % liczba współczynników MFCC z wyłączeniem zerowego współczynnika
p=30; % liczba filtrów w banku filtrów
n=480; % długość ramki w próbkach
inc=160; % przesunięcie między początkami ramek w próbkach
fl=0; % dolny kraniec pierwszego filtra z banku filtrów wyrażony jako
% ułamek fs
fh=0.5; % górny kraniec ostatniego filtra z banu filtrów wyrażony jako
% ułamek fs
melcepstrum(s,fs,w,nc,p,n,inc,fl,fh); % obliczanie współczynników MFCC.
% Wywołanie funkcji melcepstrum bez
% argumentów wyjściowych powoduje
% wizualizację macierzy
% współczynników MFCC reprezentującej
% przetwarzane słowo. Dla większej
% przejrzystości programu nie
% uwzględniono w nim procedur do
% obliczania pierwszej i drugiej
% pochodnej.
melcepstrum.m
function c=melcepstrum(s,fs,w,nc,p,n,inc,fl,fh)
%MELCEPSTRUM Oblicza melowe współczynniki cepstralne (MFCC) sygnału,
% okienkowanie sygnału w dziedzinie czasu jest dokonywane za pomocą
% okna Hamminga. Przyjęty został trójkątny kształt filtrów z banku
% filtrów w dziedzinie melowej.
%
% argumenty wejściowe:
% s sygnał mowy
% fs częstotliwość próbkowania w Hz
% nc liczba współczynników MFCC z wyłączeniem zerowego współczynnika
% n długość ramki w próbkach
% p liczba filtrów w banku filtrów
% inc przesunięcie między początkami ramek w próbkach
% fl dolny kraniec pierwszego filtra z banku filtrów wyrażony jako
% ułamek fs
% fh górny kraniec ostatniego filtra z banu filtrów wyrażony jako
% ułamek fs
%
% w jeżeli jest ustawiona wartość w='0', to zostaje uwzględniony
% zerowy współczynnik cepstralny
%
% argumenty wyjściowe:
% c macierz współczynników MFCC, indeks rzędu macierzy odpowiada
% indeksowi ramki
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
z=ramkowanie(s,hamming(n),inc); %podział sygnału na ramki z użyciem okna
% Hamminga
f=rfft(z.'); % obliczanie fft we wszystkich ramkach
% sygnału. Kod funkcji rfft znajduje się w
% pliku rfft.m znajdującym się na stronie
% internetowej:
% http://www.teledsp.kt.agh.edu.pl/rozdzial_14
[m,a,b]=melbankf(p,n,fs,fl,fh);% wywołanie funkcji obliczającej wagi
% filtrów (trójkątnych)wagi są zapisywane
% w macierzy rzadkiej m
pw=f(a:b,:).*conj(f(a:b,:)); % obliczanie kwadratu widma amplitudowego
% dla wszystkich ramek sygnału
pth=max(pw(:))*1E-20; % obliczanie minimalnej wartości jaka może
% być logarytmowana
y=log(max(m*pw,pth)); % obliczanie zlogarytmowanych
% współczynników widmowych.
c=rdct(y).'; % wykonanie dct na zlogarytmowanych współczynnikach widmowych
% Kod funkcji rdct znajduje się w pliku rdct.m znajdującym się
% na stronie internetowej:
% http://www.teledsp.kt.agh.edu.pl/rozdzial_14
nf=size(c,1); % liczba ramek sygnalu
nc=nc+1;
c(:,nc+1:end)=[];
if ~any(w=='0') % w przypadku, gdy nie jest uwzględniany zerowy
% melowy współczynnik cepstralny
c(:,1)=[];
nc=nc-1;
end
%-------------------------------------------------------------------------
% WIZUALIZACJA
%-------------------------------------------------------------------------
% Wizualizacja współczynników MFCC w przypadku braku argumentów wyj¶ciowych
if nargout<1
[nf,nc]=size(c);
t=((0:nf-1)*inc+(n-1)/2)/fs;
ci=(1:nc)-any(w=='0');
imh = imagesc(t,ci,c.');
axis('xy');
xlabel('Czas (s)');
ylabel('Numer współczynnika MFCC');
map = (0:63)'/63;
colormap([map map map]);
colorbar;
end
Przykład 14.4
Przykład (plik rozpoznawanie_DTW_melcepstrum.m) przedstawia procedurę na rozpoznawanie izolowanych słów za pomocą procedury DTW, gdzie jako cech sygnału mowy użyto współczynników MFCC. Procedurę uruchamia się za pomocą skryptu run (przykład 14.7), która oprócz rozpoznawania za pomocą DTW uruchamia również procedury do trenowania modeli HMM i do rozpoznawania mowy za ich pomocą. Procedura rozpoznająca izolowane słowa za pomocą metody DTW, oprócz skuteczności rozpoznawania oraz czasu rozpoznawania zwraca również tzw. macierz pomyłek (ang. confusion matrix). Indeks wiersza tej macierzy odpowiada numerowi klasy do rozpoznania (np. wiersz o indeksie 2 odpowiada numerowi klasy 2, czyli w naszym przypadku jest to słowo "jeden"). Natomiast indeks kolumny odpowiada numerowi rozpoznanej klasy. Z macierzy pomyłek można dowiedzieć się, ile słów zostało błędnie rozpoznanych. Przykładowo patrząc na wiersz 6 macierzy pomyłek uzyskujemy informacje, że 9 słów klasy 6 (słowo "pięć") zostało rozpoznanych dobrze (jako słowo " pięć"), a jedno słowo klasy 6 zostało rozpoznane jako należące do klasy 10 (czyli słowo "dziewięć"). Sama procedura DTW przedstawiona jest w programie 14.4 w książce na stronie 522 oraz jest wylistowana poniżej (plik dtw.m).% Przykład 14.4
% Rozpoznawanie izolowanych cyfr od 0 do 9 za pomoca DTW
function [] = rozpoznawanie_DTW_melcepstrum(experiment_name)
load(strcat(experiment_name, '/uczacy_melcepstrum.mat'));
train_feats = melfeats;
load(strcat(experiment_name, '/testowy_melcepstrum.mat'));
test_feats = melfeats;
[num_classes, num_examples] = size(train_feats);
[num_classes_test, num_examples_test] = size(test_feats);
confusion_matrix = zeros(num_classes, num_classes); % macierz pomylek
if num_classes_test ~= num_classes
error('Zbiór testowy zawiera przykłady opisujące niewidzane na etapie trenowania klasy')
end
tic();
% znajdz najmniejsza odleglosc DTW pomiedzy wszystkimi przykladami trenujacymi i testowymi.
% Klasa przykladu testowego jest zdefiniowana przez klase do ktorej nalezy przyklad uczacy
% majacy najmniejsza odleglosc DTW (bedacy najbardziej podobny)
for i=1:num_classes_test
for j=1:num_examples_test
test_wzr = cell2mat(test_feats(i,j));
if is_empty_cell(test_wzr) continue; end
min_DTW = Inf; class_idx = -1;
for k=1:num_classes
for l=1:num_examples
tr_wzr = cell2mat(train_feats(k,l));
if is_empty_cell(tr_wzr) continue; end
Cx = test_wzr; Cwzr = tr_wzr; Nwzr = [size(Cwzr,1)];
[D_DTW nr] = dtw( Cx, Cwzr, Nwzr ); % oblicz najmniejszą odległość między
% rozpoznawanym słowem a
% wzorcami słów pierwszej klasy
if D_DTW < min_DTW
min_DTW = D_DTW;
class_idx = k;
end
end
end
confusion_matrix(i, class_idx) = confusion_matrix(i, class_idx) + 1;
end
end
confusion_matrix
accuracy = sum(diag(confusion_matrix))./sum(sum(confusion_matrix)).*100;
fprintf(1, 'Skuteczność rozpoznawania za pomocą algorytmu DTW na zbiorze testowym wynosi %f %%\n', accuracy);
fprintf(1, 'Rozpoznawanie za pomocą metody DTW zajęło %f sekund.\n', toc());
%zapisz rezultaty na dysku
save_filename = strcat(experiment_name, '/', 'dtw_confusion_matrix.mat' );
if is_octave()
save('-mat', save_filename, 'confusion_matrix')
else
save(save_filename, 'confusion_matrix', '-mat')
end
dtw.m
% Implementacja metody DTW w programie MATLAB, przykład na str. 522
function [xxx nr ] = dtw( Cx, Cwzr, Nwzr)
% rozpoznawanie słowa poprzez porównanie macierzy jego wsp. cepstralnych Cx z wzorcami Cwzr metodą DTW
[ Ns, Np ] = size(Cx); % liczba wektorów współczynników cepstralnych słowa, ich długość
for numer = 1 : length(Nwzr) % porównaj Cx z Cwzr wszystkich wzorców
% Obliczenie odległości d(ns, nw) pomiędzy poszczególnymi cepstrami sygnału (ns) i sprawdzanego wzorca (nw)
Nw = Nwzr( numer ); % liczba wektorów wsp. cepstralnych wzorca
Q = round( 0.2 * max(Ns,Nw) ); % współczynnik szerokości ścieżki
d = Inf*ones(Ns,Nw); tg=(Nw-Q)/(Ns-Q); % inicjalizacja macierzy odległości, tangens kąta
for ns = 1:Ns % dla każdego cepstrum rozpoznawanego słowa
down(ns) = max( 1, floor(tg*ns-Q*tg)); % ograniczenie dolne
up(ns) = min( Nw, ceil(tg*ns+Q)); % ograniczenie górne
for nw = down(ns) : up(ns) % dla każdego cepstrum wzorca
d(ns,nw) = sqrt( sum((Cx(ns, 1:Np) - Cwzr(nw, 1:Np, numer)).^2 )); % odległość
end
end
% Obliczenie odległości zakumulowanej g()
g = d; % inicjalizacja
for ns = 2:Ns, g(ns,1) = g(ns-1,1) + d(ns,1); end % zakumuluj pierwszą kolumnę
for nw = 2:Nw, g(1,nw) = g(1,nw-1) + d(1,nw); end % zakumuluj pierwszy wiersz
for ns = 2:Ns % akumuluj w pionie (słowo)
for nw = max( down(ns), 2 ) : up(ns) % akumuluj w poziomie (wzorzec)
dd = d(ns,nw); % odległość cepstrum "ns" słowa od wzorca "nw"
temp(1) = g(ns-1,nw) + dd; % ruch do góry
temp(2) = g(ns-1,nw-1) + 2*dd; % ruch po przekątnej (do góry w prawo)
temp(3) = g(ns,nw-1) + dd; % ruch w prawo
g(ns,nw) = min( temp ); % wybierz minimalną wartość zakumulowaną
end
end
glob(numer) = g(Ns,Nw)/sqrt(Ns^2+Nw^2); % wartość zakumulowana "najkrótszej" ścieżki
end
[ xxx nr ] = min( glob ); % numer wzorca o najmniejszej wartości zakumulowanej
Przykład 14.5
Przykład (plik trenuj_HMM_melcepstrum.m) przedstawia procedurę trenowania ukrytych modeli Markowa z wykorzystaniem współczynników MFCC jako cech sygnału mowy. Procedurę trenowania ukrytych modeli Markowa można wywołać uruchamiając skrypt run.m (Przykład 14.7), z którego jest również uruchamiana procedura rozpoznawania za pomocą metody DTW oraz rozpoznawania za pomocą metody HMM. W procedurze trenowania ukrytych modeli Markowa jest wykorzytany algorytm EM (plik hmm_em.m przedstawiony poniżej). Z kolei algorytm EM korzysta z procedury obliczania macierzy prawdopodobieństw obserwacji (plik compute_pdfs.m) oraz z algorytmu "forward-backward" (plik hmm_fwdback.m przedstawiony poniżej oraz w książce w programie 14.5 na stronie 532).% Przykład 14.05
% Trenowanie izolowanych modeli HMM dla cyfr od 0 do 9
function [] = trenuj_HMM_melcepstrum(experiment_name)
% zakladamy, ze sygnal zostal zparametryzowany i zapisany w
% katalogu eksperymentu
load(strcat(experiment_name, '/uczacy_melcepstrum.mat'));
uczacy_feats = melfeats;
[num_classes, num_examples_train] = size(uczacy_feats);
F=3; % liczba stanow modelu HMM
D=size(uczacy_feats{1}, 2); % wymiarowosc cech
priors = [1 zeros(1, F-1)]';% prawdopodobienstwa poczatkowe
maxiter = 10; % maksymalna liczba iteracji EM
tic();
% utworz zbior num_classes modeli HMM (kazdy opisuje pojedyncza cyfre)
% majacych F stanow modelowanych przez pojedynczy komponent Gaussa o D wymiarach
hmms = hmm_create(num_classes, F, 1, D);
% zainicjalizuj poczatkowe wartosci rozkladow srednia i wariancja danej klasy
for i=1:num_classes
fprintf(1,'---- Initialising an HMM for the word %d \n', i-1)
[gmean, gstdev] = compute_cmvn_stats(melfeats(i,:));
hmms(i) = hmm_init(gmean, gstdev, hmms(i));
end
% estymuj parametry kazdego z modeli HMM za pomoca algorytmu EM
for i=1:num_classes
fprintf(1,'---- Training the HMM for the word %d \n', i-1)
[hmms(i), llike] = hmm_em(uczacy_feats(i,:), priors, hmms(i), maxiter);
end
fprintf(1,'Training of HMMs took %f seconds.\n', toc());
%zapisz modele na dysku
save_filename = strcat(experiment_name, '/', 'hmm_models.mat' );
if is_octave()
save('-mat', save_filename, 'hmms')
else
save(save_filename, 'hmms', '-mat')
end
end
hmm_em.m
function [hmm, llikes] = hmm_em(cfeats, priors, hmm_init, maxiter, update_mu, update_sigmas, update_transitions)
% Argumenty wejsciowe
% cfeats - cell array [1xR] zawierajaca R przykladow uczacych opisujacych dana cyfre
% priors - poczatkowe prawdopodobienstwa wejsc do stanow
% hmm_init - poczatkowy model HMM
% maxiter - maksymalna liczba iteracji algorytmu EM
% update_mu [domyslnie true] - aktualizowac wektory srednich?
% update_sigmas [domyslnie true]- aktualizowac macierze kowariancji?
% update_transitions [domyslnie true] - aktualizowac macierz przejsc?
%
% Wartosci zwracane:
% hmm - wysetymowany model HMM
% llikes - historia log-prawdopodobienstwa w poszczegolnych iteracjachifififi
llike_min_improvement = .5; % minimalna poprawa pomiedzy krokami E i M aby kontynuowac uczenie
llikes = [-Inf]; % log-prawdpodobienstwa akustyczne w kolejnych iteracjach
converged = false;
if nargin < 5 update_mu = true; end
if nargin < 6 update_sigmas = true; end
if nargin < 7 update_transitions = true; end
[xxx, R] = size(cfeats);
[F, M, D] = size(hmm_init.mu);
% inicjalizacja akumulatorow zbierajacych statystyki kroku E (uzyte w kroku M do aktualizacji parametrow)
tot_gammas = zeros(F, 1);
accs_mu = zeros(F, D);
accs_sigma = zeros(F, D, D);
accs_trans = zeros(F, F);
hmm = hmm_init; iter = 0;
while ( (iter < maxiter) && ~converged )
llike = 0.0;
RR = R; % tablica cfeats moze zawierac puste komorki, wtedy
% RR jest aktualizowane w przypadku opuszczenia przykladu
%Krok E
for r=1:R % oblicz statysytki dla kazdej z r wypowiedzi w zbiorze uczacym
o = cell2mat(cfeats(1,r));
if is_empty_cell(o) RR=RR-1; continue; end
T = size(o, 1);
% oblicz prawdopodobienstwa akustyczne z jakim aktualny model hmm
% poszczegolne ramki obserwacji wypowiedzi r
B = compute_pdfs(o, hmm);
% oblicz prawdopodobienstwa w przod i w tyl oraz prawdopodobienstwo gamma
[llike_r, alphas, betas, gammas] = hmm_fwdback(priors, B, hmm.A, true);
% zakumuluj statystyki do aktualizacji parametrow modelu w kroku M
tot_gammas = tot_gammas + sum(gammas, 2); % sum_{t=1}^Tr \gamma_j
if update_mu % dla wektorow srednich
accs_mu = accs_mu + gammas*o; % sum_{r=1}^R sum_{t=1}^T \gamma_f*o^r
end
if update_sigmas % dla macierzy kowariancji
for f=1:F
mu = reshape(hmm.mu(f,1,:), 1, D); w = zeros(D, D);
for t=1:T
w = w + gammas(f,t).*(o(t,:)-mu)'*(o(t,:)-mu); %sum_{t=1}^T \gamma_f*(o^r-mu)(o^r-mu)'
end
accs_sigma(f,:,:) = accs_sigma(f,:,:) + reshape(w, 1, D, D);
end
end
if update_transitions % dla macierzy przejsc
for t=1:T-1
accs_trans = accs_trans + hmm.A.*(alphas(:,t)*(B(:,t+1).*betas(:,t+1))');
end
end
llike = llike+llike_r;
end
llike = llike./RR;
fprintf(1, '------ Iteration %d : Log-likelihood is %g.\n', iter, llike);
%Krok M
if update_mu
new_mu = accs_mu./repmat(tot_gammas, 1, D); % rown. 14.52
hmm.mu = reshape(new_mu, F, M, D);
end
if update_sigmas
for f=1:F
new_sigma = reshape(accs_sigma(f,:,:), D , D)./repmat(tot_gammas(f), D, D); %rown. 14.53
hmm.sigma(f,1,:,:) = diag(diag(new_sigma)); % zakladamy ze estymujemy tylko wariancje
end
end
if update_transitions
hmm.A = normalise_matrix(accs_trans);
end
llike_improvement = abs(llike - llikes(end));
converged = (llike_improvement < llike_min_improvement);
iter = iter + 1;
llikes = [llikes llike];
end
end
compute_pdfs.m
function [B] = compute_pdfs(o, hmm)
% Argumenty wejsciowe
% o - macierz obserwacji [T x D]
% hmm - ukryty model Markowa
% Wartosci zwracane:
% B - macierz [FxT] zawierajÄ…ca prawdopodobienstwa akustyczne p(o(t) | q(t) = f; hmm)
[T, D] = size(o);
[F, M, D2, D2_] = size(hmm.sigma);
B = zeros(F, T);
if (D~=D2)
error('Dimensions mismatch!.')
end
% oblicz niezalezne od danych wartosci komponentow
% aby nie robic tego samego dla kazdej z T obserwacji osobno
g_const = zeros(F,M); inv_sigmas=zeros(F,M,D,D);
for f=1:F
for m=1:M
sigma = reshape(hmm.sigma(f,m,:,:),D,D);
g_const(f,m) = 1./sqrt((2*pi).^D*det(sigma));
inv_sigmas(f,m,:,:) = inv(sigma);
end
end
for t=1:T %dla kazdej ramki
for f=1:F %dla kazdego ze stanow HMM
for m=1:M %dla kazdego komponentu Gaussa w mieszaninie
mu = reshape(hmm.mu(f,m,:),1,D);
%sigma = reshape(hmm.sigma(f,m,:,:),D,D);
inv_sigma = reshape(inv_sigmas(f,m,:,:),D,D);
exponent = (o(t,:)-mu)*inv_sigma*(o(t,:)-mu)';
B(f,t) = B(f,t) + hmm.c(m).*g_const(f,m).*exp(-0.5*exponent);
%B(f,t) = B(f,t) + hmm.c(m).*mvnpdf(o(t,:),mu,sigma);
end
end
end
end
hmm_fwdback
function [loglike, alphas, betas, gammas] = hmm_fwdback(priors, B, A, back_pass)
% Argumenty wejsciowe:
% priors - wektor prawd. poczatkowych P(q=j|t=1)
% B = {b_j(o(t))}: p(o(t)|q=s_j;M), [FxT] macierz prawdopodobienstw z
% jakim dany stan q generuje obserwacje o(t)
% A = {a_ij}: P(q+1|q) macierz przejsc [FxF] modelu HMM
% back_pass - obliczac prawd. wstecz? (pomijane na etapie rozpoznawania)
% Wartosci zwracane
% loglike - log p(O;M) - log prawdopodobienstwo calkowie z
% jakim dany model M generuje obserwacje O_1:T
% alphas - macierz [FxT] prawd. w przod
% betas - macierz [FxT] prawd. w tyl lub [] jezeli back_pass=false
% gammas = macierz [FxT] prawd. przynaleznosci stanu f do obserwacji o_t
if (nargin==3) back_pass = false; end
[F, T] = size(B);
loglike = -Inf;
alphas = zeros(F, T); ascales = zeros(1,T);
betas = []; bscales = [];
gammas = [];
% oblicz p(o(1),o(2),...,o(t),q;M)
alphas(:,1) = priors.*B(:,1);
[alphas(:,1) ascales(1)] = normalise(alphas(:,1));
for t=2:T
for f=1:F
alphas(f,t) = alphas(:,t-1)'*A(:,f)*B(f,t);
end
[alphas(:,t), ascales(t)] = normalise(alphas(:,t));
end
loglike = sum(log(ascales));
% obliczenie prawodpodobienstw backward jest opcjonalne
if (back_pass == false) return; end
betas = zeros(F, T); bscales = zeros(1, T);
% oblicz p(o(t+1),o(t+2),...,o(t),q;M)
betas(:,T) = ones(F, 1);
bscales(T) = 1.0;
for t=T-1:-1:1
for f=1:F
betas(f,t) = (betas(:,t+1).*B(:,t+1))'*A(:,f);
end
[betas(:,t), bscales(t)] = normalise(betas(:,t));
end
bscales(1) = sum((betas(:,1).*B(:,1))'*priors);
beta_1=sum(log(bscales));
%oblicz statystyki zajetosci stanow (rown. 14.51)
gammas = zeros(F, T);
for t=1:T
[gammas(:,t), xxx] = normalise(alphas(:,t).*betas(:,t));
end
end
Przykład 14.6
Przykład (plik rozpoznawanie_HMM_melcepstrum.m) przedstawia procedurę rozpoznawania izolowanych cyfr za pomocą ukrytych modeli Markowa z wykorzystaniem współczynników MFCC jako cech sygnału mowy. Procedurę można wywołać uruchamiając skrypt run.m (Przykład 14.7), z którego jest również uruchamiana procedura na rozpoznawanie za pomocą metody DTW oraz procedura trenowania ukrytych modeli Markowa. W procedurze jest wykorzytany algorytm Viterbiego (plik hmm_viterbi.m przedstawiony poniżej oraz w książce w programie 14.6 na str. 534-535).% Przykład 14.0X
% Rozpoznawanie izolowanych cyfr od 0 do 9 za pomoca modeli HMM
function [] = rozpoznawanie_HMM_melcepstrum(experiment_name)
load(strcat(experiment_name, '/testowy_melcepstrum.mat'));
load(strcat(experiment_name, '/hmm_models.mat'));
feats = melfeats;
[num_classes, num_examples] = size(feats);
confusion_matrix = zeros(num_classes, num_classes); % macierz pomylek
tic();
for i=1:num_classes %sprawdz kazdy model dla kazdego przypadku testowego
for j=1:num_examples
o = cell2mat(feats(i,j));
if is_empty_cell(o) continue; end
max_llike=-Inf;
mdl_idx=-1;
for k=1:num_classes
hmm = hmms(k); F = size(hmm.mu, 1); priors = [1 zeros(1, F-1)]';
B = compute_pdfs(o, hmm);
% oblicz prawdopodobienstwo najbardziej prawdopodobnej sciezki
[llike, best_path] = hmm_viterbi(priors, B, hmm.A);
% lub klasyfikuj na podstawie prawdopodobienstwa calkowitego
%[llike, aa, bb, gg] = hmm_fwdback(priors, B, hmm.A, false);
if llike > max_llike
max_llike = llike;
mdl_idx = k;
end
end
confusion_matrix(i, mdl_idx) = confusion_matrix(i, mdl_idx) + 1;
end
end
confusion_matrix
accuracy = sum(diag(confusion_matrix))./sum(sum(confusion_matrix)).*100;
fprintf(1, 'Accuracy on test set using HMMs is %f %%\n', accuracy);
fprintf(1, 'Recognition by HMMs took %f seconds.\n', toc());
%zapisz resultaty na dysku
save_filename = strcat(experiment_name, '/', 'hmm_confusion_matrix.mat' );
if is_octave()
save('-mat', save_filename, 'confusion_matrix')
else
save(save_filename, 'confusion_matrix', '-mat')
end
end
hmm_viterbi.m
function [llike, best_path] = hmm_viterbi(priors, B, A)
% Argumenty wejsciowe:
% priors - wektor prawd. poczatkowych P(q=j|t=1)
% B = {b_j(o_t)}: p(o_t|q_t=j;M), [QxT] macierz prawdopodobienstw z
% jakim dany stan j modelu M generuje obserwacje o_t
% A = {a_ij}: P(q+1|q) macierz przejsc [QxQ] modelu HMM
% Wartosci zwracane
% llike = log max_Q p(O,Q;M) - log prawdopodobienstwo Viterbiego
% najbardziej prawdopodobnej sekwencji stanow Q*
% best_path - argmax_Q p(O,Q;M) - wektor o rozmiarze T z indeksami stanow
[F, T]=size(B);
llike = -Inf;
deltas = zeros(F, T);
psi = zeros(F,T);
best_path = zeros(1, T);
scales = ones(1,T);
t=1;
deltas(:,t) = priors.*B(:,t);
[deltas(:,t), s] = normalise(deltas(:,t));
scales(t) = 1/s;
psi(:,t) = 0;
for t=2:T
for f=1:F
[deltas(f,t), psi(f,t)] = max(deltas(:,t-1).*A(:,f));
deltas(f,t) = deltas(f,t).*B(f,t);
end
[deltas(:,t), s] = normalise(deltas(:,t));
scales(t) = 1/s;
end
[like, best_path(T)] = max(deltas(:,T));
for t=T-1:-1:1 %sledzenie sciezki wstecz
best_path(t) = psi(best_path(t+1), t+1);
end
llike = -sum(log(scales));
end
Przykład 14.7
Przykład (plik run.m) uruchamia rozpoznawanie izolowanych cyfr zarówno za pomocą metody DTW, jak również za pomocą metody HMM. Dzięki temu można porównać rezultaty rozpoznawania za pomocą obydwu metod. UWAGA! Rozpoznawanie za pomocą klasyfikatora DTW może w zależności od komputera i oprogramowania zająć od kilku do kilkudziesięciu minut.% Skrypt glowny uruchamiajacy wykonanie nastepujacych etapow:
%
% 1) ekstrakcji cech melcepstralnych (zbior uczacy i testowy)
% 2) rozpoznawania za pomoca DTW
% 3) estymacji modeli HMM
% 4) rozpoznawania za pomoca HMM
%
% Przed uruchomieniem skryptu 'run' wejdz do katalogu w wybranym
% srodowisku uruchomieniowym {matlab, octave} przez wpisanie w
% oknie komend cd 'sciezka/do/rozdzial_14' (n.p. C:\teledsp\rozdzial_14)
%
% Koniecznym moze okazac sie rowniez dodanie ww. katalogu do sciezki wybranego
% srodowiska uruchomieniowego.
%inicjalizacja ogolnych ustawien srodowiska
clear all; clc; init_environment();
% ========================= Parametry eksperymentu ========================
experiment_name='exp'; %utworz i zapisz rezulataty w katalogu o takiej nazwie
mkdir(experiment_name);
%1) ========================================================================
% dla kazdego pliku z listy uczacy.list wykonaj parametryzacje zapisujac wyniki
% w katalogu o nazwie experiment_name/uczacy_melcepstrum.mat
disp('-- EKSTRAKCJA CECH ZBIORU UCZACEGO ')
ekstrahuj_melcepstrum('uczacy.list', experiment_name)
% dla kazdego pliku z listy testowy.list wykonaj parametryzacje zapisujac wyniki
% w katalogu o nazwie experiment_name/testowy_melcepstrum.mat
disp('-- EKSTRAKCJA CECH ZBIORU TESTOWEGO ')
ekstrahuj_melcepstrum('testowy.list', experiment_name)
%2) =========================================================================
disp('-- ROZPOZNAWANIE ZA POMOCA KLASYFIKATORA DTW ')
rozpoznawanie_DTW_melcepstrum(experiment_name)
%3) =========================================================================
disp('-- TRENOWANIE HMM ')
trenuj_HMM_melcepstrum(experiment_name)
%4) =========================================================================
disp('-- ROZPOZNAWANIE ZA POMOCA HMM ')
rozpoznawanie_HMM_melcepstrum(experiment_name)
Osoby, które chciałyby bardziej dogłębnie zapoznać się z tematyką rozpoznawania mowy zachęca się do zapoznania z pakietem HTK (htk.eng.cam.ac.uk) oraz z projektem Kaldi (kaldi.sourceforge.net).