Computer Programs in Biomedicine 9 (1979) 112 © Elsevier/NorthHolland 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 spectraltype analysis for rhythmic biomedical data is presented based on autoregressive modelling techniques. The emphasis has been on portability between dataprocessing 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 prefiltering stages, model coefficients, the polezero 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 nonstationary 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 selfcontained, portable package capable of giving a high precision 'standard' (using an ICL 1906S machine) but capable of being transferred to smaller wordlength 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 leastsquares error. The resulting Filter F(z) is then factorised and the roots plotted on a 'zplane 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 'slowwaves') in the gastrointestinal tract. Spectral analysis of these rhythms has been performed interactively using FFT methods for some years [ 1]. Although these are successful for gastric and smallintestinal slowwaves which are regular and present at all times, FFT analysis is not suited for largeintestinal 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 nonwhite noise. Noise in medical signals is usually nonwhite, 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 leastsquares 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 zpolynomial factorised, the zplane 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 (1Izl}. 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 zplane and are given by the relationship:
Fig. lb. X
Fig. 1. Autoregressivemodellingfor spectral estimation. (a) Schematic diagram; (b)Typical zplane plot.
~" ~/1 + x 2 where
quadrant of the zplane diagram usually indicate significant rhythms while roots in the second and third quadrant represent the 'coloured' spectrum of the nonwhite 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 nonzero mean, prefiltering 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 zplane root Irn = Imaginary part of zplane 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 intersample 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, rolloff 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
nonzero 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 zpolynomial. After this the zplane 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 modelfit 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 gastrointestinal
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) zplane 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 zplane 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 zplane 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(ZI) :
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
ZP[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 S53133
C ~S3~ 0',
O, I i16
0 79~E09
'3 ~96E9']
O, 5
O, 29510,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) zplane 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 nonsinusoidal. 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
ZPLANE 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,o55501
O 0653
O,~3TEO0 O,02J*EO0
O.SSOEO!
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 742E01
0,9~63
0.826509
0,36oE00
~EAL
FRErJ
i  /Z/
RATIg
O,880E01 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) zplane 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.07Z7  O.03Z4 . O.OSZ~  O.02Z
.
.
.
.
b
FI~O ZPLANE
• 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.4E01
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,923E01
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) zplane 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 squaretype waveshape. In these types o f recordings either a high model must be used, or alternatively the raw data can be bandpass filtered using a fairly wide passband 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' prefiltered 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 zplane 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 gastrointestinal 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. BME2, 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 (HoldenDay, 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).