Laboratorium Cyfrowego Przetwarzania Sygnałów

Ćwiczenie 6 – Filtracja optymalna i adaptacyjna

- 1 -

Wydział Elektryczny
Zakład Automatyki

LABORATORIUM CYFROWEGO PRZETWARZANIA SYGNAŁÓW

Ćwiczenie 6

Filtracja optymalna i adaptacyjna

1. Cel ćwiczenia

• Zapoznanie się z projektowaniem filtrów optymalnych metodą minimalizacji wskaźnika jakości i

algorytmami LMS i RLS adaptacji współczynników filtrów.

• Przeprowadzenie eksperymentów symulacyjnych z filtrami adaptacyjnymi w środowisku

Simulink.

2. Podstawy teoretyczne

Zagadnienie estymacji (odtworzenia) użytecznego sygnału s(n) na podstawie innego dostępnego

sygnału x(n), który jest zaszumioną wersją s(n), tzn. x(n) = s(n) + w(n) (w(n) – niepożądane zakłócenie)
jest jednym z centralnych problemów cyfrowego przetwarzania sygnałów. Trudność polega na tym, że
widmo zakłócenia pokrywa się zwykle z widmem sygnału użytecznego i filtracja zakłócenia jest
związana z usunięciem części sygnału użytecznego. Dlatego filtrów optymalnych nie projektuje się w
dziedzinie częstotliwości, ale z wykorzystaniem statystycznych właściwości sygnału podlegającego
filtracji. W ćwiczeniu będziemy zajmować się estymacją liniową przy addytywnych zakłóceniach.

2.1. Filtracja optymalna

Struktura filtra optymalnego jest przedstawiona na rys.1. Ma on dwa sygnały wejściowe: sygnał

filtrowany x(n) oraz sygnał odniesienia d(n), i dwa sygnały wyjściowe: wynik filtracji y(n), który jest
estymatą sygnału odniesienia, tj.

ˆ

( )

( )

y n

d n

=

, oraz błąd estymacji e(n) = d(n) - y(n). Zakładamy, że

sygnały wejściowe są stacjonarne (w szerszym sensie). Ograniczymy się do filtrów nierekursywnych (o
skończonej odpowiedzi impulsowej SOI (FIR)):

1

2

0

1

2

( )

M

M

H z

h

h z

h z

h z

−

−

−

=

+

+

+

+

…

,

(6.1)

ponieważ gwarantują one stabilność i dlatego są najczęściej stosowane jako filtry adaptacyjne. Wyjście
filtra można zapisać w formie:

0

( )

(

)

( )

M

T

k

k

y n

h x n k

n

=

=

−

=

∑

h x

,

(6.2)

gdzie:

0

1

2

[ , , , ,

]

T

M

h h h

h

=

h

…

,

( ) [ ( ), (

1), (

2), , (

)]

T

n

x n x n

x n

x n M

=

−

−

−

x

…

(wskaźnik

T

oznacza

transpozycję).

Zadaniem filtra jest takie przekształcanie sygnału wejściowego x(n), aby wynik filtracji był jak

najbardziej zbliżony do sygnału odniesienia przy założonym kryterium błędu. Najczęściej stosuje się
kryterium minimalizacji błędu średniokwadratowego (mean-square error MSE)

2

[ ( )]

J

E e n

=

, (E[

⋅] oznacza wartość oczekiwaną)

(6.3)

gdzie

0

( )

( )

( )

( )

(

)

M

k

k

e n

d n

y n

d n

h x n k

=

=

−

=

−

−

∑

.

(6.4)

Laboratorium Cyfrowego Przetwarzania Sygnałów

Ćwiczenie 6 – Filtracja optymalna i adaptacyjna

- 2 -

Rys. 1. Struktura filtra optymalnego. W przypadku stacjonarnym parametry filtra H(z) są wyznaczane

jeden raz przez minimalizację kryterium błędu związanego z e(n)

Błąd dopasowania jest funkcją współczynników filtra {h

k

}:

(

)

2

2

( )

( )

( )

( )

( )

[ ( )] 2

T

T

n

T

n

dx

xx

J

E d n

n

E d n

⎡

⎤

=

−

=

−

+

⎢

⎥

⎣

⎦

h

h x

h r

h R h

,

(6.5)

gdzie

( )

[ ( ) ( )]

n

T

xx

E

n

n

=

R

x

x

jest macierzą autokorelacji sygnału wejściowego

( )

( )

,

(

)

[ (

) (

)],

,

0,1, ,

n

n

xx

xx

i j

R

i

j

E x n i x n j

i j

M

⎡

⎤ =

−

=

−

−

=

⎣

⎦

R

…

,

(6.6)

a

( )

[ ( ) ( )]

n

dx

E d n

n

=

r

x

- wektorem korelacji wzajemnej pomiędzy d(n) i x(n-i):

( )

( )

( )

[ ( ) (

)],

0,1, ,

n

n

xx

dx

i

r

i

E d n x n i

i

M

⎡

⎤ =

=

−

=

⎣

⎦

r

…

(6.7)

Optymalne wartości współczynników filtra otrzymuje się w wyniku przyrównania do zera pochodnej:

( )

( )

( )

2

2

T

n

T

n

dx

xx

J

∂

= −

+

=

∂

h

h r

h R h 0

h

,

(6.8)

co prowadzi do układu liniowego tzw. równań normalnych:

( )

( )

n

n

xx

dx

⋅ =

R

h r

(6.9)

Rozwiązanie tych równań daje optymalne współczynniki filtra:

( )

1 ( )

( )

[

]

, det

0

opt

n

n

n

xx

dx

xx

−

=

≠

h

R

r

R

(6.10)

1

( )

( )

( )

( )

0

( )

( )

( )

( )

1

( )

( )

( )

( )

(0)

(1)

( )

(0)

(1)

(0)

(

1)

(1)

( )

(

1)

(0)

( )

n

n

n

n

opt

xx

xx

xx

dx

n

n

n

n

opt

xx

xx

xx

dx

n

n

n

n

opt

xx

xx

xx

dx

M

R

R

R

M

r

h

R

R

R

M

r

h

R

M

R

M

R

r

M

h

−

⎡

⎤

⎡

⎤

⎡

⎤

⎢

⎥

⎢

⎥

⎢

⎥

−

⎢

⎥

⎢

⎥

⎢

⎥ =

⋅

⎢

⎥

⎢

⎥

⎢

⎥

⎢

⎥

⎢

⎥

⎢

⎥

−

⎢

⎥ ⎢

⎥

⎢

⎥

⎣

⎦ ⎣

⎦

⎣

⎦

Warunek odwracalności macierzy

( )

det

0

n

xx

≠

R

jest nazywany warunkiem wystarczającego pobudzenia.

Filtr (6.10) jest nazywany filtrem Wienera. Dla filtracji optymalnej zachodzi zależność:

( )

( )

2

( )

2 [ ( ) (

)] 0,

0,1, ,

k

k

J

e n

E e n

E e n x n k

k

M

h

h

⎡

⎤

∂

∂

=

= −

−

=

=

⎢

⎥

∂

∂

⎣

⎦

h

…

, (6.11)

co oznacza, że filtr usuwa korelację błędu estymacji ze wszystkimi poprzednimi (i bieżącą) próbkami
sygnału wejściowego:

( )

( )

[ ( ) (

)] 0

n

ex

r

k

E e n x n k

=

−

=

.

(6.12)

Jest to tzw. właściwość ortogonalności. Jej geometryczna interpretacja jest pokazana na rys.2.

H(z)

x(n)

błąd

e(n)

_

+

y(n)

sygnał odniesienia

d(n) = s(n)

zakłócenie

w(n)

Laboratorium Cyfrowego Przetwarzania Sygnałów

Ćwiczenie 6 – Filtracja optymalna i adaptacyjna

- 3 -

Rys. 2. Geometryczna interpretacja właściwości ortogonalności estymaty otrzymanej metodą błędu

średniokwadratowego MSE dla dwóch wymiarów (M=2). Wyjście filtra

ˆ

( )

( )

y n

d n

=

.

2.2. Estymacja współczynników filtra optymalnego z próbek

W zagadnieniach praktycznych teoretyczne kryterium (6.3) oraz korelacje (6.6) i (6.7) wyznaczane z

wartości oczekiwanych zastępuje się ich estymatami

ˆ

xx

R

i

ˆ

dx

r

obliczonymi z próbek sygnałów x(n) i d(n).

Załóżmy, że mamy do dyspozycji po N próbek sygnału wejściowego i sygnału odniesienia:

[ (0), (1), (2), , (

1)]

T

x

x

x

x N

=

−

x

…

,

[ (0), (1), (2), , (

1)]

T

d

d

d

d N

=

−

d

…

,

a próbki na wyjściu filtra określa równanie (6.2). Ponieważ zwykle liczba próbek jest większa niż długość
filtra N>M, otrzymuje się układ nadokreślony, tzn. liczba równań jest większa niż liczba niewiadomych
(poszukiwanych współczynników filtra). Istnieją dwie metody obliczeń w zależności od zakresu
sumowania we wzorze na estymatę błędu dopasowania. Każda z nich daje inne wartości estymat

ˆh

współczynników filtra optymalnego.
1. Metoda autokorelacji

Kryterium błędu:

1

2

0

( ),

( )

( )

( )

N M

n

n

J

e n

e n

d n

y n

+ −

=

=

=

−

∑

(6.13)

W metodzie tej mamy N+M równań na wartości wyjściowe y(n). Kolejne próbki y(n) są obliczane jako
suma iloczynów współczynników filtra i odpowiadających im próbek x(n), jak na poniższym rysunku. Jak
widać, ciąg próbek sygnału wejściowego jest na krańcach uzupełniany zerami. Poza pierwszymi N
równaniami w pozostałych po prawej stronie występują zera.

Metoda autokorelacji prowadzi do symetrycznej macierzy równań normalnych i gwarantuje, że filtr H(z)
jest minimalnofazowy.

x

0

x

1

x

2

x

n-M

x

n-M+1

x

n-1

x

n

x

N-1

0

0

0

0 0

h

M

h

1

h

0

.....

...... .....

...

....

h

M

h

1

h

0

.....

h

M

h

1

h

0

.....

....

d

0

d

1

d

2

d

n-M

d

n-M+1

d

n-1

d

n

d

N-1

0

0 0

..... ...

....

....

Laboratorium Cyfrowego Przetwarzania Sygnałów

Ćwiczenie 6 – Filtracja optymalna i adaptacyjna

- 4 -

2. Metoda kowariancji

Kryterium błędu:

1

2

( ),

( )

( )

( )

N

n

n M

J

e n

e n

d n

y n

−

=

=

=

−

∑

(6.14)

W metodzie tej mamy N-M równań na wartości wyjściowe y(n). Kolejne próbki y(n) są obliczane jako
suma iloczynów współczynników filtra i odpowiadających im próbek x(n), jak na poniższym rysunku.
Tworzenie równań jest ograniczone do takiego zakresu, że współczynniki filtra "nie wychodzą" poza
rekord próbek wejściowych. Zwróćmy uwagę, że nie wykorzystuje się informacji z pierwszych M próbek
sygnału odniesienia.

Metoda kowariancji nie gwarantuje minimalnofazowości H(z).

2.3. Filtracja adaptacyjna

Filtracja adaptacyjna polega na dodaniu do struktury filtra optymalnego algorytmu adaptacji

współczynników H(z) (rys. 3), które są funkcjami czasu n i są na bieżąco dopasowywane w taki sposób,
aby zapewniać optymalną, w sensie przyjętego kryterium błędu, filtrację sygnału wejściowego.

Rys. 3. Struktura filtra adaptacyjnego. Parametry h(n) filtra H

n

(z) są wyznaczane na bieżąco według

określonego algorytmu adaptacyjnego minimalizującego kryterium błędu związane z e(n)

Jeżeli wszystkie sygnały w układzie są stacjonarne, to po okresie adaptacji współczynniki filtra powinny
się ustalić i przyjąć wartości optymalne (warunek zbieżności procesu adaptacji):

lim ( )

opt

n

n

→∞

=

h

h

(6.15)

W adaptacyjnym przetwarzaniu sygnałów stosuje się metody optymalizacji wieloparametrycznej

(optymalizujemy M+1 parametrów transmitancji), które polegają poszukiwaniu minimów
deterministycznych wieloargumentowych funkcji kosztu J(

⋅). Najczęściej stosowane są metody

gradientowe, w których modyfikacja

Δh(n) jest w każdej chwili proporcjonalna do wektora ujemnego

gradientu funkcji kosztu:

( )

1

( )

(

1)

( )

(

1)

2

n

n

n

n

n

n

=

− + Δ

=

− − μ

⋅∇

h

h

h

h

W

(6.16)

0

1

( ( ))

( )

( )

( )

,

,

,

,

( )

( )

( )

( )

T

n

M

J

n

J n

J n

J n

n

h n

h n

h n

⎡

⎤

∂

∂

∂

∂

∇ =

= ⎢

⎥

∂

∂

∂

∂

⎣

⎦

h

h

…

(6.17)

H

n

(z)

x(n)

błąd

e(n)

_

+

y(n)

sygnał odniesienia

d(n)

Algorytm

adaptacji

x

0

x

1

x

n-M

x

n-M+1

x

n-1

x

n

x

N-M-1

h

M

h

1

h

0

.....

.... ......

...

....

h

M

h

1

h

0

.....

h

M

h

1

h

0

.....

...

x

M-1

x

M

x

M+1

x

N-M

x

N-2

x

N-1

d

0

d

1

d

n-M

d

n-M+1

d

n-1

d

n

d

N-M-1

.... ......

...

....

...

d

M-1

d

M

d

M+1

d

N-M

d

N-2

d

N-1

Laboratorium Cyfrowego Przetwarzania Sygnałów

Ćwiczenie 6 – Filtracja optymalna i adaptacyjna

- 5 -

W

(n)

jest dodatkową macierzą wagową poprawiającą zbieżność (np. odwrotność hesjanu

( )

2

1

[

]

n

n

−

= ∇

W

w procedurze Newtona-Raphsona). Współczynnik adaptacji

μ decyduje o szybkości przestrajania (Δh(n)

jest proporcjonalne do

μ). Ze względu na zbieżność algorytmy zachodzi warunek ograniczający:

max

2

0

< μ <

λ

,

(6.18)

gdzie

λ

max

jest największa wartością własną macierzy autokorelacji sygnału wejściowego

( )

n

xx

R

.

2.4. Filtry adaptacyjne LMS

W algorytmie adaptacyjnym LMS (Least Mean Squares) kryterium błędu

2

( )

n

J

e n

=

.

(6.19)

Zadaniem jest więc minimalizacja chwilowej (a nie oczekiwanej) wartości błędu kwadratowego, dlatego
filtrację nazywa się filtracją bez pamięci albo optymalizacją stochastyczną. Jednak asymptotycznie dla
n

→∞ algorytm ten minimalizuje średni błąd kwadratowy.

Metoda ta obejmuje szeroką rodzinę algorytmów opisanych ogólną zależnością:

( )

(

1)

( ) ( ) ( ) ( )

n

n

n

n e n

n

=

− + μ

h

h

W

x

(6.20)

Przypadki szczególne:
1. Filtr LMS

( )

(

1)

( ) ( )

n

n

e n

n

=

− + μ

h

h

x

(6.21)

( )

( )

( )

e n

d n

y n

=

−

( )

(

1) ( )

T

y n

n

n

=

−

h

x

Zaletą algorytmu jest prostota i mała złożoność obliczeniowa, a wadą – wolna zbieżność algorytmu.
2. Unormowany filtr LMS (NLMS)

( )

(

1)

( ) ( ) ( )

n

n

n e n

n

=

− + μ

h

h

x

,

(6.22)

gdzie

2

0

( )

( ) ( )

(

)

M

T

k

n

n

n

x n k

=

μ

μ

μ

=

=

γ +

γ +

−

∑

x

x

(mały parametr

γ zapobiega zerowaniu się mianownika).

Filtry NLMS charakteryzują się szybszą zbieżnością i lepszą stabilnością w porównaniu z LMS.
3. Zdekorelowany filtr LMS
Zamiast x(n) na wejście podawany jest zdekorelowany sygnał

( ) (

1)

( )

( )

( ) (

1),

( )

(

1) (

1)

T

T

n

n

n

n

n

n

n

n

n

−

=

− α

−

α

=

−

−

x

x

v

x

x

x

x

(6.23)

który ma stosunek

λ

min

/

λ

max

bliski jedności, co zwiększa szybkość adaptacji.

α(n) jest współczynnikiem

korelacji wektorów x(n) i x(n-1).

2.5. Filtry adaptacyjne RLS

Popularnym algorytmem adaptacji jest algorytm LS (Least Squares) z kryterium błędu w postaci sumy

kwadratów:

2

0

( )

n

n

k

J

e n

=

=

∑

(6.24)

lub ważonej sumy kwadratów WLS (Weighted LS):

2

0

( )

n

n k

n

k

J

e n

−

=

=

λ

∑

,

(6.25)

gdzie

λ jest współczynnikiem zapominania starych błędów o wartości z zakresu

0.95

1

≤ λ ≤

. Algorytm

z wagami umożliwia adaptację do zmieniających się charakterystyk sygnałów, kiedy dane z przeszłości
przestają być aktualne. Ponieważ kryterium zależy od sumy błędów, filtry LS nazywa się filtrami z

Laboratorium Cyfrowego Przetwarzania Sygnałów

Ćwiczenie 6 – Filtracja optymalna i adaptacyjna

- 6 -

pamięcią. W trakcie adaptacji ruch trajektorii układu w przestrzeni parametrów w kierunku minimum jest
mniej chaotyczny niż w przypadku metody LMS.

Równanie zmian współczynników można zapisać w postaci:

( )

(

1)

( ) ( )

n

n

n e n

=

− +

h

h

K

,

(6.26)

gdzie

( )

1

ˆ

( ) [

]

( )

n

xx

n

n

−

=

K

R

x

nazywa się wzmocnieniem Kalmana.

Ze względu na konieczność odwracania w każdym kroku macierzy autokorelacji algorytm LS

implementuje się w formie rekurencyjnej RLS (Recursive LS)

( )

(

1)

( ) ( )

n

n

n e n

=

− +

h

h

K

( )

( )

( )

e n

d n

y n

=

−

,

( )

(

1) ( )

T

y n

n

n

=

−

h

x

1

1

(

1) ( )

( )

( ) ( )

1

( ) (

1) ( )

T

n

n

n

n

n

n

n

n

−

−

λ

−

=

=

+ λ

−

P

x

K

P

x

x

P

x

, (gdzie

( )

1

ˆ

( ) [

]

n

xx

n

−

=

P

R

)

(6.27)

1

1

( )

[

( ) ( )] (

1)

T

n

n

n

n

−

−

= λ

− λ

−

P

I

K

x

P

, (6.28)

z warunkami początkowymi:

0

2

1

(0)

,

(0)

ˆ

x

=

=

σ

P

I

h

h

.

Algorytm RLS jest bardziej złożony obliczeniowo niż LMS, ale zwykle daje szybszą i płynniejszą
zbieżność.

2.6. Typowe zastosowania filtracji adaptacyjnej

Poniżej przedstawiono kilka typowych struktur zastosowania filtracji adaptacyjnej ilustrujących

interpretację poszczególnych sygnałów.

We wszystkich przypadkach wyjście filtra y(n) jest estymatą tylko tej składowej sygnału odniesienia

d(n), która jest skorelowana z sygnałem wejściowym x(n). Jeżeli d(n) = s(n) + x

1

(n), gdzie tylko x

1

(n) jest

skorelowane z x(n), to wyjście

1

ˆ

( )

( )

y n

x n

=

.

Rys. 4. Identyfikacja adaptacyjna. Transmitancja filtra adaptacyjnego jest estymatą transmitancji

identyfikowanego układu

ˆ

( )

( )

H z

G z

=

.

Rys. 5. Adaptacyjna predykcja sygnału

H

n

(z)

x(n)=d(n-k)

błąd predykcji

e(n)

_

+

ˆ

( )

( )

y n

d n

=

sygnał odniesienia

d(n) = s(n)

Predyktor k-krokowy

z

-k

ˆ

( )

( )

y n

d n

=

H

n

(z)

x(n)

e(n)

_

+

Filtr adaptacyjny

- model układu

G(z)

Nieznany układ

zakłócenie w(n)

d(n)

Laboratorium Cyfrowego Przetwarzania Sygnałów

Ćwiczenie 6 – Filtracja optymalna i adaptacyjna

- 7 -

Rys. 6. Adaptacyjne kasowanie szumu (usuwanie korelacji)

Rys. 6. Przykład zastosowania adaptacyjnego kasowania szumu. Zakłócenie skorelowane z szumem

własnym okrętu dociera do sonaru jako jego stłumiona i opóźniona wersja w(n)=

α⋅w

0

(n-n

0

). Zadaniem

filtracji adaptacyjnej jest wytłumienie tej składowej w sygnale sonaru.

Literatura

1. Stranneby D.: „Cyfrowe przetwarzanie sygnałów. Metody, algorytmy, zastosowania”, Wyd. BTC, 2004.
2. Zieliński T.P.: „Cyfrowe przetwarzanie sygnałów. Od teorii do zastosowań”, WKŁ, 2005.

'

#

szum własny w

0

(n)

sygnał sonaru d(n)=s(n)+w(n)

echo obiektu s(n)

H

n

(z)

x(n)=w

0

(n)

ˆ

( )

( )

e n

s n

=

_

+

d(n)=s(n)+w(n)

Filtr adaptacyjny

Źródło

szumu

Źródło

sygnału

ˆ

( )

( )

y n

w n

=

Laboratorium Cyfrowego Przetwarzania Sygnałów

Ćwiczenie 6 – Filtracja optymalna i adaptacyjna

- 8 -

3. Obliczenia komputerowe - zadania do wykonania

Do obliczeń wykorzystywana jest biblioteka symulacyjna Simulinka DSP Blockset oraz biblioteka
Matlaba Signal Processing Toolbox.

3.1. Liniowa predykcja adaptacyjna

A. Otworzyć model symulacyjny .mdl do predykcji liniowej:
>> c6_predict

Liniowa predykcja adaptacyjna

Osc

Delay

Flip

ZOH

User

Wsp. predykcji h(n)

t

Wektor chwil czasu

Szum pomiarowy

Signal

Generator

x

S.wejsciowy

d

S. odniesienia

dp

Predykcja

In

Err

Out

Taps

nLMS

LMS

Adaptive Filter

z

-2

Clock

d(n)

y(n) = dp(n)

e(n)

x(n)

• Ustawić (sprawdzić) parametry symulacji w menu Simulation | Parameters:

Stop time=1000, Solver options: Fixed-step, discrete

Ustawić (sprawdzić) parametry bloków:
Signal generator: Waveform: sine, Amplitude=1, Frequency=0.015 Hz,
Sample time=1 (okres próbkowania)
Szum: Variance=0.0 (brak zakłócenia), Sample time=1
Delay=2 (opóźnienie sygnału wejściowego filtra w próbkach, predykcja 2-krokowa)
LMS Filter: FIR length=16, Step-size mu=1.0, Use normalization: Tak

(długość filtra L, współczynnik adaptacji

μ, algorytm NLMS z normalizacją)

• Otworzyć okno oscyloskopu OSC i przeprowadzić symulację (przycisk Play  lub Simulation |

Start). Obserwować przebieg adaptacji.

• Zapisać wykres z końcową charakterystyką impulsową filtra predykcyjnego LMS h(n) (Copy Figure).

W oknie komend uruchomić program skryptowy do wykreślania sygnałów zarejestrowanych w
pamięci:

» c6_predictplot(t,d,x,dp)
Wkreślane są: sygnał odniesienia d(n), wejście filtra x(n) (opoźniona i zaszumiona wersja sygnału
odniesienia), wyjście filtra

ˆ

( )

( )

y n

d n

=

(predykcja) oraz błąd predykcji

ˆ

( )

( )

( )

e n

d n

d n

=

−

. Zapisać okno

z wykresami.
B. Powtórzyć obliczenia i zarejestrować wyniki dla pozostałych 3 kombinacji parametrów: wariancja
szumu zakłócenia Variance=0.0, 0.05, opóźnienie Delay=2, 10.
C. Porównać przebiegi adaptacji dla jednej z kombinacji parametrów z szumem dla innej wartości
współczynnika adaptacji, np. mu=0.3 (filtr LMS). Sprawdzić, jak działa adaptacja LMS bez
normalizacji (usunąć zaznaczenie opcji Use normalization w filtrze).
&

Przeanalizować uzyskane wyniki i wyciągnąć wnioski.

3.2. Estymacja opóźnienia

A. Otworzyć model symulacyjny .mdl do predykcji liniowej:

Laboratorium Cyfrowego Przetwarzania Sygnałów

Ćwiczenie 6 – Filtracja optymalna i adaptacyjna

- 9 -

>> c6_delay

Estymacja opoznienia

Flip

y

Wyjscie y(n)

User

Wsp.filtra

t

Wektor chwil czasu

x

Wejscie x(n)

Szum pomiarowy

Sine Wave

d

S.odniesienia d(n)

Results

In

Err

Out

Taps

nLMS

LMS

Adaptive Filter

1

Kx

0.5

z

-4

Delay

Clock

Echo d(n)

Error e(n)

Input x(n)

y(n)

• Ustawić (sprawdzić) parametry symulacji w menu Simulation | Parameters:

Stop time=0.5, Solver options: Fixed-step, discrete

Ustawić (sprawdzić) parametry bloków:
Sine wave: Amplitude=1, Frequency=50 Hz, Sample time=1e-3 (okres

próbkowania T

s

=1 ms)

Szum: Variance=0.2, Sample time=1e-3
Delay=4 (opóźnienie echa – składowej sygnału odniesienia)

LMS Filter: FIR length=16, Step-size mu=1.0, Use normalization: Tak

(długość filtra L, współczynnik adaptacji

μ, algorytm NLMS z normalizacją)

• Otworzyć okno oscyloskopu OSC i przeprowadzić symulację (przycisk Play  lub Simulation |

Start). Obserwować przebieg adaptacji.

• Zapisać wykres z końcową charakterystyką impulsową filtra predykcyjnego LMS h(n) - estymatę

opóźnienia (Copy Figure). W oknie komend uruchomić program skryptowy do wykreślania sygnałów
zarejestrowanych w pamięci:

» c6_delayplot(t,d,x,y)
Wkreślane są: sygnał odniesienia d(n) (wejście filtra + echo), wejście filtra x(n) (zaszumiony sygnał
harmoniczny), wyjście filtra y(n), błąd estymacji

( )

( )

( )

e n

d n

y n

=

−

. Zapisać okno z wykresami.

B. Powtórzyć obliczenia i zarejestrować wyniki dla czasu opóźnienia Delay=10.
C. Porównać przebieg adaptacji dla innej wartości współczynnika adaptacji, np. mu=0.3 (filtr LMS).
D. Powtórzyć obliczenia i zarejestrować wyniki dla samego echa jako sygnału odniesienia : parametr
Kx=0 (przerwanie gałęzi przepływu sygnału).
&

Zwrócić uwagę na estymatę opoźnienia h(n) (parametry filtra). Wyjaśnić wynik.

E. Przywrócić Kx=1. Powtórzyć obliczenia i zarejestrować wyniki przy braku szumu zakłócenia

Szum: Variance=0.
&

Zwrócić uwagę na estymatę opoźnienia h(n) (parametry filtra). Wyjaśnić wynik.

3.3. Identyfikacja układu

A. Otworzyć model symulacyjny .mdl do predykcji liniowej:
>> c6_ident

Laboratorium Cyfrowego Przetwarzania Sygnałów

Ćwiczenie 6 – Filtracja optymalna i adaptacyjna

- 10 -

Identyfikacja adaptacyjna

x

wymuszenie + szum

ZOH

User

Wsp. estymatora he(n)

t

Wektor chwil czasu

fir1

Uklad identyfikowany (FIR)

Szum pomiarowy

Signal

Generator

d

S.wyjsciowy

In

Err

Out

Taps

nLMS

LMS

Adaptive Filter

FFT

Freq Response

y

Estymacja y

e

Error e(n)

Clock

y(n) = de(n)

d(n)

x(n)

e(n)

• Ustawić (sprawdzić) parametry symulacji w menu Simulation | Parameters:

Stop time=1500, Solver options: Fixed-step, discrete

Ustawić (sprawdzić) parametry bloków:
Signal generator: square (wymuszenie prostokątne), Amplitude=1,

Frequency=0.02 Hz, Sample time=1
Szum: Variance=0.0 (brak zakłócenia), Sample time=1
Układ identyfikowany (FIR): Bandstop, Filter order=32 (filtr

pasmowozaporowy), Lower cutoff freq=0.3, Upper cutoff freq=0.7, Hamming
window (częstotliwość Nyquista = 1).

LMS Filter: FIR length=32, Step-size mu=1.0, Use normalization: Tak

(długość filtra L, współczynnik adaptacji

μ, algorytm NLMS z normalizacją)

• Otworzyć okno oscyloskopu OSC i przeprowadzić symulację (przycisk Play  lub Simulation |

Start). Obserwować przebieg adaptacji.

• Zapisać wykres z końcową charakterystyką impulsową filtra LMS h(n) - model dynamiki

identyfikowanego układu oraz jego charakterystykę amplitudową (Copy Figure). (Uwaga: w oknie
charakterystyki filtra LMS częstotliwość Nyquista = 0.5 – połowa częstotliwości próbkowania). W
oknie komend uruchomić program skryptowy do wykreślania sygnałów zarejestrowanych w pamięci:

» c6_identplot(t,d,x,y)
B. Powtórzyć obliczenia i zarejestrować wyniki z szumem zakłócenia: Szum: Variance=0.2.
C. Zmienić identyfikowany układ na pasmowoprzepustowy – wybrać z listy Filter type:
Bandpass. Przeprowadzić symulację (z szumem) i zarejestrować wyniki.
&

Czy charakterystyka częstotliwościowa filtra adaptacyjnego we właściwy sposób odwzorowuje
charakterystykę identyfikowanego układu?

D. Ustawić długi czas symulacji Simulation | Parameters: Stop time=50000.
Uruchomić symulację i obserwować przebieg adaptacji przy zmienianiu parametrów on-line. Można
zmieniać typ (np. BP

→BS→LP) i częstotliwości graniczne układu, wariancję szumu, współczynnik

adaptacji mu (np. mu=0.3) filtra adaptacyjnego, sygnał wymuszenia (wypróbować np. wymuszenie
sinusoidalne sine, zamiast prostokątnego square, przy braku szumu, tj. Variance=0.0).
&

Przeanalizować uzyskane wyniki i wyciągnąć wnioski.

E. Zamienić filtr adaptacyjny w schemacie blokowym na filtr RLS: zaznaczyć blok filtra LMS i usunąć
go, otworzyć bibliotekę bloków Simulinka Library Browser (przycisk na belce okna schematu) i wybrać

Laboratorium Cyfrowego Przetwarzania Sygnałów

Ćwiczenie 6 – Filtracja optymalna i adaptacyjna

- 11 -

DSP Blockset | Filtering | Adaptive filters | RLS Adaptive filter i przeciągnąć go w miejsce usuniętego
bloku (Uwaga: strzałki wejść i wyjść muszą się połączyć z portami bloku!). Parametry filtra RLS:
RLS Filter: FIR length=32, Memory weighting factor=1.0 (współczynnik
zapominania

λ), Initial value of taps=0.0, Initial variance estimate=0.1

Przywrócić wartości pozostałych parametrów jak w punkcie A. Powtórzyć obliczenia jak w punktach A,

B, C i zarejestrować wyniki.
&

Porównać wyniki z uzyskanymi dla filtra NLMS.

F. (*) Powtórzyć symulację z filtrem RLS z długim czasem symulacji jak w punkcie D. Zaobserwować
efekty po zmianie Memory weighting factor filtra RLS np. na 0.99.
&

Przeanalizować uzyskane wyniki i wyciągnąć wnioski.

3.4. Adaptacyjne kasowanie szumu

A. Otworzyć model symulacyjny .mdl do kasowania szumu:
>> c6_noisecanc

Adaptacyjne kasowanie szumu - filtr LMS/RLS

y(n)

h(n)

ZOH

User

Wsp. filtra

t

Wektor chwil czasu

Szum pierwotny

Switch

Signal Gen.

d

S.odniesienia (zaszumiony)

s

S. oryginalny

y

Output y(n)

Osc

In

Err

Out

Taps

nLMS

LMS

Adaptive Filter

[ts,s]

From

Workspace

FFT

Freq Response

Flip

fir1

Filtr LP szumu

e

Estymata se(n)

Clock

s(n)

d(n) = s(n) + w(n)

e(n) = se(n)

x(n) = w(n)

• Ustawić (sprawdzić) parametry symulacji w menu Simulation | Parameters:

Stop time=500, Solver options: Fixed-step, discrete

Ustawić (sprawdzić) parametry bloków:
Signal generator: sine (wymuszenie harmoniczne),

Amplitude=1,

Frequency=0.2 Hz, Sample time=1
Switch: przełącznik w pozycji Signal Gen. (dwukrotne kliknięcie zmienia położenie styku)
Szum: Variance=1.0, Sample time=1 (Uwaga: wariancja równa 0 oznacza zerowy

sygnał wejściowy filtra adaptacyjnego)

Filtr LP szumu: Lowpass, Filter order=32, Cutoff freq=0.5, Hamming
window (częstotliwość Nyquista = 1).

LMS Filter: FIR length=32, Step-size mu=0.3, Use normalization: Tak

(długość filtra L, współczynnik adaptacji

μ, algorytm NLMS z normalizacją)

W oknie komend zainicjalizować zmienne, np.:
>> s=0; ts=0;

• Otworzyć okno oscyloskopu OSC i przeprowadzić symulację (przycisk Play  lub Simulation |

Start). Obserwować przebieg adaptacji.

• Zapisać wykres z końcową charakterystyką impulsową filtra LMS h(n) - model dynamiki

identyfikowanego układu oraz jego charakterystykę amplitudową (Copy Figure). (Uwaga: w oknie

Laboratorium Cyfrowego Przetwarzania Sygnałów

Ćwiczenie 6 – Filtracja optymalna i adaptacyjna

- 12 -

charakterystyki filtra LMS częstotliwość Nyquista = 0.5 – połowa częstotliwości próbkowania). W
oknie komend uruchomić program skryptowy do wykreślania sygnałów zarejestrowanych w pamięci:

» c6_noisecancplot(t,d,s,y)
B. Rozszerzyć pasmo przepustowe filtra LP szumu zadając Cutoff freq=0.8 i powtórzyć
symulację.
&

Jak charakterystyka amplitudowa filtra adaptacyjnego ma się do widma szumu w(n) (pamiętając, że
widmo mocy P

ww

(

Ω) jest proporcjonalne do |H

LP

(e

j

Ω

)|

2

filtra szumu i że wyjście filtra ma

kompensować składową od zakłócenia w(n) w sygnale odniesienia d(n))?

C. Powtórzyć obliczenia i zarejestrować wyniki z szumem wejściowym o większej wariancji
Variance=10.0.
D. Przeprowadzić obliczenia i zarejestrować wyniki dla współczynnika adaptacji filtra LMS mu=1.0.
E. Przeprowadzić symulację i zaobserwować efekty dla filtracji LMS bez normalizacji (skasować
zaznaczenie Use normalization w bloku filtra LMS).
&

Porównać otrzymane wyniki i wyciągnąć wnioski. Jaki jest wpływ normalizacji w algorytmie LMS?

F. Przeprowadzić eksperyment z kasowaniem szumu w sygnale dźwiękowym.

• Wczytać sygnał dźwiękowy z pliku .wav o podanej nazwie, np. (*** potrzebne głosniki ***)

>> [s,ts,fs]=c6_sound('fire8'); % sygnał, czas, częstotliwość próbkowania
(sygnał jest wysyłany na urządzenie wyjściowe i rysowany jest przebieg czasowy).
W oknie komend podawana jest liczba próbek sygnału, dla fire8.wav jest to n=6588. W zależności
od długości sygnału ustawić czas symulacji w Simulation | Parameters, np. Stop time=6500.

• Przełączyć Switch w położenie From Workspace (dwukrotne kliknięcie w obszarze bloku).

Ustawić wariancję szumu Variance=3.0. Przeprowadzić symulację z wektorem sygnału x.

• Odsłuchać wyniki (*** potrzebne głosniki ***):

>> sound(d,fs); % sygnał zaszumiony
>> sound(e,fs); % błąd estymacji, estymata

ˆ( )

s n

sygnału oryginalnego po usunięciu szumu

• Zrobic i zarejestrować wykres estymaty

ˆ( )

s n

oraz różnicy

ˆ

( )

( )

s n

s n

−

:

» c6_estymplot(t,s,e);
G. Zamienić filtr adaptacyjny w schemacie blokowym na filtr RLS jak w punkcie 3.3.E (takie same
parametry filtra RLS). Powtórzyć obliczenia jak w punkcie D.

• Zmienić współczynnik zapominania λ w bloku filtra RLS Memory weighting factor z 1.0 na

np. 0.99. Przeprowadzić obliczenia i zaobserwować utratę stabilności algorytmu adaptacji.

H. (*) Przywrócić

λ=1, przełączyć Switch w położenie Signal gen. skrócić czas symulacji

(Stop time=500) i przeprowadzić symulacje z filtrem RLS z parametrami jak w punktach A,B,C.
&

Porównać wyniki uzyskane dla filtracji NLMS i RLS.

4. Opracowanie sprawozdania

W sprawozdaniu należy zawrzeć zarejestrowane wyniki eksperymentów numerycznych z odpowiednimi
opisami oraz wyjaśnieniami wskazanych problemów.