3.2 LMS prediktor, a jeho použití k odhadu frekvence

K odhadu frekvence signálu $x[n]$ použijeme LMS prediktor 2. řádu jehož struktura je na Obr. 3.2

Figure 3.2: LMS prediktor, 2. řád

3.2.1 Odvození LMS algoritmu pro prediktor

Z Obr. 3.2 jednoduše sestavíme rovnici filtrace

$$e[n] = x[n] + w_1[n]*x[n-1] + w_2[n]*x[n-2].$$ (3.1)

Zde budu používat zápisu pomocí vektorů. Důvody jsou takové, že se jednodušeji zobecňuje na vyšší řády a hlavně chci, aby jste vektorové notaci trochu více přivykli (některé věci by se nám později bez ní vysvětlovali podstatně hůře). Definujeme tedy vektory vah a zpožděných vzorků

$$\begin{gathered} \mathbf{w}^T[n] = [ w_1[n], w_2[n] ], \\ \mathbf{x}^T[n-1] = [ x[n-1], x[n-2] ]. \end{gathered}$$

S použitím právě definovaných vektorů zapíšeme rovnici filtrace (3.1) takto

$$e[n] = x[n] + \mathbf{w}^T[n] \mathbf{x}[n-1].$$ (3.2)

Jak již víte, LMS nahrazuje kritérium minimalizace $\mathrm{E}[e^2[n]]$ (které jsme rozebírali například minulé cvičení) kritériem minimalizace aktuální kvadratické chyby $e^2[n]$. Tedy váhy upravuje tak, aby se $e^2[n]$ zmenšila. Ona úprava spočívá v odečtení nějakého malého násobku gradientu $e^2[n]$ od vah (pohne tedy vahami proti směru gradientu^3.2). Nejprve s použitím (3.2) spočteme gradient $e^2[n]$ vzhledem k $\mathbf{w}$

$$\nabla_{\mathbf{w}}e^2[n] = \frac{\partial e^2[n]}{\partial \math... ...} = 2 e[n] \frac{\partial e[n]}{\partial \mathbf{w}} = 2 e[n] \mathbf{x}[n-1].$$

Nyní napíšeme rovnice pro úpravu vah tak jak bylo popisováno ^3.3

$$\mathbf{w}[n+1] = \mathbf{w}[n] - \mu \nabla_{\mathbf{w}}e^2[n] = \mathbf{w}[n] - \mu e[n] \mathbf{x}[n-1].$$ (3.3)

Poznámka: Co se tady vlastně děje? Jaký má výsledek vztah k původnímu kritériu minimalizace $\mathrm{E}[e^2[n]]$. Za předpokladu stacionarity a ergodicity (alespoň na intervalech omezené délky) můžeme tento vztah skutečně najít. Rozepíšeme-li $K$ kroků rekurze (3.3) dojdeme ke vztahu mezi vahami $\mathbf{w}[n+1]$ a $\mathbf{w}[n-K]$

$$\mathbf{w}[n+1] = \mathbf{w}[n-K] - \mu \sum_{k=0}^{K} e[n-k] \mathbf{x}[n-k-1].$$

Člen s $\mu$ má stejný tvar jako gradient z výkonu

$$\frac{1}{K+1} \sum_{k=0}^{K} e^2[n-k],$$

až na konstantu, která se dá zase schovat do $\mu$. Tento výkon je vlastně odhadem $\mathrm{E}[e^2[n]]$ za předpokladu stacionarity a ergodicity $e[n]$ na $n-K,\ldots,n$ (vzpomeňte si, jak jsme odhadovali disperzi v 1. cvičení). Je-li tedy uvedený předpoklad splněn pro dostatečně velká $K$ (aby byl odhad pokud možno přesný) minimalizuje LMS i $\mathrm{E}[e^2[n]]$.

3.2.2 Vztah vah u prediktoru 2. řádu a frekvence vstupního signálu

Jelikož se LMS snaží minimalizovat výkon na výstupu (viz. poznámka výše), umisťuje nuly poblíž těch frekvencí $x[n]$, kde je soustředěn nejvyšší výkon (v případě náhodných procesů bychom měli mluvit o maximech spektrální výkonové hustoty). V případě harmonického signálu tedy na frekvenci harmonického signálu (nikde jinde výkon není). Z toho vychází myšlenka určení frekvence vstupního signálu $x[n]$ jako polohy nuly (přesně úhlu nuly, při komplexně sdružených nulách) filtru s fixními vahami $\mathbf{w} = \mathbf{w}[n]$ ^3.4

Předně se potřebujeme podívat jak vypadá přenosová funkce takového filtru. Použijeme tedy Z-transformaci na (3.1) (s náhradou $\mathbf{w}$ za $\mathbf{w}[n]$ - fixní váhy) a vyjádříme přenosovou funkci

$$H(z) = \frac{E(z)}{X(z)} = \frac{z^2 + w_1 z + w_2}{z^2}.$$ (3.4)

Při předpokladu komplexně sdružených nul ( $r e^{j \Theta}$, $r e^{-j \Theta}$) lze čitatel $H(z)$ rovněž zapsat

$$(z-r e^{j \Theta})(z-r e^{-j \Theta}) = z^2 - 2 r z \cos ( \Theta ) + r^2.$$ (3.5)

Porovnáním čitatele (3.4) s (3.5) získáme vztah mezi vahami a úhlem nuly (odhadovaná frekvence) a jejím modulem

$$\begin{gathered}w_1 = - 2 r \cos ( \Theta ), \\ w_2 = r^2. \end{gathered}$$ (3.6)

Tyto vztah jdou jednoduše invertovat

$$\begin{gathered}r = \sqrt{w_2}, \\ \Theta = \arccos( -\frac{w_1}{2 r} ), \end{gathered}$$ (3.7)

kde $\Theta$ je náš odhad frekvence vstupního signálu $x[n]$.

3.2.3 Vlastní cvičení

Cvičení 3.1: Implementujte v Matlabu strukturu LMS prediktoru 2. řádu, jako vstupní signál použijte $x[n]$ vygenerovaný v předchozím cvičení (začněte s A - konstantní úhlová frekvence). Vyneste si průběh vah. Konvergenční konstantu $\mu$ volte tak aby došlo k ustálení.

Řešení:

x = x(:); % do sloupcoveho vektoru
mu = 0.05; % konvergencni konstanta

M = 2; % rad (nelze menit, ale je to prehlednejsi, nez psat vsude 2)
w = zeros(M,N); % sloupcove vektory vah pro jednotlive casy

w(:,3) = [ 0; 0]; % nulove poc. podminky

for n = 3:N-1
    % filtrace
    xp = - w(:,n)'*x(n-1:-1:n-M); % predikce
    e = x(n) - xp; % chyba predikce
    % uprava vah
    w(:,n+1) = w(:,n) - mu*e*x(n-1:-1:n-M);
end

Cvičení 3.2: Doplňte smyčku LMS algoritmu o odhad frekvence vstupního signálu s použitím aktuálních hodnot vah. Frekvenci odhadněte jako polohu nuly viz. (3.7). Porovnejte odhad frekvence se skutečným průběhem frekvence signálu $x[n]$. U případů B a C nastavte $\mu$ tak, aby LMS stíhal sledovat změny frekvence.

Poznámka: Před použitím (3.7) nezapomeňte otestovat, zda jsou nuly opravdu komplexně sdružené, nebo násobné.

Vynesení výsledků:

figure(2)
subplot(2,2,1);
plot(1:N,Wchp/pi,'r'); % skutecna frekvence
hold on
plot(1:N,real(Wodhad)/pi,'b'); % odhadnuta frekvence (faze nuly)
hold off
subplot(2,2,3);
plot(1:N,real(rodhad),'b'); % modul nuly
subplot(2,2,2); % prubehy vah
plot(1:N,w(1,:),'b'); 
hold on
plot(1:N,w(2,:),'r');
hold off
subplot(2,2,4); % nuly v polrnich souradnicich
polar(real(Wodhad(3:N-1)),real(rodhad(3:N-1)),'b');

figure(3) % modul prenosove fce filtru
WN = 100;
W = linspace(0,pi,WN);
H = zeros(WN,N);
for n = 1:N % pro kazde n spocteme frekvencni charakteristiku
[H(:,n)] = freqz([1 w(:,n)'],[1 0 0],W)'; 
end
surf(1:N,W/pi,abs(H)) % vykresleni modulu
colormap(jet);
shading interp;
view([-37.5,70]);

Výsledky pro lineární chirp:

Figure 3.3: Průběhy vah a nul pro lineární chirp: (zleva doprava) úhel nuly (odhad frekvence) spolu se skutečnou frekvencí vstupního signálu na čase, časový průběh vah, modul nuly na čase, nuly (přes všechny časy) v polárních souřadnicích. Všimněte si, že odhad je za skutečnou hodnotou frekvence trochu opožděn.

Figure: Modulová frekvenční charakteristika filtru (představte si, že se váhy zafixovaly na hodnotě $\mathbf{w}[n]$) na čase pro lineární chirp

Cvičení 3.3: Toto cvičení je věnováno sledování vlivu počátečních podmínek na rychlost konvergence. Jako vstupní signál použijte harmonický signál A. Počáteční podmínky volte

1. nulové.

2. $w_1 = 2$,$w_2 = 1$ (nuly na jednotkové kružnici, fáze $\pi$).

3. $w_1 = 0$,$w_2 = 1$ (nuly na jednotkové kružnici, fáze $\frac{\pi}{2}$, $-\frac{\pi}{2}$).

4. $w_1 = -1$,$w_2 = 1$ (nuly na jednotkové kružnici, fáze $\frac{\pi}{3}$, $-\frac{\pi}{3}$).

Vyneste průběhy fáze a modulu nul filtru na čase. Porovnejte časy ustálení pro jednotlivé případy.

Figure 3.5: Průběhy vah a nul filtru pro počáteční podmínky 2. (nuly v $-1$).

Figure 3.6: Průběhy vah a nul filtru pro počáteční podmínky 4. (nuly v $e^{\jmath \frac {\pi }{3}}$, $e^{- \jmath \frac {\pi }{3}}$).

Cvičení 3.4: Jako vstupní signál použijte harmonický signál A. Přidejte k $x[n]$ bílý gausovský šum, rozptyl volte $0.3^2$. Jak přidání šumu ovlivňuje odhad?

Figure 3.7: Průběhy vah a nul filtru po přidání bílého šumu

Figure 3.8: Porovnání modulových frekvenčních charakteristik filtru v ustáleném stavu před a po přidání bílého šumu