Analysis of biological rhythms using an interactive structured autoregressive modelling technique

Analysis of biological rhythms using an interactive structured autoregressive modelling technique

Computer Programs in Biomedicine 9 (1979) 1-12 © Elsevier/North-Holland Biomedical Press Section I: Programs, program packages and systems ANALYSIS O...

501KB Sizes 0 Downloads 23 Views

Computer Programs in Biomedicine 9 (1979) 1-12 © Elsevier/North-Holland Biomedical Press

Section I: Programs, program packages and systems ANALYSIS OF BIOLOGICAL RHYTHMS USING AN INTERACTIVE AUTOREGRESSIVE MODELLING TECHNIQUE

STRUCTURED

D.A. LINKENS and P.M. MORRISH Department of Control Engineering, University of SheffieM, Mappin Street, Sheffield $1 3JD, England

A method of spectral-type analysis for rhythmic biomedical data is presented based on autoregressive modelling techniques. The emphasis has been on portability between data-processing machines, so that the programmes axe written in FORTRAN 4 for the ICL 1900 series computers. The programmes are structured in blocks allowing for simple alteration to an interactive mode allowing the calculation and display of pre-filtering stages, model coefficients, the pole-zero plot, residuals, AR frequency spectrum, FFT of raw and filtered data. The method works on a small number of cycles and gives a direct read out of significant spectral components. Biological rhythms

Autoregressive modelling

1CL 1900

1. Introduction

EEG analysis [4] for the detection of non-stationary burst potentials. Initial attempts to use autoregressive modelling on gut rhythms have also been promising [5,6] and were performed using a data analysis package for industrial process identification [7]. The available options were not really suitable nor relevant to biomedical signal analysis and hence suitable sections of this package have been extracted, extended and restructured to form a self-contained, portable package capable of giving a high precision 'standard' (using an ICL 1906S machine) but capable of being transferred to smaller word-length machines for routine analysis. The basis of autoregressive modelling analysis is shown in fig. 1. It is assumed that the incoming data can be represented by a white noise signal operated on by a filter F(z). The error between the data and the filtered white noise is called the 'residuals' and this is used to readjust the parameters of the filter to give a minimum least-squares error. The resulting Filter F(z) is then factorised and the roots plotted on a 'z-plane diagram' as shown in fig. 1. Location of a pair of roots close to the unit radius circle indicates the presence of a rhythm frequency. The angle made by the roots to the origin equals wT, where ~o is the rhythm angular frequency and T is the data sampling period. Pairs of roots located in the first and fourth

The method and programmes developed and described in this paper have been motivated by the requirement for computer analysis of rhythmic potentials (called 'slow-waves') in the gastro-intestinal tract. Spectral analysis of these rhythms has been performed interactively using FFT methods for some years [ 1]. Although these are successful for gastric and small-intestinal slow-waves which are regular and present at all times, FFT analysis is not suited for large-intestinal rhythms. In the large intestine several rhythm frequencies can be present, usually one at a time, but they often only occur for a small number of cycles [2]. In such cases FFT is not suitable since the necessary frequency discrimination is not obtained, Autoregressive modelling appears to be a promising method of analysis for signals which occur for only a few cycles and are imbedded in non-white noise. Noise in medical signals is usually non-white, and in gut electrical rhythms there is a pJ'eponderance of spectral components at low frequencies in the region of 0.02 Hz. Autoregressive modelling is basically a least-squares method of data curve fitting with the model structured in difference equation or ztransform notation. It has been used extensively in industrial process modelling [3] and more recently to

1

2

D.A. Linkens, P.M. Morrish, Autoregressive modelling of biological rythms

DATA

~

:V7

FILTER -WHNoIsEITE Fig. la.

~RESlDUALS

;I

-T

F (•z) L. . . . . . . . . . .

, ~ _,

ADJUSTMENT

z-

PLAN E

PLOT

bits for the majority of the computation. The graph plotting is performed using a GHOST package developed by UKAEA, Culham Laboratory, for interfacing to a Calcomp plotter. The standard FFT calculations and matrix inversions use the NAG implementation, as does polynomial factorisation. When the AR model coefficients have been obtained and the z-polynomial factorised, the z-plane diagram is plotted giving a visual indication of spectral components. The closeness of the roots to the unit circle is related to their significance. The distance it.side the unit circle is calculated and listed as (1-Izl}. For a single pair of complex roots the relative amplitude is determined by the damping ratio (zeta) which itself is related to distance inside the unit circle. Lines of constant 'zeta' are logarithmic spirals in the z-plane and are given by the relationship:

Fig. lb. X

Fig. 1. Autoregressivemodellingfor spectral estimation. (a) Schematic diagram; (b)Typical z-plane plot.

~"- ~/1 + x 2 where

quadrant of the z-plane diagram usually indicate significant rhythms while roots in the second and third quadrant represent the 'coloured' spectrum of the non-white noise components. The method has been extended to give automatic printout of the significant rhythm components. This is a major advantage of the autoregressive technique over conventional FFT spectral analysis for which automatic peak detection represents a formidable difficulty in biomedical signals [8]. Another extension has been the calculation and display of the frequency spectrum represented by the filter F(z) which may then be compared to an FFT of the same stretch of data. Further facilities incorporate the removal of non-zero mean, pre-filtering of data, and tests for the 'whiteness' of the residuals via an autocorrelation function plot. Although all of these facilities would not be required for regular measurements using the method, they are of considerable benefit when validating the technique on a particular type of signal.

2. Computational methods The programmes have been developed on an ICL 1906S machine using single precision arithmetic of 48

x -

1 logl0(x/Re 2 + Im 2) T tan -~(Im/Re)

and Re = Real part of z-plane root Irn = Imaginary part of z-plane root. This has been used to calculate the effective zeta of each spectral component in the fitted model. An emperical criterion of significance of (1 - Izl} ~< 0.04 has been used to give automatic selection and printout of important spectral modes. It will be seen in the programme runs section that the spectral components have been ordered in ascending order of frequency and the ratio of significant components calculated to indicate any harmonic relationships. To determine the frequency spectrum of the zpolynomial, the substitution of z = e/t°T is made, where T is the inter-sample period of the data and ¢o is allowed to vary from 0 - 0 . 5 sampling frequency. The magnitude of the spectrum is plotted and it will be seen in the programme runs section that this graph can be closely related to FFT analysis.

D.A. Linkens, P.M. Morrish, Autoregressive modelling o/ biological rythms

3. Programme description

are the sampling period, roll-off of the the filter in dB/oct., lower and upper cutoff frequencies of the filter, the autoregressive model order and the number of data points. These parameters are readout to the plotter and the data plotted after the removal of any

A flow chart of the programmes is shown in fig. 2. The first segment initialises a number of parameters and reads in up to 1024 data points. The parameters

\

freq., o f order,

sampling no. p o i n t s , model parameters

fiIter

~lot ~-plane7

data

[determine

[ residuals &ACF

-

LteUe:°data 1

....

V p1Ot~~ data

]

I

~~esiduals plotACFof/

_ [sdp :cetrmimnereq o'~dei

spectrum

\plot \ filtered//

no

calculate FFT of data Ii model estimate AR coeffs.

I

"

I

[

polynomial~

~

]

testforsignificance1 Fig. 2. Flowchart of complete autoregressive modelling package.

FPlT t /

-]

3

4

D.A. Linkens, P.M. Morrish, Autoregressive modelling o f biological rythms

non-zero mean. If the filter segment is entered, the fdtered data together with its FFT is plotted. The next subroutine calculates the autoregressive model coefficients and is followed by the factorization of the z-polynomial. After this the z-plane diagram with the model roots is plotted together with a printout o f the model coefficients. The frequencies of the roots are calculated, and the parameters of the roots such as distance from the unit circle, equivalent damping ratio are listed in ascending order of spectral frequency. The next segment calculates the residuals and plots meir autocorrelation function to give an indication of the successfulness of the model-fit to the data. The harmonic frequency spectrum of the model is then calculated using the substitution z = e #°T and plotted. Finally a standard FFT is performed on the

i

1.8

I

i

i

data and plotted for comparison with the harmonic spectrum of the fitted model.

4. Typical programme runs The technique was intended to be used for the estimation of frequencies in noisy data where the signal might be present for only a few cycles. The programmes were first tested therefore on simulated data of this type. Figure 3 shows a burst of data at 11 Hz imbedded in white noise with the corresponding zplane plot for an 8th order model which gave very small residuals and an excellent frequency estimation. The programmes were then applied to various typical signals from the mammalian gastro-intestinal

i

i

i

i

@ SIGNAL CONSISTING OF

200 SAMPLE POINTS

t

1.2

t

[

1.0 0.8

'

i

0.4 r

~0.0

ii

m

~0.2 ra~

f

-

-0.6

i

-08 -1.0

~ 20

i 40

i 60

~ 80

, 1 O0

SAMPL~ TIME

Fig. 3a.

, 120

~ 140

p i 60

180

200

D.A. Linkens, P.M. Morrish,Autoregressive modellingof biologicalrythms " - - - ~ - - - T - - T 2.0

r

T

b

q

7

5

~--

Z PLANE

1.5

]NAGIN~RY

1.0

0.5 REAL

0.0

-0.5

-I.0

-I.5

-2.0 -2.0

-1.5

-1.0

-0 5

0.0

0.5

i.0

i.5

2.0

Fig. 3. White noise plus sinusoidal burst. (a) Simulated data; (b) z-plane plot.

tract and the detection technique refmed. Figures 4 a - e show the results obtained from a canine duodenal signal. In this case the data were regular and the sample comprised 200 points with about 30 cycles of information. An 8th order model was adequate to fit this data including the low frequency artefact, as can be seen from the z-plane display in fig. 4b. In this example the fundamental and second harmonic have clearly been estimated accurately, while the real axis roots account for the low frequency drift in the data. Only one additional pair of poles (frequency 0.853) was necessary in this case to account for high frequency noise in the data. The autocorrelation function of the residuals is shown in fig. 4c from which it is seen that there are no significant peaks away from the origin, indicating substantially white noise

residuals and hence a very good model fit. The frequency spectrum corresponding to the model of fig. 4b is plotted in fig. 4d and should be compared with the conventional FFT of the data shown in fig. 4e. The excellent correlation between these results is very clear, particularly bearing in mind that the FFT resolution is only 0.0154 Hz (or ~5%). Increasing the model to 12 gave identical estimates of the fundamental and second harmonic frequencies and very little change in the residuals. Examination of the z-plane diagrams showed that the additional roots were divided between the low frequency and high frequency artefacts in the signals and made practicaUy no difference to the resulting component estimation and frequency spectrum. This gives encouragement to the concept that the best model

CI

Z

SIGNAL

CONS]S

I N G OF

20'0 SAMPLE P O I N T S

~-

c~ ~ o D

-100

a :F

N ~ F- z

-200

-

300 4OO

J

O

10

20 30 40 50 60 SAMPLES PLOTTED AT 0,5000 ]NTERVALS (SECONDS)

B (Z-')'

t,00

A(Z-I) :

1.00 4- O,31Z 4 + O.1OZ"~ - O.13Z "s + O.OZ ~ - O,ooZ~

O, ) 6 Z ~

"

80

90

O'6Z'Z" - O,l°ZS

1.5

FR~S

b

70

Z-P[I, NE

ZETA

O )95E00

9.414L O(

[

/Z/

9.01 >0

REAL

[MAG

,3 r~9OFOO

"3 7£4E00

0,~90t03

~ ~r~3E Ol

0,0~22

0.~7 E09

O ~29E00

I , 0

G S531-33

C ~S3~ 0',

O, I i16

0 79~E09

'3 ~96E9']

O, 5

O, 2951-0,3

S.5IO3

I. ,3,3')0

I s~

G. 590k03

S0~22

2 503~

2~

I,~AGIh~Rv

FR~:G

O. 0

REAL

1

IZl

R~TIO

ffkRMONI C

H_kL AX[S POLLS O, 2749 0.99;6

-0.5 -I.9 - I .5 -1.5

-I.0

-0.5

O,

C

7O

015

i10

1,5

A U T O C O R R E L A T I O N OF THE R E S I D U A L S

OBTAIhED

FROM THE

INVERSE

OF THE

AUTOREGRESSIVE FILLER

6O 5O 40

,~ N

20

<

i0 l 0 F i g . 4 . a, b a n d c.

0

I0

20 30 40 50 60 70 RESIDUALS PLOTTEDAT 0.5000 INTERVALS (SECONDS)

80

90

D.A. Linkens, P.M. Morrish, A utoregressive modelling o f biological rythms

8~' ORDER ?,UTOREGRESS]VE FILTER

d

7

Sm~,]~[CAN~F~EQJ ' ~NCI~S

0

0.2937

0.5900 0.65~2

-S

-/0 ?'

-15

~

-20

~: ~-

-25

d_

50

-

O.

i

h

0.1

0.2

e

i

i

i

~

0.3 0.4 0.5 0.6 0.7 0 8 FREQUENCY (HERTZ.) RESOLUT[ON0,0013 HZ

FASI

FOURIER

TRANSFOE(N

OF SIGNAL

0.9

SIGN!FICANT FRFQUSNCIES O.2~23

0.s~6

120

1o0 80

0

i

O,O

0.1

h

0.2

i

0.3 FREQUENCY

i

0.4 (HERTZ)

0.5

0.6 RESOLUTION

0.7

0.8

0.9

0.0154 dZ

Fig. 4. Canine duodenal recording. (a) Data comprising 200 points; (b) z-plane plot, coefficients, significance indicators; (c) Autocorrelation function of residuals; (d) Frequency spectrum of the AR model; (e) FFT of raw data.

order can be found by iteratively increasing the order until no significant changes occur in the estimates, Alternatively, it appears that the use of an unnecessarily high order of model does not jeopardise the results but only gives an increased computational time. The electrical signals from the canine duodenum

are relatively sinusoidal whereas in the stomach they are very non-sinusoidal. A canine gastric recording comprising about 8 cycles of data and 200 sample points is shown in fig. 5a. The sharp spikes correspond to the basic electrical rhythm while the large pulses occurring after these spikes correspond to the periods of mechanical contraction. A 15th order

a

N

- 200

<~ d

- 400

SIGNAL

CONSISTING

OF

20'0 SAMPLE PO] ,~[i S

/~

600

-

80O

o

io

so SAMPLES

PLOTTED

;o

AT

0.5000

INTERVALS

',,o

;o

9o

B~Z )

I,O0

~,(Z ~]=

~.00 - O. ?OZ ~ * 9,03 -~ * O,08Z -~ . 0,12Z* * O,16Z ~ ~ O,C,~Z"~ - 0,07Z" - O,02Z-~ + O.05Z ° * G,22Z '= * O, I OZ u * 0.94Z ,2 - O,O6Z-" - 0 09Z t~ + 0.30Z ~

I , 5

Fi~E~ b

Z-PLANE IMAGINARY

1. 0

~ X

X

X

~

X

O. 5

0. 0



ZETA

I - /Z/

REAL

IMAC

O.Q43E09

0 267E00

9,0807

O, 753E00

0,512E00

O, 0544

O. 609E00

O, 723E00

0,o555-01

O 0653

O,~3TEO0 O,02J*EO0

O.SSOE-O!

0.!42E09

0.0197

0.190"00

9, 300r O0

O. 277E00

O, 127E00

0.455E00 O,~S6EOO

0.10 EOO

0.0o23

-0, 2~3E89

O. 874E00

O. 720u00

O. 000~-0 ~

O 0Q72

0 576E09

O, 695E00

O, S67EOO

0 742E-01

0,9~63

-0.826509

0,36oE00

~EAL

FRErJ

i - /Z/

RATIg

O,880E-01 O,Olg7

-

8o

(SECONDS~

O.

5

HARMONIC

1. 0000

I~

REAL,AXISPOLES - o. 881,0E00

X

1.5 -I.

-I.0 160

-0.5

0.0

0.5

1.0



C

,

,5 ,

,

AUTOCORQELATIONOF THE RESIDUALS OBTAINED FROM THE INVERSE OF THE A U T O R E G R E S S ] FILTER VE

140

J20 100 80 60 40

o -

20

Fig. 5. a, b and c. i

0

10

20 RESIDUALS

i

30 PLOTTED

AT

40

50

0.5000

INTERVALS

i

60 (SECONDS)

i

70

~

80

90

D.A. Linkens, P.M. Morrish, Autoregressive modelling of biological rythms

9

!

2

d

15TM ORBER AUTOREGRESS!VE FILTER

SlCNIFtC^N~~QUENCreS

0

O.~ 7 5 O, 1900

--~

0.~737

- 4

o. 45'~0 O.q~50

-~

-6

0,7175

_ ~

0. 8662

-12 -14

E g ~ ~

-16 -18

-24 26

-

0.0

0.I

0.2

0.3

0.4

0.5

FREQUENCY (HERTZ)

160

e

0.6 RESOLUTION

0.7 0,0013

0,8

0.9

HZ

FAST FOURIER TRANSFORM OF SICNAL

,o!cNi~ic*mFFctOUt~NC[kS 0 03'3~

I

40

O..3723

/

o.~6~2

0.o615 0,523

o.6.2s

120

/o0

oo

0

, 0.0

0.1

012

0.5 0,4 FREQUENCY (HERTZ)

O.S 0.6 RESOLUTION

0.7 0.9154 HZ

0.8

0.9

Fig. 5. Canine gastric recording. (a) Data comprising 200 points; (b) z-plane plot; (c) Autocorrelation function of residuals; (d) Frequency spectrum of AR model; (e) FFT of raw data.

model gave a good estimate o f the fundamental frequency as can be seen in fig. 5b. Using an empirical criterion o f significance given by {1 - Izl } ~< 0.04 only one significant rhythm was detected at 0.088 Hz. In view o f the large spikes in the data it is not suprising that the autocorrelation function o f the residuals

(fig. 5c) shows that periodic data remained in the residuals. The frequency spectrum o f the model still shows, however, the same features as the F F T of the data, as can be seen by comparing fig. 5d and 5e. The resolution o f the F F T was only about 20% and gave an estimate o f 0.0923 Hz for the fundamental fre-

.

400

N

SIGNAL

O

r

CONSISTING

OF

200

SAMPLE

POINTS

-200

,

/

§e < ~

z

-

400

,'~' -

600

-80O i

0

10

i

20 SAMPLES

30

PLOTTED

40

AT

i

50

0,5000

'[NTERVA, L 5

I. Do

AIZI; ~

/ O0 - 0.66L ~ - 0.27L "~ - O.07Z ~ + 0.08Z ~ ~ O.04Z" . O,03Z" - O.07Z-7 - O.03Z4 . O.OSZ-~ - O.02Z

.

.

.

.

b

FI~O Z-PLANE

• 0

i

70

80

90

(SECONOS)

8 ~Z i)

. ~

i

60

TM

ZETA

+ 0.1SZ 41 - 0 . 0 4 ~ I~

1

-

REAL

/Z/

IMAG

O, ~'X.4E-01

O.567EO0

0,0592

O. 926E00

O. 166E00

0,242E00

O. 434E00

O. 1675

O. 605E00

O. 574E00

OJ,62EO0

0.316E00

0,2152

0,923E-01

0.779E00

O. 627EOC

O. 21 ~ 0 0

O. 1951

-0,314E00

O. 744E00

O. 815E00

O, ~25E00

O. 1 4 8 0

-O.70oEO0

O, 473E00

0,5 FREO

I - IZl

RATIO

HARMONIC

(NONE PRESENT]

O.

0

REAL

REALAXIS POLES O. 2~36

- O.

5

-O.BS~,EOO

,0

-,.

-,.o

-oi~

o.o

C

ols

1.o

.s

AUTOCORRELATION OF THE RESIDUALS OBTAINED FROM THE INVERSE OF THE

AUTOREGRESSIVE FILTER (BROKEN LINES = ±I S.D.I 80

DO

40 o~ L~

z

20

!

-

_

..

~ . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

o r - - - , - - ,

O

i0

. . . .

20 RESIDUALS

,. . . .

.

30

40

50

0.5000

] N T k R V ~ , L c,

PLOTTED

AT

. . . .

,. . . . .

, . . . .

60 [SECONDS)

, . . . .

70



80

. . . .

. . . .

90

D.A. Linkens, P.M. Morrish, Autoregressive modelling of biological rythms

12 ~ff

'

'

.T

'

'

11

-

8 ~

-10

g

N ~

-12

~

N

<

14

-18 2O

0.0

0.1

0.2

0.5

0.4

0.5

FREOUENCY (HERTZ)

e

0.6

0.7

0.8

0.9

RESOLUTION0.0015 HZ

FAST FOURIER TRANSFORM OF SIGNAL

SIGNIFICANTFRE{I~ENCIE$ 0.0462.

140 i20

i00

g

\

80

~ 2

oO

~

40 [email protected] o

0.0

0.1

0.2

0.5

0.4

FREQUENCY (HERTZ)

0.5

0.6

~

0.7

~

0.8

0.9

RESOLUTION0.0154 HZ

Fig. 6. Human colonic recording. (a) Data comprising 200 points; (b) z-plane plot; (c) Autocorrelation function of residuals; (d) Frequency spectrum of AR model; (e) FFT of raw data.

quency. It might be noted that the estimate using the third harmonic value would have been 0.087, compared with the AR model estimate of 0.0875 Hz. A further point of comparison between fig. 5d and 5e is that both FFT and AR model have indicated the presence of a fundamental plus two harmonics with a significant separation between these three components

and the higher frequency noise artefacts. The same type of results were obtained for human gastric recordings which have a square-type waveshape. In these types o f recordings either a high model must be used, or alternatively the raw data can be band-pass filtered using a fairly wide pass-band prior to fitting the AR model. It has been found, for

12

D.A. Linkens, P.M. Morrish, A utoregressive modelling o f biological rythrns

example, that good fits with low order models can be obtained on data such as that of fig. 5a by using a falter which retains the fundamental and 2nd and 3rd harmonics, i.e., an upper cutoff frequency o f about 0.3 Hz. The same remarks apply to circadian rat locomotion data which also should be 'mildly' pre-filtered before ftting. This is not true, however, for human colonic data which is almost sinusoidal in waveshape. In this case it has been possible to give good estimates of dual rhythms when they occur simultaneously in the data. A particularly difficult stretch of data is shown in fig. 6 in which a rhythm can just be detected visually especially at the end of the recording. In spite of the poor data, an autoregressive fit was obtained although most of the z-plane roots were well within the unit circle (fig. 6b). The residuals, however, were substantially white noise (fig. 6c). The model frequency spectrum (fig. 6d) indicated a clear component at 0.051 Hz representing the best fit to the rhythm which visually occurs at about 0.05 Hz. The FFT (fig. 6e) gave an estimate of 0.046 Hz but with a poor resolution of 0.015 Hz.

5. Hardware and software specification The programmes are written in Fortran IV for use on an ICL 1906S machine running under GEORGE IV. The data have been digitised onto paper tape using an ADC attached to a GE4020 process control computer. All the gastro-intestinal signals are stored in

analogue format using 4 channel FM tape recorders. Graph plotting uses the GHOST package and matrix inversion and root factorisation are performed using standard NAG library routines.

6. Mode of availability Source programme listings including explanatory comments are available. All programmes may be available as a card deck on payment of a handling fee.

References [ll D.A. Linkens and A.E. CanneU, IEEE Trans. Biomed. Eng. BME-2, 4 (1974) 335. [2] I. Taylor, H.L. Duthie, R.H. Smallwood and D.A. Linkens, Gut 16 (1975) 808. [3] G.M. Jenkins and D.G. Watts, Spectral analysis and its applications (Holden-Day, London, 1968). [41 F.H. Lopes da Silva, A. Duquesnoy, K. van Hulten and J.G. Lommen: Quantitative Analysis of the EEG, p. 421 (Symp. Jongley sur Vevey, May 1975). [5] S.P. Datardina and D.A. Linkens, in: Random Signal Analysis, p. 82 (lEE Colloqu., London, April 1977). [61 D.A. Linkens and S.P. Datardina, IEEE C,~llf. E!ectr,.. tech. 4.2.8.1. (Venice, May 1977). [7] D.J. Batey, M.J.H. Sterling, D.J. Antcliffe and S.A. Billings, Comput. Aid. Design 7 (1975) 265. [8] R.H. Smallwood, B.H. Brown and H.L. Duthie, 5th Int. Symp. Gastrointest. Motility p. 248 (Leuven, Sept. 1975).