Zastosowanie nieliniowej filtracji Kalmana do tłumienia zakłóceń sygnału EKG

(1)

ZESZYTY NAUKOWE POLITECHNIKI ŚLĄSKIEJ Seria: AUTOMATYKA z. 103

_________1991 Nr kol. 1090

Jacek ŁĘSKI

ZASTOSOWANIE NIELINIOWEJ FILTRACJI KALMANA DO TŁUMIENIA ZAKŁÓCEŃ SYGNAŁU EKG

Streszczenie. W pracy przedstawiono nową metodą tłumienia zakłóceń sygnału EKG, których widmo częstotliwościowe nakłada się na widmo sygnału dla szerokiego zakresu częstotliwości. Sygnał EKG modelowano procesem autoregresyjnym (AR). Jednoczesna estymacja wartości próbek sygnału i nieznanych parametrów autoregresji prowadzi do problemu n i e  liniowego. Zagadnienie to zostało rozwiązane przez zastosowanie r oz

szerzonego filtru Kalmana i metody GPIA. Metoda może znaleźć zastosowa

nie w tłumieniu zakłóceń elektrokardiogramów wysiłkowych, EKG płodu i innych.

1. USTĘP

'¡Jednym z p ierwszych etapów komputerowej analizy sygnału EKG jest tłumienie zakłóceń. Zakłócenia których widmo częstotliwościowe nie nakłada się, lub nakłada się na widmo sygnału EKG dla wąskiego zakresu częstotliwości mogą być tłumione za pomocą klasycznych filtrów częstotliwościowych [2,5,26]. Metodą tą nie można bez zniekształcania sygnału (utrata informacji diagnostycznej) tłumić zakłóceń, których widmo częstotliwościowe nakłada się na widmo sygnału dla szerokiego zakresu częstotliwości [25].

Do zakłóceń takich należą np. :

- zakłócenia mięśniowe podczas elektrokardiograficznych prób wy s iłkowych [13, 16],

- zakłócenia powstające przy przesyłaniu sygnału EKG linią telefoniczną [6],

" zakłócenia EKG płodu podczas akcji porodowej [14],

-zakłócenia wysokorozdzielczych elektrokardiogramów [1,15].

Powszechnie stosowaną metodą tłumienia tych zakłóceń jest selektywne uśrednianie sygnału EKG z wielu ewolucji serca [13,25]. Metoda ta bazuje na

s

dyskusyjnym założeniu idealnej powtarzalności kształtu sygnału w każdej ewolucji serca [12].

(2)

182 J . Ł ę s k i

Stosując w praktyce tę metodę uzyskujemy uśredniony kształt sygnału. W p ewnych zastosowaniach, np. analizie obniżenia odcinka S T podczas próby wysiłkowej nie stanowi to problemu. Jednak nie jest możliwe zastosowanie tej metody np. do analizy zmienności odcinka QT podczas tej próby.

U artykule przedstawiona jest metoda tłumienia zakłóceń sygnału EKG oparta na modelowaniu sygnału procesem autoregresyjnym [7] oraz jednoczesne1 estymacji wartości próbek sygnału i nieznanych parametrów autoregresji.

Jeżeli przedstawimy ten model w postaci równań zmiennych stanu, włączając nieznane parametry do wektora stanu, to uzyskamy nieliniowy model względem wektora stanu/parametrów.

Klasycznymi metodami rozwiązania tego problemu są, np. :

- metoda porogramowania nieliniowego; estymacja parametrów metodą największą największej wiarygodności 14,21],

- estymacja adaptacyjna za pomocą przetwarzania równoległego (3), - stosowanie rozszerzonego filtru Kalmana [9],

- stosowanie dyskretnej zasady maksimum lub dyskretnego równania Eulera-La- grange'a [11],

W pracy zastosowano metodę GPIA (Generalized Partitioned Identification Algorithm) Lainiotisa [18,19], Nieznane statystyki szumów estymowano metodą Mayersa 123).

2. METODA

Sygnał EKG modelujemy procesem autoregresyjnym rzędu p [17,22]:

P

(1) 1=1

gdzie:

a. - nieznane parametry,

w - gaussowski szum biały spełniający założenia:

V Elw(k)] = 0

k

⁽²⁾

V Htw(l)w(JJJ = Q (i)3

1. j i. j'

5 - delta Kroneckera.

(3)

Zastosowanie nieliniowej filtracji Kajmana do 183

Zakłócony sygnał EKG otrzymujemy dodając gaussowski szum biały:

z(k) = x(k) + v(k) (3)

Szum v spełnia założenia:

V E [ v (k )) = 0

k (

1

)

V E[v(i )v( j ) ] - R U J S .

1, j ’J

Równania (1), (3) możemy przedstawić w formie nieliniowego modelu zmiennych stanu:

x (k+1) = f[x (k)] + G w(k)

a a a

z(k) = H x (k) + v(k) a a

i 5 a !

(5b)

gdzie-

wektor stanu/parametrów, (b)

f [ x a (k)l

x ( k ) x (k)

P P + 1

x ,( k )+ x (k)x _(k)

₁ _p _p+2

. (k)+x (k)x_ (k)

p-i p 2p

p+1

zp

r r lX

en

(4)

184 ■J. Łęski

H = [O. . ,0i:0. . .0] = [H : H J ;

a . x. 0 G =

a ( 8 )

Zakładamy ponadto:

V E [ v ( i )w (j )]=0 i. J

V E[v(i)x(j)]=0 i

(9)

V E[w(i)x(0)]=0 i

gdzie:

x(0) - stan początkowy.

Rozwijamy f[x (k)] w szereg Taylora w otoczeniu x (k/k)

a a

f[x (k)]=f[x (k/k)]+F (k)[x (k)-x (k/k)]+...

a a a a a (

10

⁾

Pomijając człony wyższego rządu i zakładając znajomość x^(k/k) otrzymujemy:

x (k+l)=F (k)x (k)+G w(k)*u(k)

a a a a

z(k)=H x (k)+v(k) a a

(lla)

(1 lb)

gdzie:

u(k)=x (k/k)-F (k)x (k/k)

a fi 3 (12)

(5)

Zastosowanie nieliniowej filtracji Kajmana do. 185

5 f F (k) =

a 5x

x (k/k) a

0 0 .

...

0 xp+l 0 1 , . , 0

X p+2

Ix P 0 . , . .

__

¹

X2p

0 I

Fx (k) F 0 (k)

(13)

Stosując do estymacji x & filtr Kalmana otrzymujemy [9,24]:

x (i/i)=x (i)+K (i)z(i)

a a a

x (i)=f[x (i-l/i-1)]

a a

z(i)=z(i)-K x (i) 3 3

K (i)=M ( 1 )HTB 1 (1)

a a a a

B (i)=H M (i)HA+R(i) a a a a

M (i)=F (i-l)P (i-l)F (i-1)T +G Q(i-1)GT

a a a a a a

P (i)=[I-K (i)H ]M (i)

a a a a

( 1 4 a )

( 14 b)

( 1 4 c )

( 14 d)

(14e)

( 1 4 f )

( 1 4 g )

gdzie: x (i/i) - estymata wartości x (i) na podstawie pierwszych 1-pomiarów,

^ XV *rv X» ^

P (i)=E{[x (i/i)-x (i)]‘[x (i/i)-x (1)1 > - macierz kowariancji błędu

a a a a a

e s t y ma c ji .

(6)

1S6

^{J. Łgsk}

W dalszej części pracy zakładamy stacjonarnosć zakłóceń:

’ Q(i)=Q

V R (i )=R 1

(14h!

(14i)

Zgodnie z koncepcją l.ainiotisa zakładamy, że stan początkowy jest sumą dwódi niezależnych wektorów gaussowskich [18]:

x(0) 0

= +

0 0(0)

(15

gdzie:

- nominalna część stanu początkowego, x r (C) - pozostała część stanu początkowego.

M acierz kowariancji stanu początkowego ma postać:

P (0)=P (0)+P (0) =

a n r

v ° ) ^{0 '} ^{' 0} ⁰

0 0 0

po e (0)

(16:

Na podstawie (14),(15),(16) widzimy, że wzmocnienie filtru Kalmana ma postać

K (i) = a

K U )

(17

Przy poczynionych założeniach filtr Kalmana jest jedynie e st ymatorem stanu 2I parametry pozostają niezmienione. Stąd równania (14a)-(14g) możem y przepisać w postaci wymagającej mniejszego nakładu obliczeniowego:

x (i/i)=x (i)+K (i)z(i)4X X X *

xx (i)=fx [xa ti_1/i" 1)]

z(i)=z(i)-(I x (i)X X

(1 Ga)

(18b

(18c

(7)

Zastosowanie nieliniowej filtracji Kalmana do. 187

K (i)=M C i)HT B _ 1 (i )

x X X X

B (i)=H M (i)H _ x x x x+R

(18d)

(18e)

Mx (i)=Fx (i-l)Px (i-l)Fx (i-l)T +G xQ Gj

Px (i)=[I-Kx (i)Hx )Mx (i)

(1 8 f )

Cisg)

Za pomocą filtru Kalmana prowadzimy "w przód" estymacją stanu, zakładając że 0 jest znane. "W tył" estymujemy nieznane parametry za pomocą wygładzania u stałym punkcie [10]:

6 =0(O)+P“g M 0 (N,))

pe e = [0e e (N-0 )+C ( 0 ) r l

(19)

( 2 0 )

gdzie:

M (N,0) = n

M x (N 0)"

M 0 (N 0) ■

z

2> - i=i e

1 )FT (i-l)HTB 1 (i)z(i)

cl 3 X (

2 1

⁾

O (N , 0) = n

0 xx(N.°)

OX 0 ('O)

O ^ N . O ) ° 0 0 (N.°)

= ^ Zn( i -1 )FT (i-l)HT B _ 1 (i)H F (i-l)Z (i-l)

a a x a a n

1=1 gdzie:

N - liczba próbek sygnału.

( 2 2 )

2 (i )=[I-K (i )H ]F (i-l)Z (i-l); 2 (0)=I

n a a a n n

(23)

Nieznane wariancje szumów estymujemy za Mayersem [23):

N R = — V1-N /

J=1

^{v(j)2 -}

N-l

N H M (j ) H

X X X (24)

(8)

188 J. Łęski

L

n 2 iN_1

w (j) - i - i r

J-l

w (J)=Gx (x(J/J)-Fx (J-l)x(J-l/j-l)]

Estymacją parametrów prowadzimy iteracyjnie podstawiając:

0

»

⁰

(

⁰

)

P ---- » P (0)

09 00

R » R

(25)

(26)

(27)

(28)

Q —

Do inicjałizacji obliczeń używamy:

x^(0) = (z(l)... z(l):

Px (0)=I*R

P ^ ( 0 ) = 0

(29)

(30)

(31)

Z astosowano heurystyczną regułą mającą na celu uzyskanie zbieżności algorytmu: w kolejnych iteracjach estymowano na przemian statystyki szumów i prametry autoregresji. Estymaty statystyk szumów wyznaczano na podstawie równań (24)-(27), gdy spełniony był warunek:

0-0(0) <Cj, - parametr. (32)

Następnie znów estymowano parametry auto-ogresji. Iteracje przerywano gdy spełnione były warunki:

¡R-R|<e,, “ ¡Q-Q|<c3 : - parametry. (33)

Zbieżność algorytmu uzyskiwano w 9 do 10 iteracjach.

(9)

Zastosowanie nieliniowej filtracji Kajmana do. 189

3. REZULTATY

Z prac Eisensteina, Vaccaro [8] oraz Lina, Changa [20] wynika, że zwiększenie rzędu autoregresji powyżej p=2 powoduje nieznaczne zmniejszenie błędów modelowania sygnału EKG. Powoduje to natomiast w przypadku stosowania algorytmu opisanego w artykule zwiększenie nakładu obliczeniowego

3

proporcjonalnego do (2p) . Z tego względu do testowania opisanej metody stosowany był model autoregresyjny rzędu 2 AR(2). Do testowania używano przebiegów E KG próbkowanych z częstotliwością fp=500 Hz i kwantowanych 12 bitowo (lpV na kwant). Z bazy danych elektrokardiograficznych wybrano 20 odcinków po 200 próbek, zawierających istotne z diagnostycznego punktu widzenia fragmenty: załamki P, zespoły QRS, odcinki S T oraz załamki T.

Przeprowadzono następujące badania symulacyjne: oryginalny sygnał EKG zakłócano szumem białym, następnie tłumiono zakłócenia i porównywano ten sygnał z oryginalnym.

Stosunek sygnał-szum definiowano jako:

p

SNR [dB] := 10 log -^5? (34)

n

gdzie:

^o r s ~ moc sysnału EKG*

P - moc zakłóceń, n

Testowano równolegle następujące metody tłumienia zakłóceń:

- "klasyczną" częstotliwościową filtrację cyfrową; filtr dolnoprzepustowy o częstotliwości granicznej 10 Hz,

- uśrednianie sygnału w dziedzinie czasu z 16 ewolucji serca, - metodę opisaną w artykule.

Rysunek 1 przedstawia zależność wyjściowego stosunku sygnału do szumu CSNRout) od wejściowego ( S N R ^ ) dla testowanych metod. Są to wyniki uśrednione d la 20 odcinków sygnału EKG używanych do testowania i 10 realizacji szumu zakłócającego dla każdego odcinka sygnału. Przykładowe tłumienia zakłóceń opisanymi metodami przedstawiają rysunki od 2 do 4 (fragment sygnału EKG z załamkami P i początkiem zespołu QP.S). Rysunki te przedstawiają w części a sygnał zakłóceniowy (SNR=3 dB), w części b sygnał po filtracji (linia przerywana) 1 oryginalny (linia ciągła).

(10)

190 J. Ł ę s k i SNftout CdB]

Rys. 1. Zależność stosunku sygnału do szumu sygnału po filtracji od stosunku sygnał-szum sygnału filtrowanego:

uśrednianie statystyczne 16 okresów sygnału, - - m etoda opisana w pracy, - filtracja cyfrowa

Fig. 1. Output signal to noise ratio as a function of input signal to noise ratio:

statistical averaging (16 cycles), - - Kalman filtering, - low-pass digital filtering

Na podstawie rysunku 1 widzimy, że metoda opisana w pracy wykazuje zależność stopnia tłumienia zakłóceń od S K R sygnału wejściowego. Metoda uśredniania tłumi wariancję zakłóceń k-razy; gdzie k jest liczbą uśrednianych okresów sygnału. Jest to oszacowanie optymistyczne wymagające spełnienia następujących założeń: stacjonarności zakłóceń oraz idealnego centrowania uśrednianych okresów. Metoda opisana w pracy nie wymaga spełnienia powyższych założeń, gdy ż estymuje wartości próbek na podstawie jednego o k re su sygnału.

"Klasyczna" częstotliwościowa filtracja cyfrowa powoduje niedopuszczalne zniekształcenia sygnału E3CG (usuwanie częstotliwości zawartych w sygnale).

(11)

Rys. 2. Tłumienia zakłóceń metodą częstotliwościowej filtracji cyfrowej:

a - sygnał zakłócony (SNR=3dB), b - sygnał oryginalny (linia ciągła), sygnał po filtracji (linia przerywana)

Fig. 2. Example of noise elimination for low-pass digital filtering:

a - noise signal (SNR=3dB), b - original and filtering signal

(12)

192 J. Łęski

Rys. 3. Tłumienie zakłóceń metodą opisaną w pracy:

a - sygnał zakłócony (SNR=3dB), b - sygnał oryginalny (linia ciągła), sygnał po filtracji (linia przerywana)

Fig. 3. Example of noise elimination for Kalman filtering:

Rys. 4. Tłumienie zakłóceń metodą uśredniania statystycznego:

a - sygnał zakłócony (SNR=3dB), b - sygnał oryginalny (linia ciągła), sygnał po filtracji (linia przerywana)

Fig. 4. Example of noise elimination for statistical averaging:

(13)

ZasLosowanie nieliniowej filtracji Kalwaria do. 193

4. PODSUMOWANIE

W pracy przedstawiona jest nowa metoda tłumienia zakłóceń sygnału EKG oparta na nieliniowej filtracji Kalmana. Jednoczesna estymacja stanu i parametrów realizowana jest za pomocą algorytmu GPIA. Sygnał EKG modelowany jest procesem autoregresyjnym rzędu 2. Metoda w przeciwieństwie do uśredniania statystycznego pozwala tłumić zakłócenia pojedynczego okresu sygnału EKG. Nie wymaga ona ponadto klasyfikacji i centrowania uśrednianych okresów sygnału. Metoda może znaleźć zastosowanie do tłumienia zakłóceń mięśniowych w elektrokardiogramach wysiłkowych, zakłóceń EKG płodu i w y s ok o  rozdzielczych elektrokardiogramów oraz innych.

I.ITERATURA

[1] H.A.M. Al-Nashash, : S. W. Kelly, D.J.E.Taylor: Beat-to-beat detection of His-Purkinje system signals using adative filters, Med & Biol. Eng. &

Comput., 1988, 26, 117-125.

[2] J.A.Alste, W . Eck, 0. E. Herrmann: ECG baseline wander reduction using linear phase filters, Computers & Biomedical Reseach 19, 417-427, (1986) [3] B. D. 0. Anderson, J.B.Moore: Filtracja optymalna, WNT, Warszawa 1984.

[4] F.Bartoli, S.Cerutti: A Kalman filter procedure for processing of electroencephalogram, Proceedings of ICASSP 1982, New York IEEE 1982, pp. 721-724.

[5] K.G. Beauchamp: Przetwarzanie sygnałów metodami analogowymi i cyfrowymi, WNT, Warszawa 1878.

16] A. S. Bersor.: Telephone transmission of electrocardiograms, Proceedings of the 2nd IFIP TC-4 Working Conference on Computer Application on ECG &

VCG analysis, Hannover, 11-14 October 1971, pp. 83-97.

[7] G.E. Box, G.M. Jenkins: Analiza szeregów czasowych - Prognozowanie i sterowanie, PWN, Warszawa 1983.

[8] B. A. Eisenstein, R.J.Vaccaro: Feature extraction by system identyfication, IEEE Trans.Sys.Men Cyb. SMC-12, N o . 1. Jan 1982, 15-23.

[9] 4.E.Elbert: Estimation and Control of Systems, Van Nostrand Reinhold Company Inc., New York 1984.

[10] B. J.Eulrich, D.Andrisani, D. G. Lainiotis: Partioning Identification Algorithms, IEEE Trans. Automatic Control, Vol. AC-25, No. 3, June 1980.

(14)

194 J . Łęski

[11] P.Eykoff: Identyfikacja w układach dynamicznych, PWN, Warszawa 1980.

[12] J.Groeben, R. S. Crow, A. D. Forbes: Some observations on ECG wavefors variability, Proceedings of the 2nd IFIP TC-4 Working Conference or Computer Application on ECG & VCG analysis, Hannover, 11-14 October 1971.

[13] P.Hsiaet al: An automated system for ST segment and arrhythmia analysis in exercise radionuclide ventriculography, IEEE Trans. Biomed. Eng., i BME-33, No. 6, June 1986, pp. 585-593.

[14] W.Jarisch, J. S. Detwiler: Statistical modeling of fetal heart rate variability, IEEE Trans. Biomed. Eng., BME-27, No. 10, Oct. 1980, pp.

582-589.

[15] S. Jesus, H.Rix: High resolution ECG analysis by an improved signal averaging method and comparison with a beat-to-beat approach, J. Biomed.

Eng. , Vol. 10, Jan. , pp. 25-32.

[16] N. B. Jones, et al: Spectral analysis and the interferance EMG, IEEE Proc., V o l . 129, No.9, Dec 1982, 987-992.

[17] R.L.Kashyap: Optimal feature selection and de cision rule in c lassification problem with time series, IEEE Trans. Inf. Theory 3, 281-307, (1978).

[18] D. G. Lainiotis: Partitioning: A unifying framework for adaptive systems, I: Estimation, Proc. IEEE, Vol.64, Aug. 1976.

[19] D. G. Lainiotis: Partitioning: A unifying framework for adaptive systems, II: Control, Proc. IEEE, Vol.64, Aug .1976.

120] K. P. Lin, W. H. Chang: ECG signal analysis by linear predictive method IEEE Ninth Annual Conference of the Engineering in Me dicine and Biolog Society 1987, pp. 557-558.

[21] K.Mańczak, Z.Nachorski: Komputerowa identyfikacja obiektów dynamicznych PWN, Warszawa 1983.

[22] J.P. Marques de Sa, C. Abren-Lima: A new ECG Classifier based on linear p rediction techniques, Computers & Biomedical R esearch 19, 213-22, (1986).

[2.3] R. A. Mayers, B. D.Tapley: Adaptive seqential estimation with unknown nois<

statistics. IEEE Trans. Automatic Control, Vol. AC-21, No .4, 1980.

A

[24] J.5.Meditch: Estymacja i sterowanie statystycznie optymalne w układać!

liniowych, WNT, Warszawa 1975.

(15)

[25] 0. Pahlm, L.Sornmo: Data processing of exercise E C G ’s, IEEE Trans.

Blomed. E n g . , BME-34, No2, Feb. 1987, pp. 158-165.

[26] A. Wojtkiewicz: Elementy syntezy filtrów cyfrowych, WNT, Warszawa 1984.

Recenzent: Prof, dr hab.inż. Ryszard Tadeusiewicz

Wpłynęło do Redakcji 30.03.1989 r.

nPHM EHEHilE H EJM H EłlH O il SHJIbTPAUH H no KAJIbMAHy nOAABJIEHHH II0M EX CRTHAjlA 3 K r

P e 3 b m e

B

c i a i h e onHcaHO npHM eHenne H ejiH H efiH oii (pm iBTpaijHH no K aiiB M aH y ą j u i.n o n a B jie H H a noM ex c x r H a jia

3Kr.

IIoM exH

hmcbt

c n e K T p , K O io p tiS H aK iia A H B ae T C s H a c n e K T p c a r H a j ia , H e c y m e ro c y ą e c i B e H H y B c £ H a rH o c - THHeCKOM TOHKH SpeHHH. H H $ 0 pM aiiHB H He HOTyT ÓHTB yC T paH eH H KJiaCCH

’ŁaCKHMH U,H$ipOBHMH $HJIBTpaM H . j U K HfleHTH$HKaHHH MO^ejIH C ilT H ajia n pH ueneH a a B T o p e r p e e c H B H a n m o a e jiB * IIo M e x a M o aejm p o BaH H ÓejiUM ra y c c o B H M myMOM.» OflHOBpeMeHHaH o n e iiK a cootohhhji n a p a M e ip o B s e ^ e T k HexHHeiiHoft n p o ó Jie M e . 3 i a npoOxeM a pem eHa npn n o M o ąa npHueneHHH paom upeH Horo ę n jiB T p a K a x B M a H a h M e ro n a G P IA . IIp e n c T a B jie H o c p a B H e a HO Boro w e io f la c x x a c c x n e c K H M ycp e n H e H n e M h H H S H e n p o n y c K H o t

$ H JiB T p a iy ie li. M e io fl MOKeT H a fiiH npuM eH eH H e, H a n p n M ep , b H CCJieno- BaHHHx $ H 3 in e c K o x n a r p K 3 0 K

, 3Kr

n jio .u a , h . t . a .

(16)

196 J. t q s k i

APPLICATION OF NONLINEAR KALMAN FILTERING FOR ECG NOISE REDUCTION

S u m m a r y

In the paper a nonlinear Kalman filter procedure is appliead for a reduction of the ECG noise. Such a noise has a bandwidth which overlaps the signal carrying the information useful from the clinical stand point and can not be removed by means of classical digital filtering. An autoregressive model is used for identification of the signal model. Noise is modeled by white gaussian sequence. Combined state and parameter estimation leads to a nonlinear problem.' The nonlinear problem is solved by the e xtended Kalman filter and generalized partitioned identification algorithm. The paper presents a comparison of the new method with classical averaging and low-pass filtering methods. For example the method may be applied in the ECG exercise, fetal ECG and others.