FILTRACJA CYFROWA

FILTRY O SKOŃCZONEJ ODPOWIEDZI IMPULSOWEJ SOI (FIR)

Operacja splotu dyskretnego

y(n) = h(n)



x(n)

h(n)

– odpowiedź impulsowa układu (filtru) na wymuszenie



(t)

Równanie filtru typu FIR

b

k

-

współczynniki filtru – próbki odpowiedzi impulsowej filtru w odstępach

czasu równych okresowi próbkowania

Struktura bezpośrednia

z

-1

+

z

-1

z

-1

+

+

z

-1

y(n)

x(n)

x

x

x

x

b

0

b

1

b

2

b

N-1

 

















1

0

N

k

k

k

n

x

b

n

y

Z

-1

– element opóźniający – jeden okres próbkowania

x(n-1)

x(n-2)

x(n-N-1)

Wymagania na realizację operacji splotu dyskretnego



Przechowywanie N próbek sygnału wejściowego



Przechowywanie N próbek odpowiedzi impulsowej filtru (współczynników)



Konieczność organizacji 2 buforów o rozmiarze N



Wymaganie N mnożeń i N dodawań



Wymaganie N-

1 przesunięć próbek w buforze sygnału



Obliczenia wykonywane w sposób ciągły, dla każdej próbki wejściowej

algorytm realizacji

x(0)







+

+

y(n)

=

h(0)

h(1)

h(2)

h(N-1)

+

+



x(n) x(n-1)

x(n-2)

x(n-N-1)

…

…

…

x(1)

x(2)

…

Równanie filtru II rzędu

y(n) = h(0)



x(n) + h(1)



x(n-1) + h(2)



x(n-2)

Odpowiedź na wyjściu filtru dla kolejnych chwil czasu n

y(0) = h(0)



x(0) +

h(1)



0 +

h(2)



0

y(1) = h(0)



x(1) +

h(1)



x(0) +

h(2)



0

y(2) = h(0)



x(2) +

h(1)



x(1) +

h(2)



x(0)

y(3) = h(0)



x(3) +

h(1)



x(2) +

h(2)



x(1)

Organizacja pamięci danych – bufory liniowe



Po każdym cyklu obliczeń przesunięcie próbek w buforze sygnału



Najstarsza próbka usuwana z bufora



Bieżąca próbka sygnału x(n) wchodzi na pozycję najmłodszej



Pozycja próbek sygnału w buforze wzg. kolejności starszeństwa stała



Najmłodsza, najstarsza próbka sygnału pod tym samym adresem

x(n-2)

x(n-1)

z

-1

+

z

-1

+

y(n)

x(n)

x

x

x

h(0)

h(1)

h(2)

float R_in[128]; /* Input samples R_in[0] most recent, R_in[127] oldest. */

float h[]=

/* Impulse response of FIR filter*/

{

-0.0001,0.0003,0.0004,0.0002,-0.0002,-0.0005,-0.0004,0.0001,

0.0006,0.0006,0.0001,-0.0006,-0.0010,-0.0005,0.0005,0.0013,

0.0010,-0.0002,-0.0016,-0.0017,-0.0003,0.0016,0.0025,0.0012,

-0.0013,-0.0032,-0.0025,0.0006,0.0036,0.0039,0.0007,-0.0036,

-0.0054,-0.0027,0.0029,0.0068,0.0052,-0.0012,-0.0076,-0.0081,

…
-0.0015,0.0074,0.0112,0.0055,-0.0059,-0.0139,-0.0108,0.0026,

0.0159,0.0173,0.0033,-0.0165,-0.0254,-0.0129,0.0145,0.0355,

0.0291,-0.0075,-0.0507,-0.0623,-0.0141,0.0897,0.2093,0.2890,

0.2890,0.2093,0.0897,-0.0141,-0.0623,-0.0507,-0.0075,0.0291,

0.0355,0.0145,-0.0129,-0.0254,-0.0165,0.0033,0.0173,0.0159,

}

float fir_filter (float sample)

{

int i;

float acc=0;

float prod;

R_in[0] = (float) sample;

/* Update most recent sample */

acc = 0;

/* Zero accumulator */

for (i=0; i<128; i++)

/* 128 taps */

{

prod = (h[i]*R_in[i]);

acc = acc + prod;

}

for (i=127; i>0; i--)

/* Shift delay samples */

R_in[i]=R_in[i-1];

return (acc);

}

organizacja programu

Przejście do buforów kołowych



Pomysł!

Zamiast przesuwać próbki w buforze liniowym wykonać przemieszczenie

próbek w buforze wzg. kolejności starszeństwa



Nie zmienia się fizyczne położenie próbki w buforze, ale zmienia się jej

kolejność starszeństwa



Przy stałym rozmiarze przez bufor przechodzi zbiór N próbek, gdzie liczba

poziomów starszeństwa jest stała



Zmiana kolejności starszeństwa próbek w buforze musi być wykonana w

zamkniętym cyklu rotacji bufor kołowy w miejsce liniowego



Zmiana kolejności starszeństwa poprzez manipulację wskaźnikami

adresowymi



Manipulacja wskaźnikami adresowymi z wykorzystaniem adresowania

kołowego

Rotacja bufora kołowego



Bieżąca próbka sygnału wchodzi na pozycję najmłodszej w buforze

w danym cyklu obliczeń n



Przy rotującym buforze jest to miejsce najstarszej próbki w buforze

w poprzednim cyklu obliczeń n-1



Po zapisaniu bieżącej próbki inkrementacja wskaźnika adresowego próbek



Przed rozpoczęciem obliczeń, po zapisaniu bieżącej próbki sygnału,

wskaźnik ustawiony na adres najstarszej próbki w buforze

cykl obliczeń

algorytm realizacji z adresowaniem

kołowym

.start "response",0x809800

.start "input",0x809804

.start "code",0x809810

.sect "response"
H

.float 1.0 ; h(2)

.float 0.5 ; h(1)

.float 0.25 ; h(0)

.sect

"input"

X

.float 0.0 ; x(n-2)

.float 0.0 ; x(n-1)

.float 0.0 ; x(n)

.sect "code"

ADR_H

.word

H

ADR_X

.word

X + 2

; x(n)

N

.set 3

SAMPLE .set 0x809000
OUTPUT .set 0x809C01

LDI N, BK

LDI @ADR_H, AR0

; AR0 h(2)

LDI @ADR_X, AR1

; AR1 x(n)

FILTER LDF @SAMPLE, R3

STF R3, *AR1++%

; AR1 x(n-2)

LDF 0, R0

LDF 0, R2

RPTS

N

MPYF *AR0++%, *AR1++%, R0

|| ADDF R0, R2, R2

ADDF R0, R2

STF R2, @OUTPUT

B

FILTER

FILTRY O NIESKOŃCZONEJ ODPOWIEDZI IMPULSOWEJ NOI (IIR)

Równanie filtru typu IIR




Współczynniki filtru :

b(n)

– zera transmitancji filtru,

a(m)

– bieguny transmitancji filtru.

Struktura bezpośrednia

Równanie filtru II rzędu

y(n) = b(0)



x(n) + b(1)



x(n-1) + b(2)



x(n-2) +

+ a(1)



y(n-1) + a(2)



y(n-2)

wymagania realizacji

x(n)

+

b

0

+

b

1

+

b

2

+

b

N-1

+

b

N

+

+

+

+

+

a

1

a

2

a

M-1

a

M

z

-1

y(n)

z

-1

z

-1

z

-1

z

-1

z

-1

x(n-N)

x(n-2)

x(n-1)

y(n-1)

y(n-2)

y(n-M)

 

  



  



















M

k

N

k

k

n

y

k

a

k

n

x

k

b

n

y

1

0

Równania filtru - struktura zmodyfikowana







Struktura sekcji bikwadratowej

Równania różnicowe sekcji

d(n) = x(n)

+

a(1)



d(n-1)

+

a(2)



d(n-2)

y(n) = b(0)



d(n) +

b(1)



d(n-1)

+

b(2)



d(n-2)

wymagania realizacji

   















N

k

k

k

n

d

a

n

x

n

d

1

 













N

k

k

k

n

d

b

n

y

0

+

Z

-1

+

+

+

X

b(0)

X

b(1)

X

b(2)

X

a(1)

X

a(2)

Z

-1

x(n)

y(n)

d(n)

d(n-1)

d(n-2)

Równania różnicowe filtru - struktura zmodyfikowana (sekcja bikwadratowa)

d(n) = x(n)

+

a(1)



d(n-1)

+

a(2)



d(n-2)

y(n) = b(0)



d(n) +

b(1)



d(n-1)

+

b(2)



d(n-2)

Odpowiedź na wyjściu filtru dla kolejnych chwil czasu n

d(0) = x(0)

+

a(1)



0

+

a(2)



0

y(0) = b(0)



d(0) +

b(1)



0

+

b(2)



0

d(1) = x(1)

+

a(1)



d(0) +

a(2)



0

y(1) = b(0)



d(1) +

b(1)



d(0) +

b(2)



0

d(2) = x(2)

+

a(1)



d(1) +

a(2)



d(0)

y(2) = b(0)



d(2) +

b(1)



d(1) +

b(2)



d(0)

d(3) = x(3)

+

a(1)



d(2) +

a(2)



d(1)

y(3) = b(0)



d(3) +

b(1)



d(2) +

b(2)



d(1)

Organizacja pamięci danych – bufory liniowe



Po każdym cyklu obliczeń przesunięcie próbek w buforze opóźnień d



Najstarsza próbka usuwana z bufora



Bieżąca próbka d(n) wchodzi na pozycję najmłodszej



Pozycja próbek d w buforze wzg. kolejności starszeństwa stała



Najmłodsza, najstarsza próbka pod tym samym adresem

n = 0

n = 1

n = 2

n = 3

Realizacja operacji filtru

Bufor węzłów opóźnień dla dalszych chwil czasu

x(n)

a(1)

a(2)

d(n) d(n-1)

d(n-2)

b(0)

b(1)

b(2)





+

+

d(n)

=





+

=



+

y(n)

d(2)

d(1)

d(0)

d(3)

d(2)

d(1)

d(0)

n = 2

n = 3

x(0)

a(1)

a(2)

d(0)

0

0

b(0)

b(1)

b(2)





+

+

d(0)

=





+

=



+

y(0)

n = 0

x(1)

a(1)

a(2)

d(1)

d(0)

0

b(0)

b(1)

b(2)





+

+

d(1)

=





+

=



+

y(1)

n = 1

algorytm realizacji

float a[a1, a2];
float b[b0, b1, b2];
float d[d0, d1, d2];

float IIR_filter (float sample)
{

float xn, y0;

y0 = 0;

d[0] = 0;

xn = (float) sample;

d[0] = xn + d[1]*a[1] + d[2]*a[2];

// calculate current d[0]

y0 = d[0]*b[0] + d[1]*b[1] + d[2]*b[2];

// calculate current output

d[2] = d[1];

// shift delay samples

d[1] = d[0];

return (y0);

}

Przejście do buforów kołowych



Pomysł!

Zamiast przesuwać próbki w buforze liniowym wykonać przemieszczenie

próbek w buforze wzg. kolejności starszeństwa



Nie zmienia się fizyczne położenie próbki w buforze, ale zmienia się jej

kolejność starszeństwa



Przy st

ałym rozmiarze przez bufor przechodzi zbiór N próbek, gdzie liczba

poziomów starszeństwa jest stała



Zmiana kolejności starszeństwa próbek w buforze musi być wykonana w

zamkniętym cyklu rotacji bufor kołowy w miejsce liniowego



Zmiana kolejności starszeństwa poprzez manipulację wskaźnikami

adresowymi



Manipulacja wskaźnikami adresowymi z wykorzystaniem adresowania

kołowego

Rotacja bufora kołowego



Przed rozpoczęciem obliczeń wskaźnik bufora

d

ustawiony na adres

najstarszej próbki w buforze



Bieżąca próbka d, obliczona w tym cyklu, wchodzi na pozycję najmłodszej

w buforze w danym cyklu obliczeń n



Przy rotującym buforze jest to miejsce najstarszej próbki w buforze

w poprzednim cyklu obliczeń n-1



Po zakończeniu cyklu obliczeń dekrementacja wskaźnika adresowego

bufora

d

Organizacja pamięci danych – bufor cyrkularny

Algorytm realizacji z adresowaniem

kołowym:

1. Ustaw AR0, AR1

2. Pobierz

próbkę sygnału x(n)

3. Wyznacz a(2)



d(n-

2), zwiększ AR0

4. Wyznacz b(2)



d(n-

2), zwiększ AR0 i AR1

5. Wyznacz a(1)



d(n-1)

Wyznacz [a(2)



d(n-

2)] + x(n), zwiększ AR0

6. Wyznacz b(1)



d(n-1)

Wyznacz d(n), zwiększ AR0 i AR1

7. Zapisz d(n) do bufora cyrkularnego

8. Wyznacz b(0)



d(n)

Wyznacz [b(2)



d(n-2)] + [b(1)



d(n-1)

],

zmniejsz AR1

9. Wyznacz y(n)

cykl obliczeń

a(2)

b(2)

a(1)

AR0

0

0

d(0)

AR1

d(1)

0

d(0)

AR1

d(1)

d(2)

d(0)

AR1

d(1)

d(2)

d(3)

AR1

n = 0

n = 1

n = 2

n = 3

d(n-2)

d(n-1)

d(n)

AR1

AR0

– rejestr indeksowy bufora współczynników

w

artość początkowa



adres a(2).

AR1

– rejestr indeksowy bufora węzłów opóźnień

w

artość początkowa



adres d(n-2).

d(n)

b(1)

b(0)

.start "response”, 0x809800

.start "delay", 0x809808

.start "code", 0x809815

.sect "response"
H

.float 0.95

;a(2)

.float 0.7

;b(2)

.float 0.85

;a(1)

.float 0.6

;b(1)

.float 0.5

;b(0)

.sect

"delay"

D

.float 0.0

;d(n-2)

.float 0.0

;d(n-1)

.float 0.0

;d(n)

.sect "code"

ADR_H

.word

H

ADR_D

.word

D

SAMPLE .set 0x809000
OUTPUT .set 0x809C01

LDI 3, BK1

LDI 5, BK2

LDI @ADR_D, AR1

; d(n-2)

LDI @ADR_H, AR0

; a(2)

FILTER LDF @SAMPLE, R2

; x(n)

MPYF *AR0++%, *AR1, R0

; a(2)



d(n-2)

MPYF *AR0++%, *AR1++%, R1 ; b(2)



d(n-2)

MPYF *AR0++%, *AR1, R0

; a(1)



d(n-1)

|| ADDF R0, R2, R2

; a(2)



d(n-2)+x(n)

MPYF *AR0++%, *AR1++%, R0 ; b(1)



d(n-1)

|| ADDF R0, R2, R2 ; a(2)



d(n-2)+x(n)+a(1)



d(n-1)

STF R2, *AR1

; zapisz d(n)

MPYF *AR0++%, *AR1--%, R0 ; b(0)



d(n)

|| ADDF R0, R1, R2 ; b(2)



d(n-2)+b(1)



d(n-1)

ADDF R2, R0, R0

; y(n)

STF R0, @OUTPUT

; zapisz d(n)

B

FILTER

Dla N sekc

ji bikwadratowych połączonych kaskadowo

Równania różnicowe filtru

Dla kolejnych sekcji i

i = 0

d(0,n) = x(n)

+

a1(0)



d(0,n-1) + a2(0)



d(0,n-2)

y(0,n) = b0(0)



d(0,n) + b1(0)



d(0,n-1) + b2(0)



d(0,n-2)

i = 1

d(1,n) = y(0,n)

+

a1(1)



d(1,n-1) + a2(1)



d(1,n-2)

y(1,n) = b1(1)



d(1,n) + b1(1)



d(1,n-1) + b2(1)



d(1,n-2)

…

i = N-1

d(N-1,n) = y(N-2,n) + a1(N-1)



d(N-1,n-1) + a2(N-1)



d(N-1,n-2)

y(N-1,n) = b1(N-1)



d(N-1,n) 0 + b1(N-1)



d(N-1,n-1) + b2(N-1)



d(N-1,n-2)

SB(0)

x(1,n)

x(n)

y(0,n)

x(N-1,n)

y(1,n)

SB(1)

SB(N-1)

y(N-1,n)

. . .

Organizacja pamięci danych – bufory cyrkularne



























AR0

– rejestr indeksowy bufora współczynników;

Wartość początkowa



adres a2(0).

AR1

– rejestr indeksowy kolejnych buforów węzłów opóźnień;

Wartość początkowa



początek pierwszego bufora - adres d(0,n-2).

Dla kolejnych buforów i (0, N-1)

adres d(i,n-2) = adres d(0,n-2) + 4



i

a2(0)

AR0

b2(0)

a1(0)

b1(0)

b0(0)

.

.

.

a2(N-1)

b2(N-1)

a1(N-1)

b1(N-1)

b0(N-1)

d(0,n-2)

AR1

d(0,n)

d(0,n-1)

d(0,n)

d(N-1,n-2)

AR1

d(N-1,n)

d(N-1,n-1)

d(N-1,n)

.

.

.

adr0+4(N-1)

adr0

0

0

d(0,0)

AR1

n = 0

d(0,1)

0

d(0,0)

AR1

n = 1

d(N-1,0)

AR1

n = 0

0

0

d(N-1,0)

AR1

n = 1

0

d(N-1,1)

Algorytm działania:

1.

Pobierz próbkę sygnału wejściowego

2.

Ustaw AR0 na początek bufora współczynników

3.

Ustaw AR1 na początek bufora kołowego sekcji 0

– adres węzła opóźnień d(0,n-2) = adres 0

4. Wyznacz d(0,n) oraz y(0,n)

– obliczenia dla sekcji 0

5. Obliczenia dla sekcji i-tej

6.

Ustaw AR1 na początek bufora kołowego sekcji i
– adres węzła opóźnień d(i,n-2) = adres 0 + 4



i

7. Wyznacz d(i,n) oraz y(i,n)

8. Czy wykonane obliczenia dla wszystkich sekcji

9.

Zakończone obliczenia y(N-1,n) – powrót do pkt. 1.

N

T

FILTRY ADAPTACYJNE

Filtr adaptacyjny, to filtr,

którego współczynniki uaktualniane są przez

algorytm adaptacyjny

dla optymalizacji odpowiedzi filtru ze względu na

przyjęte kryterium.

Ogólna budowa filtru adaptacyjnego

Algorytm adaptacyjny zmienia współczynniki filtru (wagi) dla minimalizacji
błędu e(n) pomiędzy sygnałem wyjściowym y(n) a sygnałem d(n).

Równanie filtru klasycznego FIR

y(n) = h(0)



x(n) + h(1)



x(n-1) + h(2)



x(n-2) + ... + h(N-1)



x(n-(N-1))

Równanie filtru adaptacyjnego

y(n) = h(n,0)



x(n) + h(n,1)



x(n-1) + h(n,2)



x(n-2) + ... + h(n,N-1)



x(n-(N-1))

Popularny algorytm adaptacyjny - LMS Least Mean Square.

Przy zastosowaniu algorytmu LMS współczynniki filtru są uaktualniane
według równania

h(n+1,i) = h(n,i) + 2mu



e(n)



x(n-i)

mu

– krok adaptacji.

ALGORYTM

ADAPTACYJNY

STRUKTURA

FILTRU

e(n) = d(n)

– y(n)

Dla

przykładowego filtru 2 rzędu

h(n+1,2) = h(n,2) + 2mu



e(n)



x(n-2)

h(n+1,1) = h(n,1) + 2mu



e(n)



x(n-1)

h(n+1,0) = h(n,0) + 2mu



e(n)



x(n-0)

Najprostszy algorytm działania:

1.

Pobierz próbkę sygnału wejściowego x(n).

2.

Wyznacz sygnał wyjściowy filtru y(n).

3. Wy

znacz nowe współczynniki filtru wg. algorytmu

adaptacyjnego LMS.

4.

Zapisz nowe współczynniki do bufora cyrkularnego.

Organizacja pamięci danych

AR0

x(n-2)

x(n-1)

x(n)

AR1

x(n)

h(n,2)

h(n,1)

h(n,0)

h(n+1,2)

h(n+1,1)

h(n+1,0)

n

n+1

.start "response",0x809800

.start "input",0x809804

.start "code",0x809810

.sect "response"
H

.float 1.0

.float 0.5

.float 0.25

.sect

"input"

X

.float 0.0

.float 0.0

.float 0.0

.sect "code"

ADR_H

.word

H

ADR_X

.word

X + 2

N

.set

3

DWAXMU .float 0.1
SAMPLE .set

0x809000

OUTPUT .set

0x809C01

LDP @0x800000

LDI N, BK

LDI @ADR_H, AR0

LDI @ADR_X, AR1

FILTER
;
; czytaj próbki x(n) oraz d(n)
;

LDF @SAMPLE, R7

STF R7, *AR1++%

LDF @SAMPLE + 1, R7

;
; wyznacz odpowiedź y(n)
;

LDF

0, R2

RPTS

N-1

MPYF3

*AR0++%, *AR1++%, R0

|| ADDF3

R0, R2, R2

ADDF

R0, R2

STF

R2, @OUTPUT

;

; wyznacz i zapisz współczynniki h(n+1,i))

;

SUBF

R2, R7

;R7=e(n)=d(n)-y(n)

STF

R7, @OUTPUT + 1

MPYF

DWAXMU, R7

;R7=2*mu*e(n)

LDF

0, R2

MPYF3

*AR1++%, R7, R0 ;R0=x(n-2)*2mue(n)

LDI

N-2, RC

RPTB

LMS

MPYF3

*AR1++%, R7, R0 ;R0=x(n-i-1)*2mue(n)

|| ADDF3

*AR0, R0, R2

;R2=h(n,i)+

; +x(n-i)*2mue(n)

LMS

STF

R2, *AR0++%

;h(n+1,i)=R2

ADDF3

*AR0, R0, R2

STF

R2, *AR0++%

B

FILTER