2.2 Odhad koeficientů AR modelu

2.2.1 Odhad AR koeficientů obecně (pro nestacionární proces)

Náhodný proces $x[n]$ je AR (autoregressní) proces řádu $M$, jestliže platí

$$x[n] = -\sum_{m=1}^{M} a_m x[n-m] + v[n],$$ (2.1)

kde $v[n]$ je stacionární bílý gusovský proces s nulovou střední hodnotu a variancí $\sigma^2$. Typicky se nazývá excitační nebo budící

$$v[n] \sim \mathcal{N}(0,\sigma^2).$$ (2.2)

Koeficienty $a_1,\ldots,a_M$ se nazývají AR koeficienty. Rovnice vlastně říkají jak je proces $x[n]$ pomocí budícího procesu $v[n]$ atd. namodelován. Proto se používá také označení AR model. Napříkad lze říci, že proces $x[n]$ je generován AR modelem řádu $M$, lze mluvit o odhadu koeficientů AR modelu atd.. Takový AR model je znázorněn na obrázku Obr. 2.2 (nakreslili jsme pouze filtr posaný rovnicí (2.1)).

Figure 2.2: Modelování AR procesu

Určitě jste už řadu AR procesů potkali. Například proces $x[n]$ na kterém jste zkoušeli odhady momentů v prvním cvičení byl po částech AR řádu 1.

Nyní předpokládejme, že máme k dispozici realizaci procesu $x[n]$ a známe řád AR procesu. Naším úkolem bude odhadnout jeho koeficienty případně $\sigma^2$.

Zde stručně uvedeme odhad metodou nejmenších čtverců, který byl uveden (podrobně) na přednášce ^2.2. Nejprve nebudu předpokládat nic o stacionaritě případně ergodicitě $x[n]$, tedy $x[n]$ může být nyní nestacionární (AR koeficienty případně variance budícího procesu $\sigma^2$ se mohou měnit v čase).

Rovnice (2.1) vyjádříme budící proces

$$v[n] = x[n] + \sum_{m=1}^{M} a_m x[n-m].$$ (2.3)

Náš odhad je založen na minimalizaci disperze procesu $v[n]$. Tedy AR koeficienty odhadneme tak, aby byl minimalizován moment $\mathrm{E}[v^2[n]]$. Kdybychom $a_m$ znali můžeme z (2.3) takto budící proces skutečně vypočítat. Bohužel budící proces $v[n]$ a koeficienty $a_m$ jsou nám nyní neznámé, jsou předmětem odhadu. Doslova $a_m$ měníme a hledáme jejich optimální hodnotu s tím se mění i $v[n]$. Aby bylo zřejmé že nyní $v[n]$ a $a_m$ mohou být odlišné od skutečných, označím je jinak. Místo $v[n]$ použiji označení $e[n]$ a místo $a_m$ použiji $a_m'$. Rovnice (2.3) s použitím nového značení nyní vypadá takto

$$e[n] = x[n] + \sum_{m=1}^{M} a_m' x[n-m],$$ (2.4)

viz. rovněž Obr. 2.3.

Figure 2.3: Analýza AR procesu

Filtr na obrázku se nazývá analyzující (odhadujeme parametry signálu, tedy analyzujeme signál). Na obrázku je uveden signál $x_p[n]$, nazveme ho predikce

$$x_p[n] = - \sum_{m=1}^{M} a_m' x[n-m].$$ (2.5)

S použitím $x_p[n]$ lze (2.4) zapsat

$$e[n] = x[n] - x_p[n].$$ (2.6)

Odtud název pro $e[n]$ chyba predikce. Dříve uvedené kritérium na odhad $a_m$ má nyní význam minimalizace disperze chyby predikce $e[n]$. Hledáme tedy $a_m'$ taková, aby byla minimalizována disperze $e[n]$ ^2.3.

Nyní již samotné odvození. Označme vektor koeficientů a vektor zpožděných vzorků signálu $x[n]$

$$\begin{gathered} \mathbf{a}^T = [ a_1, \ldots, a_M ], \\ \mathbf{x}^T[n-1] = [ x[n-1], \ldots, x[n-M] ]. \end{gathered}$$

Rovnici (2.6) můžeme nyní zapsat

$$e[n] = x[n] + \mathbf{a}^T \mathbf{x}[n-1].$$

Vyjádříme $e^2[n]$ a aplikujeme operátor $\mathrm{E}[.]$, dospějeme tak k disperzi

$$\begin{aligned}\mathrm{E}[e^2[n]] &= \mathrm{E}[x^2[n] + 2 x[n]... ...{r}_p + \mathbf{a}^T \mathbf{R}_{xx} \mathbf{a}, \\ \end{aligned}$$

kde jsme v posledním kroku označili autokorelační matici $\mathbf{R}_{xx}$ a vektor autokorelací $\mathbf{r}_p$

$$\begin{gathered}\mathbf{R}_{xx} = \mathrm{E}[ \mathbf{x}[n-1]... ... \mathbf{r}_p = \mathrm{E}[x[n]\mathbf{x}[n-1]]. \end{gathered}$$ (2.8)

Hledáme hodnoty $\mathbf{a}$ pro které je disperze (2.7) minimální. Zderivujeme tedy (2.7) podle složek vektoru $\mathbf{a}$ výsledky uspořádáme do sloupcového vektoru (gradient) ^2.4

$$\nabla_{\mathbf{a}} (\mathrm{E}[e^2[n]]) = \frac{\partial \mathrm... ...2[n]]}{\partial \mathbf{a}} = 2 \mathbf{r}_p + 2 \mathbf{R}_{xx} \mathbf{a} .$$

Nyní položíme gradient rovný $\mathbf{o}$ (nulový vektor) a dostáváme

$$\mathbf{a}_{opt} = - \mathbf{R}_{xx}^{-1}\mathbf{r}_p.$$ (2.9)

Pro disperzi chyby predikce v případě $\mathbf{a}=\mathbf{a}_{opt}$ dostáváme

$$\begin{aligned}\mathrm{E}[e^2[n]]_{min} &= \mathrm{E}[x^2[n]] +... ...athrm{E}[x^2[n]] + \mathbf{r}_p^T \mathbf{a}_{opt}. \end{aligned}$$

Z kritéria pro odhad zároveň plyne, že tato disperze je minimální odtud označení $\mathrm{E}[e^2[n]]_{min}$.

Nyní již víme jak z momentů $x[n]$ AR koeficienty odhadnout otázkou je jak odhadnout momenty $x[n]$. To jsme ale probíraly na minulém cvičení. Například víte-li, že signál lze považovat za stacionární a ergodický na úsecích nějaké fixní délky, můžete použít blokový odhad atd..

2.2.2 LPC vokodér

V následujících úlohách se stručně podíváme na problematiku komprese řečového signálu. Řečový signál lze zjednodušeně popsat jako AR proces jehož parametry se v čase mění viz. Obr. 2.2. Přenosu řečového signálu naznačuje Obr. 2.4

Figure 2.4: Přenos řečového signálu

Řečový signál putuje od mluvčího přes mikrofon A/D převod atd. na vstup analyzujícího filtru. Zpracování probíhá obvykle po blocích. V rámci každého bloku se odhadnou AR koeficienty. AR koeficienty spolu s informacemi o chybě predikce (např. $\sigma^2$) jdou na vstup přenosové linky. Na druhém konci se z AR koeficientů sestaví syntetizující filtr (inverzní k analyzujícímu) viz. Obr. 2.2 a z informací o chybě predikce se sestaví budící signál (například bílý šum s variancí $\sigma^2$). Přivedením budícího signálu na vstup syntetizujícího filtru obdržíme jeden blok syntetizovaného řečového signálu.

Stáhněte si nějaký řečový signál například tento a načtěte si jej do MATLABu

x = load("rec.asc");

Cvičení 2.1: Vaším úkolem je blokově odhadnout AR koeficienty (2.9) a varianci budícího procesu (2.10) pro řečový signál $x[n]$. Předpokládejte fixní řád AR modelu řeči $M=8$ a dále stacionaritu a ergodicitu řeči na úsecích délky $L=256$ (použijte bloky délky $L=256$).

Poznámka: K odhadu autokorelační matice použijte vnějšího součinu jako v (2.8)

x = x(:);       % udela z x sloupcovy vektor
N = length(x);
M = 10;         % rad
L = 256;        % delka bloku
B = floor(N/L); % pocet bloku

y = zeros(N,1); e = zeros(N,1); u = zeros(N,1);

for b = 1:B; %B %4
  % kodovani
  x_blok = x( (b-1)*L+1:L*b );
  range = 1+M:L;
  X = zeros(L-M,M+1);
  for k = 1:M+1
    X(:,k) = x_blok(range-k+1);
  end;
  Rxx_ = X'*X/(L-M);   % spocteme autokorelacni matici Rxx_ o rozmeru M+1
  Rxx = Rxx_(1:M,1:M); % vybereme autokorelacni matici Rxx o rozmeru M
  rp = Rxx_(2:M+1,1);  % a vektor rp
  a = -inv(Rxx)*rp;
  e_vykon = Rxx(1,1) + rp'*a; % vykon chyby predikce
end

Cvičení 2.2: Pro jeden blok^2.5Vyneste predikci $x_p[n]$ spolu s původním signálem $x[n]$ a chybou predikce $e[n]$.

  an_A = [1; zeros(M,1)]; % koeficienty analyzujícího filtru
  an_B = [1;a];           
  xp_blok = -filter([0;a],an_A,x_blok); % predikce
  e_blok = x_blok - xp_blok;            % chyba pridikce
                                        % e_blok = filter(an_B,an_A,x_blok);

Figure 2.5: Původní signál $x[n]$ a jeho predikce $x_p[n]$

Cvičení 2.3: Nyní se pokuste řeč syntetizovat. V každém bloku z odhadnutých AR koeficientů vytvořte syntetizující filtr. Tento filtr vypadá jako na Obr. 2.2 je inverzní k analyzujícímu. Jako jeho buzení $u[n]$ použijte:

1.

Chybou predikce $u[n] = e[n]$. Výstup synetizujícího filtru $y[n]$ by měl být v tomto případě identický s původním řečovým signálem.

2.

Bílý stacionání gausovský proces s nulovou střední hodnotou a variancí rovnou výkonu chyby predikce (2.10).

3.

Posloupnost periodicky se opakujících jednotkových pulsů s periodou $P = 38$ (Periodu odhadněte podle Vašeho signálu $x[n]$, stačí orientačně pohledem na modulové DFT spekrum).

Pro každý případ buzení porovnejte výsledek syntézy $y[n]$ s řečovým signálem $x[n]$.

Pro každý případ buzení si pro jeden blok vyneste: původní signál $x[n]$, výsledek syntézy $y[n]$, chybu predikce $e[n]$ a buzení syntetizujícího filtru $u[n]$.

Nevynášejte pouze časové průběhy, ale rovněž výkonová spektra v decibelech (lépe odhady PSD), rovněž si vyneste kvadrát modulu frekvenční charakteristiky syntetizujícího filtru v decibelech, abyste ho mohli porovnat s výkonovým spektrem $x[n]$, případně $y[n]$.

Při porovnávání výsledku syntézy $y[n]$ s původním řečovým signálem $x[n]$ je dobré si signály poslechnout (např. sound(y)), a vynést. Vyneste si rovněž spektrogramy obou signálů (např. specgram(y)^2.6), jde z nich vyčíst více než z časových průběhů.

Figure 2.6: Průběh $x[n]$, spektrogram $x[n]$, spektrogram $y[n]$ (syntetizující filtr je buzen směsí šumu a impulzů)

Skript pro syntézu $y[n]$ zde neuvádím. Zmíním se pouze o vynášení výkonových spekter.

figure
minp = 1e-3; 
W = linspace(0,2*pi,100); % frekvence
% vykonove spektrum v [dB]
plot((0:L-1)/L*2, 10*log10(abs(fft(x_blok)).^2/L+minp), 'b' );
hold on
% modul |H|^2 vynasobeny vykonem chyby predikce v [dB]
plot(W/pi,10*log10(abs(freqz(sy_B*sqrt(e_vykon),sy_A,W)).^2) , 'r');
hold off
title('|X|^2 (modre), |H|^2*sigma')
ylabel('[dB]')
xlabel('norm. frekv. (W/pi)')

Figure 2.7: Průběhy signálů $x[n]$, $y[n]$, $e[n]$, $u[n]$ a jejich výkonová spektra pro jeden blok, syntetizující filtr je buzen chybou predikce.

Figure 2.8: Průběhy signálů $x[n]$, $y[n]$, $e[n]$, $u[n]$ a jejich výkonová spektra pro jeden blok, syntetizující filtr je buzen bílým gausovským šumem s variancí rovnou výkonu chyby predikce.

Figure 2.9: Průběhy signálů $x[n]$, $y[n]$, $e[n]$, $u[n]$ a jejich výkonová spektra pro jeden blok, syntetizující filtr je buzen posloupností jednotkových pulsů.