Adaptacyjne filtry ortogonalne
Autorzy: Tomasz P. Zieliński
Pliki: rozdzial_08.zip
Przykład 8.1
W przykładzie 8.1 (plik przyklad_8_1.m)jest zaprezentowany program, w którym zaimplementowano przykładowy filtr FIR w postaci kratowej. Należy zwrócić uwagę, że błędy e+k [n] i e?k [n] oznaczono w nim jako eak[n] i ebk[n] (odpowiednio zmienne ea(k) i eb(k)), czyli znaki + i – zastąpiono literami a i b.% Przykład 8.1
% Autor: Tomasz Zieliński, tzielin@agh.edu.pl
clear all; clf;
% Podaj współczynniki transmitancji H(z)=A(z) filtra (UWAGA! a(1)=1, M<=N)
a = [ 1 -0.9 0.64 -0.576 ];
N = length(a);
aa = a(2:N); M=N-1; % odrzucamy pierwszy element równy 1
% Oblicz współczynniki "odbicia" gamma(i) filtra kratowego - tabela 8.1,
% równania (8.21)(8.22)
aij(M,1:M)=aa(1:M);
for i=M:-1:1
gamma(i)=-aij(i,i);
for j=1:i-1
aij(i-1,j) = (aij(i,j)+gamma(i)*aij(i,i-j))/(1-gamma(i)^2);
end
end
gamma, pause
% Odtwórz współczynniki {a} na podstawie wsp. {gamma} - tabela 8.1, (8.16),(8.17)
aa = zeros(M,M);
aa(1,1) = -gamma(1);
for i=2:M
aa(i,i) = -gamma(i);
for j=1:i-1
aa(i,j) = aa(i-1,j) - gamma(i)*aa(i-1,i-j);
end
end
aa = [ 1 aa(M,:) ], pause
% Sygnał filtrowany
Nx=100; n=0:Nx-1; x=sin(2*pi*n/12+pi/4);
% Filtracja kratowa "tylko zera" (rys. 8.1b i 8.5a):
% H1(z)=A(z)=[ 1 -0.9 0.64 -0.576 ];
y1 = filtrK_Z( x, gamma);
plot(y1); title('Sygnał wyjściowy - filtr: tylko zera'); pause
% Porównanie ze zwykła filtracją
y1z = filter(a,1,x);
plot(y1-y1z); title('Różnica z filter(a,1,x)');
%-----------------------------------------------------------------------------------
% function y = filtrK_Z(x,gamma)
% Filtr kratowy - tylko zera
% Autor: Tomasz Zieliński, tzielin@agh.edu.pl
% N=length(x);
% M=length(gamma);
% eb(1:M) = zeros(1,M);
% for n = 1 : N
% ea(1) = x(n); % kolejna próbka
% for k = 1 : M % kolejne sekcje filtra (8.2)
% ea(k+1) = ea(k) - gamma(k)*eb(k);
% ebnext(k+1) = eb(k) - gamma(k)*ea(k);
% end
% eb = [ x(n) ebnext(2:M) ];
% y(n) = ea(k+1);
% end
Przykład 8.2
W przykładzie 8.2 (plik przyklad_8_2.m) jest zaprezentowany program, w którym zaimplementowano w postaci kratowej przykładowy filtr IIR typu tylko-bieguny. Należy zwrócić uwagę, że błędy ek+[n] i ek-[n] oznaczono w nim jako eak[n] i ebk[n] (odpowiednio zmienne ea(k)i eb(k)), czyli + i – zastąpiono literami a i b.% Przykład 8.2
% Autor: Tomasz Zieliński, tzielin@agh.edu.pl
clear all; clf;
% Podaj współczynniki transmitancji H(z)=1/A(z) filtra (UWAGA! a(1)=1, M<=N)
a = [ 1 -0.9 0.64 -0.576 ];
% Podaj współczynniki transmitancji H(z)=A(z) filtra (UWAGA! a(1)=1, M<=N)
a = [ 1 -0.9 0.64 -0.576 ];
N = length(a);
aa = a(2:N); M=N-1; % odrzucamy pierwszy element równy 1
% Oblicz współczynniki "odbicia" gamma(i) filtra kratowego - tabela 8.1,
% równania (8.21)(8.22)
aij(M,1:M)=aa(1:M);
for i=M:-1:1
gamma(i)=-aij(i,i);
for j=1:i-1
aij(i-1,j) = (aij(i,j)+gamma(i)*aij(i,i-j))/(1-gamma(i)^2);
end
end
gamma, pause
% Sygnał filtrowany
Nx=100; n=0:Nx-1; x=sin(2*pi*n/12+pi/4);
% Filtracja kratowa "tylko bieguny" (rys. 8.2b i 8.5b):
% H2(z)=1/A(z)=1/[ 1 -0.9 0.64 -0.576 ];
y2 = filtrK_B( x, gamma);
plot(y2); title('Sygnał wyjściowy - filtr: tylko bieguny'); pause
% Porównanie ze zwykła filtracją
y2z = filter(1,a,x);
subplot(111); plot(y2-y2z); title('Różnica z filter(1,a,x)');
%-----------------------------------------------------------------------------------
% function y = filtrK_B(x,gamma)
% Filtr kratowy - tylko bieguny
% N=length(x);
% M=length(gamma);
% eb(1:M) = zeros(1,M);
% for n = 1 : N
% ea(M+1) = x(n); % kolejna próbka
% for k = M : -1 : 1 % kolejne sekcje filtra (8.27)
% ea(k) = ea(k+1) + gamma(k)*eb(k);
% eb(k+1) = eb(k) - gamma(k)*ea(k);
% end
% eb = [ ea(1) eb(2:M+1) ];
% y(n) = ea(1);
% end
Przykład 8.3
W przykładzie 8.3 (plik przyklad_8_3.m) zaimplementowano filtr kratowy IIR typu zera-bieguny.% Przykład 8.3
% Autor: Tomasz Zieliński, tzielin@agh.edu.pl
clear all; clf;
% Podaj współczynniki transmitancji H(z)=B(z)/A(z) filtra (UWAGA! a(1)=1, M<=N)
b = [ 1 3 3 1 ]; a = [ 1 -0.9 0.64 -0.576 ];
P=length(b); N=length(a);
aa=a(2:N); M=N-1; % odrzucamy pierwszy element równy 1
% Oblicz współczynniki "odbicia" gamma(i) filtra kratowego - tabela 8.1
% równania (8.21)(8.22)
aij(M,1:M)=aa(1:M);
for i=M:-1:1
gamma(i)=-aij(i,i);
for j=1:i-1
aij(i-1,j) = (aij(i,j)+gamma(i)*aij(i,i-j))/(1-gamma(i)^2);
end
end
gamma, pause
% Oblicz współczynniki c(i) filtra kratowego - równanie (8.40)
for k = P:-1:1
s=0; for(j=k+1:P) s = s + c(j)*aij(j-1,j-k); end
c(k) = b(k) - s;
end
c, pause
% Sygnał filtrowany
Nx=100; n=0:Nx-1; x=sin(2*pi*n/12+pi/4);
% Filtracja "kratowa"
% "zera i bieguny" (rys. 8.4e i 8.6b) : H3(z)= [ 1 3 3 1 ] / [1 -0.9 0.64 -0.576 ]
y3 = filtrK_ZB( x, gamma, c);
% Rysunki
plot(y3); title('Sygnał wyjściowy - filtr: zera i bieguny'); pause
% Porównanie ze zwykła filtracją
y3z = filter(b,a,x);
plot(y3-y3z); title('Różnica z filter(b,a,x)'); pause
%-----------------------------------------------------------------------------------
% function y = filtrK_ZB(x,gamma,c)
% Filtr kratowy - zera i bieguny
% N=length(x);
% M=length(gamma);
% eb(1:M) = zeros(1,M);
% for n = 1 : N
% ea(M+1) = x(n); % kolejna próbka
% for k = M : -1 : 1 % kolejne sekcje filtra (8.27)
% ea(k) = ea(k+1) + gamma(k)*eb(k);
% eb(k+1) = eb(k) - gamma(k)*ea(k);
% end
% eb = [ ea(1) eb(2:M+1) ];
% y(n) = sum( c.*eb(1:M+1) );
% end
Przykład 8.4
W przykładzie 8.4 (plik corrx.m) przedstawiono programy w języku MATLAB, realizujące wszystkie obliczenia z projektowania kratowego filtra Wienera. Napisano je na podstawie programów w językach Fortran i C, zamieszczonych w literaturze [11]. W kolejności są to funkcie:przykład 8.4a = corrxy.m, przykład 8.4b = lev.m, przykład 8.4c = firw.m, przykład 8.4d = filtrLWF.m.
function [r] = corrxy(x,y,M) % Przykład 8.4a % Obliczanie M-próbek funkcji korelacji wzajemnej sygnałów x[n] i y[n] % Oblicz rxx lub rdx, w zależności od sygnałów wejściowych N = length(x); for k = 0 : M r(k+1) = sum( x(1+k:N) .* y(1:N-k) )/N; endlev.m
function [L,E] = lev(M,rxx)
% Przykład 8.4b
% Obliczanie współczynników filtrów predykcji liniowej rzędu od 1 do M
% na podstawie wektora rxx autokorelacji sygnału x[n]
% metodą Durbina-Levinsona - inny zapis niż w tabl. 8.2, taki jak w [11]
% Autor implementacji w Matlabie: Tomasz Zieliński, tzielin@agh.edu.pl
L(1,1)=1; L(2,2)=1; L(2,1)=-rxx(2)/rxx(1);
E(1)=rxx(1); E(2)=E(1)*(1-L(2,1)^2);
for k = 3:M+1
delta = sum( rxx(2:k) .* L(k-1,1:k-1) ); % równanie (8.58b)
gamma = delta/E(k-1); % równanie (8.68)
L(k,1) = - gamma;
for i = 2:k-1
L(k,i) = L(k-1,i-1) - gamma*L(k-1,k-1-i+1); % równanie (8.69)
end
L(k,k)=1;
E(k)=E(k-1)*(1-gamma^2); % równanie (8.70)
end
firw.m
function [L,E,g,h] = firw(M,rdx,rxx)
% Przykład 8.4c
% Projektowanie współczynników filtra FIR Wienera na podstawie [11]
% Autor implementacji w Matlabie: Tomasz Zieliński, tzielin@agh.edu.pl
[L,E] = lev(M,rxx); % algorytm Levinsona
for k=1:M+1
g(k) = sum( L(k,1:k) .* rdx(1:k) )/E(k); % g=Dinv*L*rdx, równanie (8.83)
end
for k=1:M+1
h(k) = sum( (L(k:M+1,k).') .* g(k:M+1) ); % h=L'*g, równanie (8.83)
end
filtrLWF.m
function [y,e] = filtrLWF(M,g,gamma,d,x)
% Przykład 8.4d
% Nieadaptacyjny kratowy filtr Wienera (Lattice Wiener Filter), patrz rys. 8.7b
% Na podstawie [11]
% Autor implementacji w Matlabie: Tomasz Zieliński, tzielin@agh.edu.pl
w = 0.000001*randn(1,M+1); % wektor roboczy, stany wewnętrzne każdej sekcji
for n = 1 : length(x)
ea(1) = x(n); % wejście do pierwszej sekcji filtra a priori
eb(1) = x(n); % j.w.
y(n) = g(1)*eb(1); % estymata y[n] bazująca na sekcji 1
e(n) = d(n) - y(n); % estymata e[n] bazująca na sekcji 1
for k = 2:M+1 % kolejne k-te sekcje filtra kratowego (8.2), rys.8.7b
ea(k) = ea(k-1) - gamma(k)*w(k); % błąd ek+[n]
eb(k) = w(k) - gamma(k)*ea(k-1); % błąd ek-[n]
w(k) = eb(k-1); % stan wewnętrzny k-tej sekcji
y(n) = y(n) + g(k)*eb(k); % estymata a posteriori xest[n] rzędu "k", równ. (8.84)
e(n) = e(n) - g(k)*eb(k); % błąd a posteriori e[n] rzędu "k", równanie (8.86)
end
end
Przykład 8.5
W przykładzie 8.5 (plik filtrGLLPF.m) przedstawiono implementację adaptacji wag filtra kratowego typu FIR pracującego jako liniowy predyktor.function [y] = filtrGLLPF(M,x,lambda,beta)
% Przykład 8.5 na podstawie [11]
% Adaptacyjny gradientowy kratowy liniowy predyktor.
% (Gradient Lattice Linear Prediction Filter)
% Rys. 8.7b tylko z adaptacją wag gamma filtra
% Autor implementacji w Matlabie: Tomasz Zieliński, tzielin@agh.edu
w = 0.000001*randn(1,M+1); % wektor roboczy, stany wewnętrzne każdej sekcji
gamma = w; da = w;
for n = 1 : length(x) % dla kolejnych próbek sygnału x[n]
ea(1) = x(n); % wejście do pierwszej sekcji filtra
eb(1) = x(n); % wejście do pierwszej sekcji filtra
for k = 2:M+1 % kolejne sekcje filtra
ea(k) = ea(k-1) - gamma(k)*w(k); % równanie (8.2)
eb(k) = - gamma(k)*ea(k-1) + w(k); % czyli "motylek"
da(k) = lambda*da(k) + ea(k-1)*ea(k-1) + w(k)*w(k); % nowe d[k], (8.95)(8.96)
gamma(k) = gamma(k) + beta*(ea(k)*w(k)+eb(k)*ea(k-1))/da(k); % nowe gamma[k], (8.97)
w(k) = eb(k-1); % aktualizacja stanu wewnętrznego k-tej sekcji filtra kratowego
end
y(n) = ea(k); % sygnał wyjściowy y[n]
end
Przykład 8.6
W przykładzie 8.6 (plik filtrGLWF.m) przedstawiono implementację adaptacyjnego gradientowego kratowego filtra Wienera typu FIR–LMS.function [y,e] = filtrGLWF(M,d,x,lambda,beta)
% Przykład 8.6
% Adaptacyjny gradientowy filtr kratowy Wienera typu FIR-LMS
% (Gradient Lattice Wiener Filter)
% Rys. 8.7b tylko z adaptacją wag filtra (współczynnik gamma (k) i g(k))
% Na podstawie [11]
% Autor implementacji w Matlabie: Tomasz Zieliński, tzielin@agh.edu.pl
w = 0.000001*randn(1,M+1); % stany wewnętrzne każdej sekcji
g = w; gamma = w; da = w; db = w;
init = 0;
for n = 1 : length(x) % dla kolejnych próbek sygnału x[n]
ea(1) = x(n); % wejście do pierwszej sekcji filtra
eb(1) = x(n); % wejście do pierwszej sekcji filtra
y(n) = g(1)*eb(1); % estymata y[n] bazująca na sekcji 1
e(n) = d(n) - y(n); % estymata e[n] bazująca na sekcji 1
db(1) = lambda*db(1) + eb(1)*eb(1); % uaktualnienie db(1), równanie (8.116)
g(1) = g(1) + beta*e(n)*eb(1)/db(1); % uaktualnienie g(1), równanie (8.115)
for k = 2:M+1 % kolejne sekcje filtra
ea(k) = ea(k-1) - gamma(k)*w(k); % równanie (8.2), czyli "motylek"
eb(k) = - gamma(k)*ea(k-1) + w(k); %
da(k) = lambda*da(k) + ea(k-1)*ea(k-1) + w(k)*w(k); % nowe d[k],(8.95)(8.96)
gamma(k) = gamma(k) + beta*(ea(k)*w(k)+eb(k)*ea(k-1))/da(k); % nowe gamma[k],(8.97)
y(n) = y(n) + g(k)*eb(k); % estymata x[n] rzędu "k", równ. (8.84)
e(n) = e(n) - g(k)*eb(k); % bład e[n] estymaty rzędu "k", równanie (8.86)
db(k) = lambda*db(k)+eb(k)*eb(k); % uaktualnienie dk [k] (8.116)
if(k<=init)
g(k) = g(k) + beta*e(n)*eb(k)/db(k); % uaktualnienie g[k] (8.115)
end
w(k) = eb(k-1); % aktualizacja stanu wewnętrznego k-tej sekcji filtra kratowego)
end
init = init+1;
end
Przykład 8.7a
W przykładzie 8.7a (plik filtrFDAF_overlap_save.m) przedstawiono implementację w języku MATLAB filtra adaptacyjnego FDAF – overlap-save.function [y, e, Hest, H ] = filtrFDAF_overlap_save(M,d,x,nr,mi,alfa,delta,gamma)
% Przykład 8.7a
% Autor: Tomasz Zieliński, tzielin@agh.edu.pl
lambda = 1-alfa; % współczynniki wykorzystywane podczas
P = delta*ones(1,M); % nakładania ograniczenia na gradient funkcji kosztu
bx = zeros(1,2*M); % inicjowanie bufora sygnału wejściowego x
H = fft( [ones(1,M)/M zeros(1,M) ] ); % inicjowanie wag filtra
y = []; e = []; Hest = []; % inicjowanie
for k = 0 : length(x)/M-1
x1 = x(k*M+1:k*M+M); % pobierz kolejny blok próbek x[n]
bx = [ bx(M+1:2*M) x1 ]; % włóż je do bufora
X = fft(bx); % FFT
by = real(ifft(X.*H)); % oblicz y[n]
y = [ y by(M+1:2*M)]; % oddziel poprawny blok próbek y[n]
e1 = d( k*M+1:k*M+M ) - y( k*M+1:k*M+M ); % oblicz blok próbek błędu e[n]
e = [ e e1]; % zapamiętaj jako sygnał wyjściowy
ee = [ zeros(1,M) e1 ]; % dodaj zera
E = fft(ee); % FFT
XE = conj(X) .* E; % iloczyn widm
if(1) % korekta gradientu f. kosztu w czasie (z tym dokładniej, bez tego - szybciej)
delta = ifft(XE ); % przejście z gradientem do dziedziny czasu
if(1) delta = [ delta(1:M) zeros(1,M)]; end % OK
if(0) delta = [ delta(1:M)/((x1*x1')+gamma) zeros(1,M)]; end % test
if(0)
P = lambda*P + alfa*abs(X(1:M)).^2; % test
delta = [ delta(1:M)./(P) zeros(1,M)];
end
XE = fft(delta);
end % powrót z gradientem do dziedziny częstotliwości
H = H + 2*mi*XE; % adaptacja wag filtra w dziedzinie częstotliwości
Hest = [ Hest; H( nr ) ]; % zapamiętanie wybranych wag filtra
end
Przykład 8.7b
W przykładzie 8.7b (plik filtrFDAF_overlap_add.m) przedstawiono implementację w języku MATLAB filtra adaptacyjnego FDAF – overlap-add.function [y, e, Hest, H ] = filtrFDAF_overlap_add(M,d,x,nr,mi,alfa,delta,gamma)
% Przykład 8.7b
% Autor: Tomasz Zieliński, tzielin@agh.edu.pl
xx = zeros(1,2*M); % inicjowanie bufora sygnału wejściowego x[n]
H = zeros(1,2*M); % inicjowanie wag filtra
X2 = zeros(1,2*M); % inicjowanie wartości poprzedniego widma sygnału x[n]
J = (-1).^(0:2*M-1); % macierz kompensacji czasowej poprzedniego widma sygnału x[n]
lambda = 1-alfa; % alfa = 0.01; współczynniki wykorzystywane podczas
P = delta*ones(1,M); % delta = 0.25 nakładania ograniczenia na gradient f. kosztu
y = []; e = []; Hest = [];
for k = 0 : length(x)/M-1
x1 = x(k*M+1:k*M+M); % pobranie kolejnego bloku próbek x[n]
xx = [ x1 zeros(1,M) ]; % dodanie zer
X1 = fft(xx); % FFT
X = X1 + J.*X2; % suma widm FFT: obecnego (X1) i poprzedniego (X2)
X2 = X1; % zapamiętanie obecnego widma
yy = ifft(X.*H); % obliczenie y[n]
y = [ y real(yy(1:M))]; % wydzielenie poprawnego bloku próbek y[n] i ich zapamiętanie
e1 = d( k*M+1:k*M+M ) - y( k*M+1:k*M+M ); % obliczenie bloku próbek sygnału błędu e[n]
e = [ e real(e1)]; % zapamiętanie go
ee = [ e1 zeros(1,M) ]; % dodanie zer
E = fft(ee); % FFT sygnału błędu e[n]
delta = ifft( conj(X) .* E); % obliczenie gradientu funkcji kosztu, przejście do dz. czasu
if(0) % # ograniczenie gradientu
delta = [ delta(1:M)/((x1*x1')+gamma) zeros(1,M) ]; % # w dziedzinie czasu
else % #
P = lambda*P + alfa*abs(X(1:M)).^2; % #
delta = [ delta(1:M)./P zeros(1,M)]; % #
end % #
Delta = fft(delta); % powrót z gradientem funkcji kosztu do dz. częstotliwości
H = H + 2*mi*Delta; % adaptacja wag filtra w dziedzinie częstotliwości
Hest = [ Hest; H( nr ) ]; % zapamiętanie wybranych wag filtra
end
Przykład 8.7c
W przykładzie 8.7c (plik filtrFDAF_overlap_save_part.m) przedstawiono implementację w języku MATLAB filtra adaptacyjnego PBFDAF – overlap-save z podziałem odpowiedzi impulsowej na mniejsze fragmenty (Partitional Block FDAF).function [y,e,Hest,H] = filtrFDAF_overlap_save_part(Nh,M,d,x,nr,mi,alfa,delta,gamma)
% Przykład 8.7c
% Autor: Tomasz Zieliński, tzielin@agh.edu.pl
Lb = length(x)/M; % liczba sekcji (bloków) próbek wejściowych
Lh = Nh/M; % liczba sekcji (bloków) wag filtra
lambda = 1-alfa; % alfa = 0.01; % współczynniki wykorzystywane podczas
P = delta*ones(1,M); % delta = 0.25 % ogranicz. gradientu funkcji kosztu
for k=1:Lh
H(k,1:2*M) = fft( [ones(1,M)/M zeros(1,M)] ); % inicjowanie macierzy H wag filtr
end
y = []; e = []; Hest = []; % inicjowanie
bx = zeros(1,2*M); % inicjowanie bufora bx na próbki x[n]
X = zeros(Lh,2*M); % inicjowanie macierzy widm X
for k = 0 : Lb-1
x1 = x(k*M+1:k*M+M); % pobranie kolejnej sekcji próbek wejściowych x[n]
bx = [ bx(M+1:2*M) x1 ]; % zmodyfikowanie zawartości bufora tych próbek
X = [ fft(bx); X(1:Lh-1,:) ]; % FFT bufora, aktualizacja macierzy widm sekcji sygnału x[n]
XH = X.*H; % filtracja blokowa sygnału x[n] w dz. częstotliwości
by = real(ifft(sum(XH))); % złożenie wyników filtracji (sum), powrót do dz. czasu
y = [ y by(M+1:2*M) ]; % zapisanie sygnału wyjściowego y[n]
e1 = d( k*M+1 :k*M+M ) - y( k*M+1 : k*M+M ); % obliczenie próbek błędu e[n]
e = [ e real(e1)]; % zapamiętanie ich
ee = [ zeros(1,M) e1 ]; % utworzenie wektora próbek błędu e[n]
E = fft(ee); % przejście do dziedziny częstotliwości
for k=1:Lh
XE(k,:) = conj(X(k,:)) .* E; % obliczenie gradientu funkcji kosztu w dz. częstotliwości
end
if(1) % opcjonalne ograniczenie gradientu funkcji kosztu
XE = XE.'; % transpozycja bez sprzężenia
delta = ifft(XE ); % przejście z gradientem do dziedziny częstotliwości
if(1) delta = [ delta(1:M,:); zeros(M,Lh)]; end % OK
if(0) delta = [ delta(1:M,:)/((x1*x1')+gamma); zeros(M,Lh)]; end % test
if(0) P = lambda*P + alfa*abs(X(1:M)).^2; delta = [ delta(1:M,:)./(P); % test
zeros(M,Lh)]; end
XE = fft(delta); % powrót z gradientem do dziedziny czasu
XE = XE.'; % transpozycja bez sprzężenia
end
H = H + 2*mi*XE; % adaptacja wag filtra w dziedzinie częstotliwości
Hest = [ Hest; H( nr ) ]; % zapamiętanie wybranych wag filtra
end