Podstawy filtracji adaptacyjnej
Autorzy: Tomasz P. Zieliński
Pliki: rozdzial_07.zip
Przykład 7.1
W przykładzie 7.1 (plik filtrNLMS.m) są zaimplementowane niektóre odmiany filtrów adaptacyjnych LMS, opisanych powyżej. Wielkości wejściowe i wyjściowe oraz zmienne wewnętrzne mają takie same oznaczenia jak w tekście. Parametrem wejściowym ialg (numer algorytmu) wybiera się rodzaj filtra: 1 = LMS, 2 = NLMS, 3 = AP-NLMS, 4 = zdekorelowany LMS, 5 = zdekorelowany NLMS oraz 6 = LMS Newtona (LMS-N). Wszystkie algorytmy mogą być użyte w wersji „z wyciekiem”, jeśli wartość parametru leaky jest większa od zera. Zdekorelowane filtry ortogonalne będą omówione w rozdziale 8.function [y, e] = filtrNLMS(M,L,d,x,mi,gamma,alfa,delta,leaky,ialg)
% Przykład 7.1
% Autor: Tomasz Zieliński, tzielin@agh.edu.pl
bx = zeros(M,1); % inicjowanie bufora bx sygnału wejściowego x[n]
h = zeros(M,1); % inicjowanie wag filtra
coef = 1-mi*leaky; % współczynnik pamiętania wag h, zwykle 1, jeśli leaky=0
y = []; e = []; % wyzerowanie sygnałów wyjściowych
Rinv = delta*eye(M,M); % LMS-N: inicjowanie odwrotności macierzy autokorelacji sygnału x[n]
X = zeros(L,M); % AP-NLMS: macierz na próbki x[n]
bd = zeros(L,1); % AP-NLMS: inicjowanie bufora bd sygnału wejściowego d[n]
for n = 1 : length(x)
if(n==1) bxx = eps*ones(M,1); else bxx=bx; end % zapamiętaj poprzednią wartość
bx = [ x(n); bx(1:M-1) ]; % pobierz nową próbkę x[n] do bufora
if(ialg==3) % algorytm AP-NLMS
X = [ bx'; X(1:L-1,1:M) ]; % pobierz nowy wiersz do macierzy X
bd = [ d(n); bd(1:L-1) ]; % pobierz nową próbkę d[n] do bufora bd
yL = X*h; y(n) = yL(1); % oblicz yL, y[n]
eL = bd - yL; e(n)=eL(1); % oblicz eL, e[n]
else % algorytmy POZOSTAŁE
y(n) = h' * bx; % oblicz y[n] = sum( x .* bx)
e(n) = d(n) - y(n); % oblicz e[n]
end
if (ialg==1) % LMS
h = coef*h + mi * e(n) * bx;
end
if (ialg==2) % NLMS
energia = bx' * bx;
h = coef*h + mi/(gamma+energia) * e(n) * bx; % modyfikacja wag
end
if (ialg==3) % AP-NLMS
h = coef*h + mi* X'*inv( X*X'+gamma*eye(L) )*eL; % modyfikacja wag
end
if (ialg==4) % zdekorelowany LMS
alfa = (bx'*bxx) / (bxx'*bxx);
v = bx - alfa*bxx; %
h = coef*h + mi * e(n) * v; % modyfikacja wag
end
if (ialg==5) % zdekorelowany NLMS
alfa = (bx'*bxx) / (bxx'*bxx);
v = bx - alfa*bxx; %
h = coef*h + mi/(gamma+bx'*bx) * e(n) * v; % modyfikacja wag wer. 1
% h = coef*h + mi/(gamma+bx'*v) * e(n) * v; % modyfikacja wag wer.2
end
if (ialg==6) % filtr LMS-N (Newton); poniżej modyfikacja Rinv
Rinv = (Rinv - Rinv*bx*bx'*Rinv / ((1-alfa)/alfa+bx'*Rinv*bx)) / (1-alfa);
h = coef*h + mi*Rinv * bx * e(n); % modyfikacja wag
end
end % of for
Przykład 7.2
W przykładzie 7.2 (plik filtrWRLS1.m) zaimplementowano filtry WLS i WRLS. Wielkości wejściowe i wyjściowe oraz zmienne wewnętrzne są oznaczone tak, jak w powyższym wyprowadzeniu.function [y, e] = filtrWRLS1(M,d,x,lambda,delta)
% Przykład 7.2
% Autor: Tomasz Zieliński, tzielin@agh.edu.pl
bx = zeros(M,1); % inicjowanie bufora sygnału wejściowego x
h = zeros(M,1); % inicjowanie wag filtra
P = delta*eye(M,M); % inicjowanie odwrotności macierzy autokorelacji
% sygnału x[n]: P = R^{-1}
y = []; e = []; % wyzerowanie sygnałów wyjściowych
for n = 1 : length(x)
bx = [ x(n); bx(1:M-1) ]; % pobierz nową próbkę do bufora
y(n) = h' * bx; % oblicz y[n] = sum(h .* bx)
e(n) = d(n) - y(n); % oblicz e[n]
if(0) % wersja 1 - filtr RLS
P = (P + P*bx*bx'*P/(lambda+bx'*P*bx))/lambda; % zmodyfikuj macierz P
h = h + P * bx * e(n); % modyfikacja wag filtra
else % wersja 2 - obserwator RLS z [8], [22]
K = P*bx/(lambda+bx'*P*bx); % wzmocnienie Kalmana
h = h + K * e(n); % modyfikacja wag filtra
P = (eye(M)-K*bx')*P/lambda; % zmodyfikuj macierz P
end
end
Przykład 7.3
Ponieważ adaptacyjne filtry WRLS oraz adaptacyjna estymacja WRLS stanowią inny opis tego samego problemu, równania (7.68) (7.72) możemy wykorzystać do inicjowania początkowych wag filtra h[-1] oraz macierzy Rxx[-1] w algorytmie adaptacyjnego filtra WRLS. Szczegóły programowe przedstawiono poniżej (plik filtrWRLS2.m).function [y, e] = filtrWRLS2(M,d,x,lambda,delta,ialg)
% Przykład 7.3
% Autor: Tomasz Zirliński, tzielin@agh.edu.pl
% INICJOWANIE
if(ialg==0) % start od zera
bx = zeros(M,1); % inicjowanie bufora sygnału wejściowego x
h = zeros(M,1); % inicjowanie wag filtra
P = delta*eye(M,M); % inicjowanie odwrotności macierzy autokorelacji
% sygnału x[n]; P=R 1
n1 = 1; % indeks startowy pętli iteracyjnej
else % start od pierwszej estymaty obliczonej w sposób blokowy
w = lambda.^(M-1:-1:0); w=w'; % wektor wag
W = diag(w,0); % macierz wag
X = toeplitz( x(M:2*M-1), x(M:-1:1)); % macierz sygnału filtrowanego
d1 = d(M:2*M-1)'; % wektor sygnału referencyjnego
P = pinv(X'*W*X)*X'*W; % inicjowanie odwrotności macierzy autokorelacji
h = P*d1; % inicjowanie wag filtra
bx = x(2*M-1:-1:M).'; % inicjowanie bufora sygnału wejściowego x
n1 = 2*M; % indeks startowy pętli iteracyjnej
end
y = []; e = []; % wyzerowanie sygnałów wyjściowych
for n = n1 : length(x)
bx = [ x(n); bx(1:M-1) ]; % pobierz nową próbkę do bufora
y(n) = h' * bx; % oblicz y[n] = sum( x .* bx)
e(n) = d(n) - y(n); % oblicz e[n]
if(1) % wersja 1 - podobna do LMS Newtona
P = (P - P*bx*bx'*P/(lambda+bx'*P*bx))/lambda; % zmodyfikuj macierz P = inv(R)
h = h + P * bx * e(n); % modyfikacja wag
else % wersja 2 - ze wzmocnieniem Kalmana [8],[22]
K = P*bx/(lambda+bx'*P*bx); % wzmocnienie Kalmana
h = h + K * e(n); % modyfikacja wag
P = (eye(M)-K*bx')*P/lambda; % zmodyfikuj macierz P = inv(R)
end
end