Med. & Biol. Eng. & Comput., 1977, 15, 656-665
Autocorrelation techniques for measuring avian heart rates J. S. B. C l a r k
J.H.
Filshie
Agricultural Research Council Poultry Research Centre, King's Buildings, West Mains Road, Edinburgh EH9 3JS, Scotland
A b s t r a c t - - O n l i n e computer methods have been developed to derive the heart rate of a bird horn its telemetered electrocardiogram (e.c.g.). A PDP-8/ L minicomputer, with real time clock and simple xy-display, is o4 ,=d for this purpose. Autocorrelation methods have been developed to deal with signals which are obscured by the electromyogrem and electrical interference. These techniques produce reliable results in circumstances where more conventional methods are not dependable.
1 Introduction HEART rate is an important index of biological function. It is particularly useful in helping to define the behavioural state of a bird in relation to its environment. In the course of ethological experiments on the measurement of stress it was required to measure heart rate with a 2~o accuracy on the active, unrestrained bird. It was necessary to average the heart rate over the range of 3 s periods at 5 s intervals to 45 s periods at 300 s intervals, each experiment lasting between 45 min and 24 h. In one experiment, birds were subjected individually to stimulation such as the inflation and deflation of a balloon, and heart-rate measurements were made as the birds responded with, in many cases, alarm calls and escape manoeuvres. In addition to the ethological experiments, studies on chick embryo heart rates in the last week of ncubation (LhroaI-IZIN et al., 1975) have suggested that information on the heart rate at an earlier stage would lead to a better understanding of some hatching problems. The chick-embryo electrocardiogram (e.c.g.) is almost fully obscured by noise prior to the 12th day of incubation and can only be determined after signal processing. It is possible to determine heart rate either mechanically, by means of sound or pressure transducers, or electrically, from the e.c.g. Either method is satisfactory when the bird is quiescent. However, such measurements are of interest only in establishing the basal heart rate and the measurements of real interest are those derived from the unrestrained and active bird. I n such situations, intravascular techniques for measuring heart sounds are undesirable and externally placed microphones are both vulnerable and subject to much interference due to skin movement. First received 11th November, 1976 and in final form 7th January 1977
656
The telemetered e.c.g, presents a more reliable indication of heart rate. Nevertheless it is still prone to interference from the electromyogram (e.m.g.) during vigorous activity and from radio-interference in electrically noisy situations. Although the former can be reduced by careful placement of the electrodes and the latter by shielding the system, it is never possible to eliminate interference entirely. It is only under the most noise-free conditions that a simple heart-rate meter can be used. If there is only a small amount of interference on the signal, it is possible to precondition the signal with an automatic-gain-controlled amplifier and Schmitt trigger. This method is of little use if the interference consists of noise spikes of the same or greater amplitude than the QRS complex of the e.c.g. Another possible method of rate measurement is peak counting with removal of erroneous peaks and replacement of missing ones. This is a patternrecognition technique. Its implementation is more suitable for visual work and, moreover, it is capable of extracting only the dominant pattern. Signal averaging can be used to extract the e.c.g. pattern from noise provided that a synchronising signal is present to align the corresponding parts of the e.c.g. This requires a low noise level to guarantee recognition of a unique part of the pattern for synchronisation. It is inherently not a rate-measuring method as the rate is a derivative of the reconstructed signal. Thus the method is either retrospective or demands some foreknowledge of the expected rates. Much of this processing will be redundant since the rate, not the complete waveform, is required. The e.c.g, presents analytical difficulties since it is a quasistationary signal which, although repetitive, exhibits both slow and rapid rate changes. Moreover, interfering signals, such as e.m.g, and electrical noise, are typically in the form of bursts of largeamplitude spikes. Also electrical mains interference can be very troublesome in avian e.c.g, work, more
Medical & Biological Engineering & Computing
November 1977
especially since the half-period of 10 ms coincides approximately with the period of the QRS complex of the waveform. Triggering techniques can be severely restricted by these forms of interference, since they depend on some dominant part of the signals. Fourier analysis and autocorrelation can both be used to extract rate information and have a significant advantage over methods employing triggering, since they deal with the whole waveform. While both these transforms contain frequency information, the Fourier transform also retains phase data; so Fourier analysis requires more processor time than autocorrelation. Since the phase information is of little interest in this application, autocorrelation is to be preferred. 1.1 COmputers and autocorrelation The choice of the digital word length employed depends upon considerations such as probable signal amplitude, signal-to-noise ratio, processor word length and execution time. This can be particularly critical in the type of real-time application considered here. If the processor can perform bit manipulation there are advantages in using 1 bit digitisation, often termed polarity quantisation (JESPERS, 1962; VELTMAN and Bos, 1963). This method saves both storage space and manipulation time since each sample of the signal is represented by a single bit. The performance of both 5 bit and 1 bit autocorrelation was observed. The availability of microprocessors at relatively low prices encouraged the investigation of autocorrelation. It is anticipated that a microprocessor could be dedicated to performing autocorrelation and feeding back rate data to the minicomputer. The 1 bit autocorrelation process could be optimised also, using complete large-scale integrated (1.s.i.) circuit chips designed to carry out this calculation (JORDAN and KELLY, 1976) and again relegating the processor to subsequent handling tasks. FAVRET and CAPUTO (1966) report on the use of autocorrelation to extract foetal e.c.g.s from signals obscured by noise and relatively large amplitude maternal e.c.g.s. They find that certain characteristics in the shape of the foetal e.c.g, lead to destructive interference of the autocorrelogram peak. However, since samples of the avian e.c.g, commonly seen show little trace of these characteristics a n d there is no maternal signal to be eliminated, it seemed worthwhile to investigate the performance of autocorrelation in the case of the avian e.c.g. 2 Apparatus During initial experiments to determine rate by autocorrelation technique, it was found difficult to estimate the variance of the rate. The input signal was therefore idealised to enable this process to be Medical & Biological Engineering & Computing
better understood. The experimental configuration for comparison of multi- and single-bit correlation was a Hewlett-Packard 3722A pseudorandom noise generator, a Hewlett-Packard sweep-frequency oscillator and a P D P - 8 / L processor. A simple laboratory-built video display was used to monitor results in real time. Data was printed out on an A S R 33 Teletype or punched on a Facit 4070 high-speed paper-tape punch for subsequent analysis. Once the response of the autocorrelation process was understood with idealised experimental signals, observations were made under field conditions. The apparatus further required was a Farnell f.m. analogue data tape recorder for convenient playback of the radio-telemetered avian e.c.g., a laboratorybuilt first-order bandpass passive filter, with a bandwidth of 2-80 Hz, and a 50 Hz operational amplifier notch filter for use on very low-amplitude embryo e.c.g.s. 3 Autocorrelation Mathematically, the autocorrelation function of a process x(t) is defined as
j. x(t)x(t+r)dt T
Rxx(r) = lim ~ ar~ T
(1)
0
as in BENDAX and PIERSOL (1966a), where Rx~(r) is the autocorrelation between x(t) and x(t+z). In a digital system, this average can be approximated by sampling the signal every At seconds and summing a finite number of N of the sample products. Division by N normalises the average: 1
N-1
nxx(r) = ~ - Y~
x(Kat)x(gat+o
(2)
K=O
The autocorrelation function can be obtained by successive calculations of the above equation using increasing values of r. The delay increment A~ can be made equal to the sampling interval At. This provides a useful simplification. A program to carry out eqn. 2 in real time will have to perform N multiplications, divisions and additions in every At. Some economies can be made. F o r real-time visual display of the autocorrelation function, relative rather than absolute values of Rxx(z) are required. The normalisation process of division by N may be omitted during display, and performed if required on the coefficients during hard-copy output when the additional time is not important. Since 'exclusive-OR' (XOR) is the logic equivalent of sign multiplication, it can replace multiplication if the incoming data is reduced merely to polarity. Another strategy could be to reduce only one member of the multiplication pair to polarity, while retaining full accuracy for the other. This gives obvious simplification to the multiplication but does November 1977
657
not give savings on storage and buffer manipulation. KORN (1966) gives a relation between 1 bit correlation and the true correlation function Rxx(r) for Gaussian data with zero mean E(sign x(t) sign x ( t - z)} 2
= ~
. arcsm
R~(O
,@[E{xZ(t) X2(t -
./.)}]
where E{fn} is the expected value o f f n . However, the signals dealt with here are not, in general, Gaussian.
4. Program construction 4.1 Multibit autocorrelation program An initial series of questions by the program allows the user to type in the sampling rate R in milliseconds and the number of sampled points N to be held in storage. The time delay over which the autocorrelogram will be displayed is N x R since Az is made equal to At. A memory buffer of length N words is, as a result of the initial dialogue, set up for storage of the incoming sampled data. This acts as a circulating loop buffer with the most recent data word overwriting the oldest, thus sirrulating a shift register. A doubleprecision buffer to hold N accumulating sums is also reserved, and can be displayed on an oscilloscope screen during sampling to form a developing picture of the autocorrelation. The timing routine, which derives its time standard from a crystal oscillator, is initialised to give sampling at the correct intervals. The buffers are cleared. When a fresh signal sample is read, it is stored in the circulating buffer. Crossmultiplications are performed between this 'last sample' and every data word in the circulating buffer, the result of each being cumulatively added into successive positions in the double-precision buffer. This process must be completed and program 'housekeeping' carried out before the next sample is taken if real-time working is to be achieved. The program data collection is started by manual control at the processor teletypewriter at any time after the initial parameters have been entered. Two modes of starting are provided. The immediatestart mode commences with the data sampling and correlation calculations simultaneously. However, this procedure causes a distortion of the autocorrelogram for short averaging times. Since the buffers are initially cleared, the first point (stored in display form) in the double-precision buffer is formed from the sum of S passes, whereas the last point representing maximum delay is the result of ( S - N ) passes, unless S < N, when it will be zero amplitude. The distortion is reduced by increasing S to give a longer averaging time. Following eqn. 1, to make S infinitely large would remove the distortion. Since short averaging times are of interest, a 658
prefill mode of starting is provided. This allows the circulating buffer to fill completely with signal data before the correlation calculations are commenced. Another alternative considered was normalisation by the factor ( S - J) on the Jth point to give a better average of each accumulated sum. This method gave anomalous results for shorter averaging times, and was not found to be as useful as the prefill method. Data collection is halted either after a selected number of samples have been taken, or at any time by operator action. The current pass of calculations is completed to allow the autocorreiogram to be updated. The coefficients can be obtained on the computer printout together with the number of passes and period of averaging. The shortest intersample period found practicable for avian e.c.g, work using the 5 bit program was 8 ms. F r o m spectrum analyses of the e.c.g, waveform, as described later, it was thought that closer sampling might be desirable. To provide this an interrupted-pass mode was included, which halts data taking automatically when the buffer is full, updates calculations meanwhile and continues data taking as soon as possible. This mode results in a much slower autocorrelation integration, but was found useful in verifying some conclusions about sampling rate. 4.2 Single-bit autocorrelation program The polarity-quantised autocorrelation algorithm replaces the sampled digital word by a simple comparator decision on the polarity of the signal. The result of this decision is stored in a single bit, so that 12 signal samples may be stored and manipulated in a single PDP8 word. The multiplication process is replaced by the logic operaor XOR. Although the PDP8 must execute this operation in software, it can be performed considerably faster than multiplication. The X O R process is carried out on complete computer words so that 12 signal samples are crossmultiplied with the latest sample simultaneously, thus saving further time. Maintenance of the circulating buffer is aided by the logical shift instruction which operates on the complete computer word. Since the increments in the cumulative addition process are restricted to a single bit, a double-precision buffer is not necessary unless long integration times are required. The saving in processor time on the overhead of maintaining the double-precision buffer is considerable and the single-precision buffer has been found to be sufficient for dealing with avian e.c.g.s in practice. 4.3 Buffer construction The 'circular' buffer required for incoming data should be readily implemented in principle by maintaining a pointer to indicate the oldest sample in the store and overwriting this with the newest sample. However, the limited instruction repertoire and simple architecture of the PDP-8 require a
Medical & Biological Engineering & Computing
November 1977
heavy overhead to maintain the pointer and to form the wrap-around, converting the linear memory stack to a circular operation. Since the access required is systematic and consecutive, it is convenient and economical to shift all the data down the buffer when creating space for a new sample. When each location is addressed for crossmultiplication, it is also shifted by one place down the buffer. The oldest sample of data at the bottom of the buffer is discarded. The 1-bit version uses a similar strategy, although
it can be implemented more compactly. The link bit attached to the accumulator can be used as an overflow and carry-over register. The updating of the buffer is 2-stage. The first computer word of the buffer containing 12 samples is rotated one place left. The new data sample contained in the link is automatically loaded into the accumulator while the overflow bit is left in the link. The XOR operation can now be performed on this word. The pointer is shifted to the next complete buffer word and a similar rotation performed. This process is
DOCRLN [enter subroutine with latest bit in (AC)
l set (L)
1
and TMPDAT as the latest bit
I
finitia[ise MWDCNT=-number of words in bit store Iinitia[ise PNTRT=start of bit store BBUFST |initiGlise PNTCRL=stQrt of summing store SMBFST
I
r
-~--WOSTEp/get current data word (pointed at by PNTRT I.
1
rotate 1 bit left, replace in bit s(ore using PNTRT
I initiotise bit counter MBTCNT=-12
[ save (L) in
I
LINK 2
I
I
' XNOR TMPDAT with current data word and place resu{t in accumulator
-BTSTEP[ rotate,,, a bit into the link
~
YES
-~
increment current word in summing store pointed to by PNTCRL
.
>its
I reptace LINK 2 into link
I
NO+
I increment current data word pointer PNTRT
--~.
(using MWDCNT)
Medical 8t Biological Engineering & Computing
November 1977
Fig. I Subroutine to perform 1 bit autocorrelation on a P D P - 8 / L
659
repeated until the complete buffer is shifted. The operation of the singl~-bit correlation algorithm is shown in Fig. l, which describes the flow diagram of the appropriate subroutine DOCRLN. 4.4. Sampling rate The choice of sampling rate in any digital data system depends principally on the power spectrum of the signal involved. The Nyquist criterion (BENDAT and PIERSOL, 1966b) requires a sampling frequency at least twice the value of the highest significant frequency present in the power spectrum.
indictate the regions within which the noise envelope was measured. The sampling rate was the same for both cases, with 115 sampled data in a pseudorandom sequence length. The length of signal sample used, corresponding to the integration time, was three sequence lengths in both cases.
100
fS f f oo
g E
2o1 0-0
25
50 frequency, Hz
75
100
Fig. 2 Power spectrum and cumulative power curve of an averaged avian e.c.g.
Power spectra were taken of both individual and averaged e.c.g, samples. They were remarkably consistent. A typical spectrum is shown in Fig. 2. The cumulative power curve shows that only 7 0 ~ of the spectrum is contained in the range up to 50 Hz. This would indicate a sampling frequency of greater than 100 Hz. A regular harmonic structure in the spectrum appears to persist up to 80 Hz, with an accumulated power content of over 8 0 ~ and a Nyquist frequency of 160 Hz. The QRS complex of the avian e.c.g, is approximately l0 ms in duration. An intersample period of less than 5 ms is required to ensure adequate representation of this dominant part of the e.c.g.
5 Investigation of autocorrelation using pseudorandom noise 5.1 Resolving power The pseudorandom noise generator and a continuous band-limited white-noise source were used with an operational-amplifier mixer to provide an input signal which could be set to different signalto-noise ratios. The resolving power of the autocorrelation process was measured as the ratio of the correlation peak lying at one pseudorandom sequence length from the origin, to the average width of the noise envelope lying on the baseline. Fig. 3 shows typical autocorrelograms for 5 and 1 bit processes, respectively. The cursor lines 660
Fig. 3 Autocorrelation of a pseudorandom sequence mixed with Gaussian noise (a) 5 bit process (b) I bit process
The relationship of resolution to integration time is shown in Fig. 4 for signal-to-noise ratios of one and two. The results for 5 and 1 bit autocorrelation sampled at 115 points per signal sequence length and for the 1 bit process sampled at twice this rate are shown in each case. Several curve-fitting programs were used on the data obtained, and the best fits, all close to inverse-square curves, are shown. The r.m.s, error in each case was less than 0.15. The curves show that 1 bit autocorrelation can approach the same resolving power as the multibit process by using its potential for a higher sampling rate, provided that the input signal has the appropriately high frequency content to justify this increase. The harmonic structure of the avian e.c.g. appears on the Nyquist criterion to justify intersample periods of less than 6 ms while the 5 and 1 bit autocorrelation can support a sampling rate equivalent to 8 and 2 ms, respectively, as implemented in the present programs. Thus it is possible
Medical & Biological Engineering & Computing
November 1977
that the l bit process may not suffer a severe penalty of lowered e.c.g, resolution in discarding amplitude information if the sampling rate is increased.
widths remain reasonably constant, reflecting the bandwidth of the signal. Intuitively, the width of
5.2 Response to change of rate The response of autocorrelation to change in the repetition rate of the input-signal pattern was evaluated, using a swept-frequency oscillator as an external clock to change the bit rate of the pseudorandom noise generator linearly. Both 5 and 1 bit autocorrelations exhibit the same type of pattern as in Fig. 5, which shows the autocorrelogram of a pseudorandom sequence which has changed its base frequency by 7 ~ over the averaging period. The successive peak heights diminish but the peak ~2{ Y=5 84x
-o
c
P
y=,~ 92x0 56
Fig, 5 Five bit autocorrelation of a pseudorandom sequence with a linear increase in basic frequency of 7% over the integration period
o
~c
Q
.E
0
1
2 3 4 5 6 7 integration time, signat cydes
the peak can be seen to depend on the bandwidth or mean zero-crossing rate of the signal, since this determines the frequency change or timeshift required to destroy the signal similarity. A plot of fractional change in correlation peak height against integration time, during linear increase of the clock frequency, is shown in Fig. 6. The third parameter is bandwidth. Variation of this is achieved by use of different lowpass filters. The dependence of correlation peak height on frequency drift and on signal bandwidth is clearly seen,
_l 8
a m
30 Y=910x 055 c~
o 2~
L~
6 Autocorrelation on avian e.c.g.
A noise-free sample of e.c.g, was mixed with random Gaussian noise to a QRS peak/r.m.s, value of 2/1 and varying sample lengths autocorrelated. The results are shown in Fig. 7. The peak height
~3
~ 2c c
0"20,
Ck
OIC
~
0 I0 ~
0
1
2
3
4
5
6
7
/ ~
I
fO~q
8
intergration time, signa[ cycte 0
Fig. 4 Resolving power of 5 and I bit autocorrelation (a) S I N = 1 (b) S/N = 2 Curve P 5 bit at 115 points per sequence length Curve Q 1 bit at 230 points per sequence length Curve R I bit at 115 points per sequence length
Medical & Biological Engineering & Computing
500
1000 1500 2000 number of sweeps
2500 3000
Fig. 6 Fractional change in correlation peak height ~R/Ro against number of passes N s for a pseudorandom sequence linearly increasing in frequency, passed through a Iowpass filter with cutoff frequency fo
November 1977
661=
and slope of the sides increased with longer integration time, as expected, although quantitative assessment of this was difficult for the weaker autocorrelograms. For very long integration periods, greater than about 30 s for a quiescent bird, the pattern began to lose resolution due to frequency drift of the e.c.g. In this case, a decrease in the amplitude of the peak indicated the loss of resolution. However, the background noise on the baseline continued to diminish also. The 1 bit autocorrelation process yielded similar results, while requiring an integration about four times longer than the 5 bit process to give similar resolution. For data with lower signal-to-noise
Fig, 7 Autocorrelogrnms of e.c.g, mixed with random noise showing the effect on signal detection of increasing integration time (a) low-noise e.c.g. ; 2 s integration (b) e.c.g, and random noise, ratio 2 . 1 , 3 s integrat/on (c) e,c.g, and random noise, ratio 2 , I ; l o s integration
6_62
ratio, the 1 bit process was less effective, requiring proportionately longer signal samples than the multibit method. Fig. 8 shows 5 bit autocorrelograms of a relatively clean e.c.g, which, in b and c, is mixed with Gaussian noise in the peak-amplitude ratios of 3:1 and 3:2, respectively. Note that the first point on the correlogram in both Figs. 8b and c has been displaced by display buffer overflow to provide a reasonable scale for the e.c.g, period peak. As expected, the response is a reasonably linear reduction in the sharpness of the peak, with regard to both height and width. The performance of the autocorrelation process on an e.c.g, signal undergoing a marked change in rate is shown in Fig. 9. Fig. 9a, taken over the period of change from 280 to 405 beats/rain, shows the
Fig. 8 Five bit autocorre/ograms, integrated over 3s of e.c.g mixed with random noise (a) low-noise e.c g, (b) e.c.g and zancJorn noise, raUo 3,7 (c) e.c.g, and random noise, ratio ,2 2
Medical & Biological Engineering & Computing
November t977
extreme flattening and broadening of the peak. Because the e.c.g, signal remained of characteristic shape, obscured only briefly by some muscle noise, the peak could be distinguished readily from any background noise. The flat top of the correlogram also indicated the linear change of the e.c.g, rate. Figs. 9b and c show the stabilised patterns corresponding to the rates immediately before and after the change of rate. These contrast with Fig. 9a and show that the discriminating power of the correlation process is considerably reduced during the rate transition. In practice, such gross changes are encountered rarely, and the autocorrelograms remain usable in most conditions. 1 and 5 bit autocorrelations respond similarly to changes of e.c.g, rate. The importance of sampling rate is indicated by
the series of autocorrelograms in Fig. 10. These were formed from a 5 bit autocorrelation of a very noisy embryo e.c.g, just barely recognisable as a periodic signal. The deterioration in peak resolution moving from 2 ms sampling rate to 7 ms is evident. F o r these, the interrupted-pass mode had to be used. Fig. 10c, with a sampling rate of 8 ms, was formed in the continuous sampling mode and shows an improvement in peak resolution compared to the intermediate rates, Fig. 10b. It seems likely that the interrupted-pass mode, with its loss of continuity of the sampled signal, and, therefore, loss of some correlation information, is considerably less efficient than the other mode. However, the higher sampling rates outweigh this disadvantage and result in a nett improvement in resolution. Some evidence of a triphasic QRS complex was
Fig, 9 Five bit autocorrelograrn o/ e.c.g, during a rate change (a) during rate inclease from 280 to 405 beats/rnm (b) before rate change (c) temporarily stabi/ised after rate increase
Fig 10 Five b/t autocorre/ograms of noisy e c g. h o m a 12 day embryo chick integrated over 20 s (a) sampling rate 2 ms /b) sampling fate 7 his (c) sampl/ng rate 8 ,Ti,s
Medical & Biological Engineering & C o m p u t i n g
N o v e m b e r 1977
663
noted in a part of the chick embryo e.c.g, used to form these correlograms. Fig 10a shows the effect of this as a negative-correlation peak which occurs immediately before the e.c.g, period peak, as mentioned by FAWET and CAPUTO (1966). However, the phenomenon was never important enough to cause significant cancellation of the correlation peak as in their work on foetal e.c.g.s. Straightforward facilities to automate the rate measurement for data-logging purposes have been devised. A maximum-seeking subroutine is used to search the correlogram buffer between preset limits for the largest amplitude maximum and interpolate the peak position using the positive and negative slopes. It is found useful to provide flexible limits to the searching process for rapidly varying e.c.g.s. The limits are modified according to trends established from the previous two results. This allows the limits to be set more closely together, providing a narrower window which gives greater immunity to erroneous peaks as well as saving execution time. The x-ordinate position of the correlation peak is multiplied by the intersample period to give the average e.c.g, period. This can be typed or punched out immediately, or retained in a data buffer. From the buffer, the rate data can be subsequently displayed to identify trends visually and subjected to other statistical operations at the end of data collection.
activity, such as teeding, when e.m.g, and other disturbing signals were present. 7 Conclusion Autocorrelation was found a very useful technique in detecting and logging heart rates in both relatively noise-free and noisy situations. The 1 bit process is useful for less noisy signals where sporadic bursts of large-amplitude noise may occur, since it is less affected by these than the 5 bit process, yet can be executed faster and achieve higher sampling rates. The 5 bit process performs well on generally noisy signals and those which cannot be visually detected as periodic. B o t h processes give promise of good potential in this field using faster proglams implemented on more sophisticated minicomputers, dedicated microprocessors and custom 1.s.i. chips now available.
450
•
E400
~
rT
350 L)
~soo 250
j
o
~200 0
4
8
F
I
12 16 20 24 28 32 36 40 44 48
time, rain Fig. 11 Visual stimulation of this broiler occurred at each vertical cursor. Heart-rate readings were derived from 3 s autocorrelation processes taken at 5 s intervals
Programs using the 5 and 1 bit algorithms and the above strategies have been used successfully for the analyses of a large number of e.c.g, records in several ethological experiments. 3 s averages, taken continuously or at 5 s intervals, were required. The birds were subjected to periodic auditory or visual stimulation, such as the inflating of a balloon. Fig. 11 shows a plot of the output of one such experiment, with the stimulation events m a r k e d by vertical cursor lines. Accurate and consistent logging was achieved through these periods of experimental stimulation of the poultry and during normal 664
References BENDAT,J. S. and PIERSOL,A. G. (1966a) Measurement and analysis o f random data. Wiley, 19. BENDAT, J. S. and PmRSOL, A. G. (1966b) Measurement and analysis o f random data. Wiley, 279. FAVRET, A. G. and CAPOTO, A. F. (1966) Evaluation of autocorrelation techniques for detection of the fetal electrocardiogram. IEEE Trans., BME--13, 37-43. JESPERS, P. (1962) .4 new method to compuw correlation Junctions. Proceedings of the IEEE Symposium on Information Theory: Brussels. JOROAN, J. R. and KELLY,R. G. (1976) Integrated circuit correlator for flow measurement. Meas. & Control 9, 267-270. KORN, G. A. 0966) Random process simulation and measurement. McGraw-Hill, 135. I~AUGHLtN,K. F., LUNDY,H. and TAIT, J. A. (1975) A method for monitoring embryonic heart rate during the last week of incubation. J. Physiol. 244, 8P-9P. VELTMAN, B. P. Th. and VAN DEN BOS, (1963) On the applicability o f the relay-correlator and the polaritycoincidence correlator in automatic control. Proceedings
of the Second International Congress IFAC, Basel.
Medical & Biological Engineering & Computing
November 1977
Traitement direct des donn6es de signaux electrocardiographiques aviens et les techniques utilis6es pour I'obtention des signaux par conditions de bruit dlev6 Somraaire---Certaines m6thodes de traitement direct par ordinateur ont et6 mises au point pour obtenir le r6gime cardiaque d'un oiseau ~. partir de son 61ectrocadiogramme t616mesur6 (e.c.g.). A cet effet, on utilise un mini-calculateur P D P - 8 / L avec horloge en temps r6el et un simple affichage du type-xy. Des m6thodes d'auto-corr61ation ont 6galement 6t6 d6velopp6es pour le traitement des signaux susceptibles d'6tre obscurcis par l'61ectromyographe ou les interf&ences 61ectriques. Ces techniques ont pour effet de produire des r6sultats fiables sous certaines conditions ou les m6thodes plus conventionnelles ne le sont pas.
Rechnerabhi~ngige Verarbeitung der ekg-Signale von V6geln und Verfahren zur Ableitung dieser Signale bei starkem Rauschen Zusammenfassung--Rechnergesteuerte Computer-Verfahren wurden entwickelt, um den Herzschlag von V6geln aus FermeB-EKG-Daten abzuteiten. Ein Minicomputer vom Typ P D P - 8 / L mit Realzeituhr und einfacherxy-Anzeige wird zu diesem Zweck verwendet. Es wurden automatische Korrelationsverfahren entwickelt, die Signale verarbeiten, die durch Elektromyogramm- und Stromst6rungen iiberlagert werden. Diese Verfahren zeigen zuverl~issige Ergebnisse unter Umst~inden, wo konventionellere Verfahren unzuverlfissig sind.
Medical & Biological Engineering & Computing
November 1977
665
Autocorrelation techniques for measuring avian heart rates J. S. B. C l a r k
J.H.
Filshie
Agricultural Research Council Poultry Research Centre, King's Buildings, West Mains Road, Edinburgh EH9 3JS, Scotland
A b s t r a c t - - O n l i n e computer methods have been developed to derive the heart rate of a bird horn its telemetered electrocardiogram (e.c.g.). A PDP-8/ L minicomputer, with real time clock and simple xy-display, is o4 ,=d for this purpose. Autocorrelation methods have been developed to deal with signals which are obscured by the electromyogrem and electrical interference. These techniques produce reliable results in circumstances where more conventional methods are not dependable.
1 Introduction HEART rate is an important index of biological function. It is particularly useful in helping to define the behavioural state of a bird in relation to its environment. In the course of ethological experiments on the measurement of stress it was required to measure heart rate with a 2~o accuracy on the active, unrestrained bird. It was necessary to average the heart rate over the range of 3 s periods at 5 s intervals to 45 s periods at 300 s intervals, each experiment lasting between 45 min and 24 h. In one experiment, birds were subjected individually to stimulation such as the inflation and deflation of a balloon, and heart-rate measurements were made as the birds responded with, in many cases, alarm calls and escape manoeuvres. In addition to the ethological experiments, studies on chick embryo heart rates in the last week of ncubation (LhroaI-IZIN et al., 1975) have suggested that information on the heart rate at an earlier stage would lead to a better understanding of some hatching problems. The chick-embryo electrocardiogram (e.c.g.) is almost fully obscured by noise prior to the 12th day of incubation and can only be determined after signal processing. It is possible to determine heart rate either mechanically, by means of sound or pressure transducers, or electrically, from the e.c.g. Either method is satisfactory when the bird is quiescent. However, such measurements are of interest only in establishing the basal heart rate and the measurements of real interest are those derived from the unrestrained and active bird. I n such situations, intravascular techniques for measuring heart sounds are undesirable and externally placed microphones are both vulnerable and subject to much interference due to skin movement. First received 11th November, 1976 and in final form 7th January 1977
656
The telemetered e.c.g, presents a more reliable indication of heart rate. Nevertheless it is still prone to interference from the electromyogram (e.m.g.) during vigorous activity and from radio-interference in electrically noisy situations. Although the former can be reduced by careful placement of the electrodes and the latter by shielding the system, it is never possible to eliminate interference entirely. It is only under the most noise-free conditions that a simple heart-rate meter can be used. If there is only a small amount of interference on the signal, it is possible to precondition the signal with an automatic-gain-controlled amplifier and Schmitt trigger. This method is of little use if the interference consists of noise spikes of the same or greater amplitude than the QRS complex of the e.c.g. Another possible method of rate measurement is peak counting with removal of erroneous peaks and replacement of missing ones. This is a patternrecognition technique. Its implementation is more suitable for visual work and, moreover, it is capable of extracting only the dominant pattern. Signal averaging can be used to extract the e.c.g. pattern from noise provided that a synchronising signal is present to align the corresponding parts of the e.c.g. This requires a low noise level to guarantee recognition of a unique part of the pattern for synchronisation. It is inherently not a rate-measuring method as the rate is a derivative of the reconstructed signal. Thus the method is either retrospective or demands some foreknowledge of the expected rates. Much of this processing will be redundant since the rate, not the complete waveform, is required. The e.c.g, presents analytical difficulties since it is a quasistationary signal which, although repetitive, exhibits both slow and rapid rate changes. Moreover, interfering signals, such as e.m.g, and electrical noise, are typically in the form of bursts of largeamplitude spikes. Also electrical mains interference can be very troublesome in avian e.c.g, work, more
Medical & Biological Engineering & Computing
November 1977
especially since the half-period of 10 ms coincides approximately with the period of the QRS complex of the waveform. Triggering techniques can be severely restricted by these forms of interference, since they depend on some dominant part of the signals. Fourier analysis and autocorrelation can both be used to extract rate information and have a significant advantage over methods employing triggering, since they deal with the whole waveform. While both these transforms contain frequency information, the Fourier transform also retains phase data; so Fourier analysis requires more processor time than autocorrelation. Since the phase information is of little interest in this application, autocorrelation is to be preferred. 1.1 COmputers and autocorrelation The choice of the digital word length employed depends upon considerations such as probable signal amplitude, signal-to-noise ratio, processor word length and execution time. This can be particularly critical in the type of real-time application considered here. If the processor can perform bit manipulation there are advantages in using 1 bit digitisation, often termed polarity quantisation (JESPERS, 1962; VELTMAN and Bos, 1963). This method saves both storage space and manipulation time since each sample of the signal is represented by a single bit. The performance of both 5 bit and 1 bit autocorrelation was observed. The availability of microprocessors at relatively low prices encouraged the investigation of autocorrelation. It is anticipated that a microprocessor could be dedicated to performing autocorrelation and feeding back rate data to the minicomputer. The 1 bit autocorrelation process could be optimised also, using complete large-scale integrated (1.s.i.) circuit chips designed to carry out this calculation (JORDAN and KELLY, 1976) and again relegating the processor to subsequent handling tasks. FAVRET and CAPUTO (1966) report on the use of autocorrelation to extract foetal e.c.g.s from signals obscured by noise and relatively large amplitude maternal e.c.g.s. They find that certain characteristics in the shape of the foetal e.c.g, lead to destructive interference of the autocorrelogram peak. However, since samples of the avian e.c.g, commonly seen show little trace of these characteristics a n d there is no maternal signal to be eliminated, it seemed worthwhile to investigate the performance of autocorrelation in the case of the avian e.c.g. 2 Apparatus During initial experiments to determine rate by autocorrelation technique, it was found difficult to estimate the variance of the rate. The input signal was therefore idealised to enable this process to be Medical & Biological Engineering & Computing
better understood. The experimental configuration for comparison of multi- and single-bit correlation was a Hewlett-Packard 3722A pseudorandom noise generator, a Hewlett-Packard sweep-frequency oscillator and a P D P - 8 / L processor. A simple laboratory-built video display was used to monitor results in real time. Data was printed out on an A S R 33 Teletype or punched on a Facit 4070 high-speed paper-tape punch for subsequent analysis. Once the response of the autocorrelation process was understood with idealised experimental signals, observations were made under field conditions. The apparatus further required was a Farnell f.m. analogue data tape recorder for convenient playback of the radio-telemetered avian e.c.g., a laboratorybuilt first-order bandpass passive filter, with a bandwidth of 2-80 Hz, and a 50 Hz operational amplifier notch filter for use on very low-amplitude embryo e.c.g.s. 3 Autocorrelation Mathematically, the autocorrelation function of a process x(t) is defined as
j. x(t)x(t+r)dt T
Rxx(r) = lim ~ ar~ T
(1)
0
as in BENDAX and PIERSOL (1966a), where Rx~(r) is the autocorrelation between x(t) and x(t+z). In a digital system, this average can be approximated by sampling the signal every At seconds and summing a finite number of N of the sample products. Division by N normalises the average: 1
N-1
nxx(r) = ~ - Y~
x(Kat)x(gat+o
(2)
K=O
The autocorrelation function can be obtained by successive calculations of the above equation using increasing values of r. The delay increment A~ can be made equal to the sampling interval At. This provides a useful simplification. A program to carry out eqn. 2 in real time will have to perform N multiplications, divisions and additions in every At. Some economies can be made. F o r real-time visual display of the autocorrelation function, relative rather than absolute values of Rxx(z) are required. The normalisation process of division by N may be omitted during display, and performed if required on the coefficients during hard-copy output when the additional time is not important. Since 'exclusive-OR' (XOR) is the logic equivalent of sign multiplication, it can replace multiplication if the incoming data is reduced merely to polarity. Another strategy could be to reduce only one member of the multiplication pair to polarity, while retaining full accuracy for the other. This gives obvious simplification to the multiplication but does November 1977
657
not give savings on storage and buffer manipulation. KORN (1966) gives a relation between 1 bit correlation and the true correlation function Rxx(r) for Gaussian data with zero mean E(sign x(t) sign x ( t - z)} 2
= ~
. arcsm
R~(O
,@[E{xZ(t) X2(t -
./.)}]
where E{fn} is the expected value o f f n . However, the signals dealt with here are not, in general, Gaussian.
4. Program construction 4.1 Multibit autocorrelation program An initial series of questions by the program allows the user to type in the sampling rate R in milliseconds and the number of sampled points N to be held in storage. The time delay over which the autocorrelogram will be displayed is N x R since Az is made equal to At. A memory buffer of length N words is, as a result of the initial dialogue, set up for storage of the incoming sampled data. This acts as a circulating loop buffer with the most recent data word overwriting the oldest, thus sirrulating a shift register. A doubleprecision buffer to hold N accumulating sums is also reserved, and can be displayed on an oscilloscope screen during sampling to form a developing picture of the autocorrelation. The timing routine, which derives its time standard from a crystal oscillator, is initialised to give sampling at the correct intervals. The buffers are cleared. When a fresh signal sample is read, it is stored in the circulating buffer. Crossmultiplications are performed between this 'last sample' and every data word in the circulating buffer, the result of each being cumulatively added into successive positions in the double-precision buffer. This process must be completed and program 'housekeeping' carried out before the next sample is taken if real-time working is to be achieved. The program data collection is started by manual control at the processor teletypewriter at any time after the initial parameters have been entered. Two modes of starting are provided. The immediatestart mode commences with the data sampling and correlation calculations simultaneously. However, this procedure causes a distortion of the autocorrelogram for short averaging times. Since the buffers are initially cleared, the first point (stored in display form) in the double-precision buffer is formed from the sum of S passes, whereas the last point representing maximum delay is the result of ( S - N ) passes, unless S < N, when it will be zero amplitude. The distortion is reduced by increasing S to give a longer averaging time. Following eqn. 1, to make S infinitely large would remove the distortion. Since short averaging times are of interest, a 658
prefill mode of starting is provided. This allows the circulating buffer to fill completely with signal data before the correlation calculations are commenced. Another alternative considered was normalisation by the factor ( S - J) on the Jth point to give a better average of each accumulated sum. This method gave anomalous results for shorter averaging times, and was not found to be as useful as the prefill method. Data collection is halted either after a selected number of samples have been taken, or at any time by operator action. The current pass of calculations is completed to allow the autocorreiogram to be updated. The coefficients can be obtained on the computer printout together with the number of passes and period of averaging. The shortest intersample period found practicable for avian e.c.g, work using the 5 bit program was 8 ms. F r o m spectrum analyses of the e.c.g, waveform, as described later, it was thought that closer sampling might be desirable. To provide this an interrupted-pass mode was included, which halts data taking automatically when the buffer is full, updates calculations meanwhile and continues data taking as soon as possible. This mode results in a much slower autocorrelation integration, but was found useful in verifying some conclusions about sampling rate. 4.2 Single-bit autocorrelation program The polarity-quantised autocorrelation algorithm replaces the sampled digital word by a simple comparator decision on the polarity of the signal. The result of this decision is stored in a single bit, so that 12 signal samples may be stored and manipulated in a single PDP8 word. The multiplication process is replaced by the logic operaor XOR. Although the PDP8 must execute this operation in software, it can be performed considerably faster than multiplication. The X O R process is carried out on complete computer words so that 12 signal samples are crossmultiplied with the latest sample simultaneously, thus saving further time. Maintenance of the circulating buffer is aided by the logical shift instruction which operates on the complete computer word. Since the increments in the cumulative addition process are restricted to a single bit, a double-precision buffer is not necessary unless long integration times are required. The saving in processor time on the overhead of maintaining the double-precision buffer is considerable and the single-precision buffer has been found to be sufficient for dealing with avian e.c.g.s in practice. 4.3 Buffer construction The 'circular' buffer required for incoming data should be readily implemented in principle by maintaining a pointer to indicate the oldest sample in the store and overwriting this with the newest sample. However, the limited instruction repertoire and simple architecture of the PDP-8 require a
Medical & Biological Engineering & Computing
November 1977
heavy overhead to maintain the pointer and to form the wrap-around, converting the linear memory stack to a circular operation. Since the access required is systematic and consecutive, it is convenient and economical to shift all the data down the buffer when creating space for a new sample. When each location is addressed for crossmultiplication, it is also shifted by one place down the buffer. The oldest sample of data at the bottom of the buffer is discarded. The 1-bit version uses a similar strategy, although
it can be implemented more compactly. The link bit attached to the accumulator can be used as an overflow and carry-over register. The updating of the buffer is 2-stage. The first computer word of the buffer containing 12 samples is rotated one place left. The new data sample contained in the link is automatically loaded into the accumulator while the overflow bit is left in the link. The XOR operation can now be performed on this word. The pointer is shifted to the next complete buffer word and a similar rotation performed. This process is
DOCRLN [enter subroutine with latest bit in (AC)
l set (L)
1
and TMPDAT as the latest bit
I
finitia[ise MWDCNT=-number of words in bit store Iinitia[ise PNTRT=start of bit store BBUFST |initiGlise PNTCRL=stQrt of summing store SMBFST
I
r
-~--WOSTEp/get current data word (pointed at by PNTRT I.
1
rotate 1 bit left, replace in bit s(ore using PNTRT
I initiotise bit counter MBTCNT=-12
[ save (L) in
I
LINK 2
I
I
' XNOR TMPDAT with current data word and place resu{t in accumulator
-BTSTEP[ rotate,,, a bit into the link
~
YES
-~
increment current word in summing store pointed to by PNTCRL
.
>its
I reptace LINK 2 into link
I
NO+
I increment current data word pointer PNTRT
--~.
(using MWDCNT)
Medical 8t Biological Engineering & Computing
November 1977
Fig. I Subroutine to perform 1 bit autocorrelation on a P D P - 8 / L
659
repeated until the complete buffer is shifted. The operation of the singl~-bit correlation algorithm is shown in Fig. l, which describes the flow diagram of the appropriate subroutine DOCRLN. 4.4. Sampling rate The choice of sampling rate in any digital data system depends principally on the power spectrum of the signal involved. The Nyquist criterion (BENDAT and PIERSOL, 1966b) requires a sampling frequency at least twice the value of the highest significant frequency present in the power spectrum.
indictate the regions within which the noise envelope was measured. The sampling rate was the same for both cases, with 115 sampled data in a pseudorandom sequence length. The length of signal sample used, corresponding to the integration time, was three sequence lengths in both cases.
100
fS f f oo
g E
2o1 0-0
25
50 frequency, Hz
75
100
Fig. 2 Power spectrum and cumulative power curve of an averaged avian e.c.g.
Power spectra were taken of both individual and averaged e.c.g, samples. They were remarkably consistent. A typical spectrum is shown in Fig. 2. The cumulative power curve shows that only 7 0 ~ of the spectrum is contained in the range up to 50 Hz. This would indicate a sampling frequency of greater than 100 Hz. A regular harmonic structure in the spectrum appears to persist up to 80 Hz, with an accumulated power content of over 8 0 ~ and a Nyquist frequency of 160 Hz. The QRS complex of the avian e.c.g, is approximately l0 ms in duration. An intersample period of less than 5 ms is required to ensure adequate representation of this dominant part of the e.c.g.
5 Investigation of autocorrelation using pseudorandom noise 5.1 Resolving power The pseudorandom noise generator and a continuous band-limited white-noise source were used with an operational-amplifier mixer to provide an input signal which could be set to different signalto-noise ratios. The resolving power of the autocorrelation process was measured as the ratio of the correlation peak lying at one pseudorandom sequence length from the origin, to the average width of the noise envelope lying on the baseline. Fig. 3 shows typical autocorrelograms for 5 and 1 bit processes, respectively. The cursor lines 660
Fig. 3 Autocorrelation of a pseudorandom sequence mixed with Gaussian noise (a) 5 bit process (b) I bit process
The relationship of resolution to integration time is shown in Fig. 4 for signal-to-noise ratios of one and two. The results for 5 and 1 bit autocorrelation sampled at 115 points per signal sequence length and for the 1 bit process sampled at twice this rate are shown in each case. Several curve-fitting programs were used on the data obtained, and the best fits, all close to inverse-square curves, are shown. The r.m.s, error in each case was less than 0.15. The curves show that 1 bit autocorrelation can approach the same resolving power as the multibit process by using its potential for a higher sampling rate, provided that the input signal has the appropriately high frequency content to justify this increase. The harmonic structure of the avian e.c.g. appears on the Nyquist criterion to justify intersample periods of less than 6 ms while the 5 and 1 bit autocorrelation can support a sampling rate equivalent to 8 and 2 ms, respectively, as implemented in the present programs. Thus it is possible
Medical & Biological Engineering & Computing
November 1977
that the l bit process may not suffer a severe penalty of lowered e.c.g, resolution in discarding amplitude information if the sampling rate is increased.
widths remain reasonably constant, reflecting the bandwidth of the signal. Intuitively, the width of
5.2 Response to change of rate The response of autocorrelation to change in the repetition rate of the input-signal pattern was evaluated, using a swept-frequency oscillator as an external clock to change the bit rate of the pseudorandom noise generator linearly. Both 5 and 1 bit autocorrelations exhibit the same type of pattern as in Fig. 5, which shows the autocorrelogram of a pseudorandom sequence which has changed its base frequency by 7 ~ over the averaging period. The successive peak heights diminish but the peak ~2{ Y=5 84x
-o
c
P
y=,~ 92x0 56
Fig, 5 Five bit autocorrelation of a pseudorandom sequence with a linear increase in basic frequency of 7% over the integration period
o
~c
Q
.E
0
1
2 3 4 5 6 7 integration time, signat cydes
the peak can be seen to depend on the bandwidth or mean zero-crossing rate of the signal, since this determines the frequency change or timeshift required to destroy the signal similarity. A plot of fractional change in correlation peak height against integration time, during linear increase of the clock frequency, is shown in Fig. 6. The third parameter is bandwidth. Variation of this is achieved by use of different lowpass filters. The dependence of correlation peak height on frequency drift and on signal bandwidth is clearly seen,
_l 8
a m
30 Y=910x 055 c~
o 2~
L~
6 Autocorrelation on avian e.c.g.
A noise-free sample of e.c.g, was mixed with random Gaussian noise to a QRS peak/r.m.s, value of 2/1 and varying sample lengths autocorrelated. The results are shown in Fig. 7. The peak height
~3
~ 2c c
0"20,
Ck
OIC
~
0 I0 ~
0
1
2
3
4
5
6
7
/ ~
I
fO~q
8
intergration time, signa[ cycte 0
Fig. 4 Resolving power of 5 and I bit autocorrelation (a) S I N = 1 (b) S/N = 2 Curve P 5 bit at 115 points per sequence length Curve Q 1 bit at 230 points per sequence length Curve R I bit at 115 points per sequence length
Medical & Biological Engineering & Computing
500
1000 1500 2000 number of sweeps
2500 3000
Fig. 6 Fractional change in correlation peak height ~R/Ro against number of passes N s for a pseudorandom sequence linearly increasing in frequency, passed through a Iowpass filter with cutoff frequency fo
November 1977
661=
and slope of the sides increased with longer integration time, as expected, although quantitative assessment of this was difficult for the weaker autocorrelograms. For very long integration periods, greater than about 30 s for a quiescent bird, the pattern began to lose resolution due to frequency drift of the e.c.g. In this case, a decrease in the amplitude of the peak indicated the loss of resolution. However, the background noise on the baseline continued to diminish also. The 1 bit autocorrelation process yielded similar results, while requiring an integration about four times longer than the 5 bit process to give similar resolution. For data with lower signal-to-noise
Fig, 7 Autocorrelogrnms of e.c.g, mixed with random noise showing the effect on signal detection of increasing integration time (a) low-noise e.c.g. ; 2 s integration (b) e.c.g, and random noise, ratio 2 . 1 , 3 s integrat/on (c) e,c.g, and random noise, ratio 2 , I ; l o s integration
6_62
ratio, the 1 bit process was less effective, requiring proportionately longer signal samples than the multibit method. Fig. 8 shows 5 bit autocorrelograms of a relatively clean e.c.g, which, in b and c, is mixed with Gaussian noise in the peak-amplitude ratios of 3:1 and 3:2, respectively. Note that the first point on the correlogram in both Figs. 8b and c has been displaced by display buffer overflow to provide a reasonable scale for the e.c.g, period peak. As expected, the response is a reasonably linear reduction in the sharpness of the peak, with regard to both height and width. The performance of the autocorrelation process on an e.c.g, signal undergoing a marked change in rate is shown in Fig. 9. Fig. 9a, taken over the period of change from 280 to 405 beats/rain, shows the
Fig. 8 Five bit autocorre/ograms, integrated over 3s of e.c.g mixed with random noise (a) low-noise e.c g, (b) e.c.g and zancJorn noise, raUo 3,7 (c) e.c.g, and random noise, ratio ,2 2
Medical & Biological Engineering & Computing
November t977
extreme flattening and broadening of the peak. Because the e.c.g, signal remained of characteristic shape, obscured only briefly by some muscle noise, the peak could be distinguished readily from any background noise. The flat top of the correlogram also indicated the linear change of the e.c.g, rate. Figs. 9b and c show the stabilised patterns corresponding to the rates immediately before and after the change of rate. These contrast with Fig. 9a and show that the discriminating power of the correlation process is considerably reduced during the rate transition. In practice, such gross changes are encountered rarely, and the autocorrelograms remain usable in most conditions. 1 and 5 bit autocorrelations respond similarly to changes of e.c.g, rate. The importance of sampling rate is indicated by
the series of autocorrelograms in Fig. 10. These were formed from a 5 bit autocorrelation of a very noisy embryo e.c.g, just barely recognisable as a periodic signal. The deterioration in peak resolution moving from 2 ms sampling rate to 7 ms is evident. F o r these, the interrupted-pass mode had to be used. Fig. 10c, with a sampling rate of 8 ms, was formed in the continuous sampling mode and shows an improvement in peak resolution compared to the intermediate rates, Fig. 10b. It seems likely that the interrupted-pass mode, with its loss of continuity of the sampled signal, and, therefore, loss of some correlation information, is considerably less efficient than the other mode. However, the higher sampling rates outweigh this disadvantage and result in a nett improvement in resolution. Some evidence of a triphasic QRS complex was
Fig, 9 Five bit autocorrelograrn o/ e.c.g, during a rate change (a) during rate inclease from 280 to 405 beats/rnm (b) before rate change (c) temporarily stabi/ised after rate increase
Fig 10 Five b/t autocorre/ograms of noisy e c g. h o m a 12 day embryo chick integrated over 20 s (a) sampling rate 2 ms /b) sampling fate 7 his (c) sampl/ng rate 8 ,Ti,s
Medical & Biological Engineering & C o m p u t i n g
N o v e m b e r 1977
663
noted in a part of the chick embryo e.c.g, used to form these correlograms. Fig 10a shows the effect of this as a negative-correlation peak which occurs immediately before the e.c.g, period peak, as mentioned by FAWET and CAPUTO (1966). However, the phenomenon was never important enough to cause significant cancellation of the correlation peak as in their work on foetal e.c.g.s. Straightforward facilities to automate the rate measurement for data-logging purposes have been devised. A maximum-seeking subroutine is used to search the correlogram buffer between preset limits for the largest amplitude maximum and interpolate the peak position using the positive and negative slopes. It is found useful to provide flexible limits to the searching process for rapidly varying e.c.g.s. The limits are modified according to trends established from the previous two results. This allows the limits to be set more closely together, providing a narrower window which gives greater immunity to erroneous peaks as well as saving execution time. The x-ordinate position of the correlation peak is multiplied by the intersample period to give the average e.c.g, period. This can be typed or punched out immediately, or retained in a data buffer. From the buffer, the rate data can be subsequently displayed to identify trends visually and subjected to other statistical operations at the end of data collection.
activity, such as teeding, when e.m.g, and other disturbing signals were present. 7 Conclusion Autocorrelation was found a very useful technique in detecting and logging heart rates in both relatively noise-free and noisy situations. The 1 bit process is useful for less noisy signals where sporadic bursts of large-amplitude noise may occur, since it is less affected by these than the 5 bit process, yet can be executed faster and achieve higher sampling rates. The 5 bit process performs well on generally noisy signals and those which cannot be visually detected as periodic. B o t h processes give promise of good potential in this field using faster proglams implemented on more sophisticated minicomputers, dedicated microprocessors and custom 1.s.i. chips now available.
450
•
E400
~
rT
350 L)
~soo 250
j
o
~200 0
4
8
F
I
12 16 20 24 28 32 36 40 44 48
time, rain Fig. 11 Visual stimulation of this broiler occurred at each vertical cursor. Heart-rate readings were derived from 3 s autocorrelation processes taken at 5 s intervals
Programs using the 5 and 1 bit algorithms and the above strategies have been used successfully for the analyses of a large number of e.c.g, records in several ethological experiments. 3 s averages, taken continuously or at 5 s intervals, were required. The birds were subjected to periodic auditory or visual stimulation, such as the inflating of a balloon. Fig. 11 shows a plot of the output of one such experiment, with the stimulation events m a r k e d by vertical cursor lines. Accurate and consistent logging was achieved through these periods of experimental stimulation of the poultry and during normal 664
References BENDAT,J. S. and PIERSOL,A. G. (1966a) Measurement and analysis o f random data. Wiley, 19. BENDAT, J. S. and PmRSOL, A. G. (1966b) Measurement and analysis o f random data. Wiley, 279. FAVRET, A. G. and CAPOTO, A. F. (1966) Evaluation of autocorrelation techniques for detection of the fetal electrocardiogram. IEEE Trans., BME--13, 37-43. JESPERS, P. (1962) .4 new method to compuw correlation Junctions. Proceedings of the IEEE Symposium on Information Theory: Brussels. JOROAN, J. R. and KELLY,R. G. (1976) Integrated circuit correlator for flow measurement. Meas. & Control 9, 267-270. KORN, G. A. 0966) Random process simulation and measurement. McGraw-Hill, 135. I~AUGHLtN,K. F., LUNDY,H. and TAIT, J. A. (1975) A method for monitoring embryonic heart rate during the last week of incubation. J. Physiol. 244, 8P-9P. VELTMAN, B. P. Th. and VAN DEN BOS, (1963) On the applicability o f the relay-correlator and the polaritycoincidence correlator in automatic control. Proceedings
of the Second International Congress IFAC, Basel.
Medical & Biological Engineering & Computing
November 1977
Traitement direct des donn6es de signaux electrocardiographiques aviens et les techniques utilis6es pour I'obtention des signaux par conditions de bruit dlev6 Somraaire---Certaines m6thodes de traitement direct par ordinateur ont et6 mises au point pour obtenir le r6gime cardiaque d'un oiseau ~. partir de son 61ectrocadiogramme t616mesur6 (e.c.g.). A cet effet, on utilise un mini-calculateur P D P - 8 / L avec horloge en temps r6el et un simple affichage du type-xy. Des m6thodes d'auto-corr61ation ont 6galement 6t6 d6velopp6es pour le traitement des signaux susceptibles d'6tre obscurcis par l'61ectromyographe ou les interf&ences 61ectriques. Ces techniques ont pour effet de produire des r6sultats fiables sous certaines conditions ou les m6thodes plus conventionnelles ne le sont pas.
Rechnerabhi~ngige Verarbeitung der ekg-Signale von V6geln und Verfahren zur Ableitung dieser Signale bei starkem Rauschen Zusammenfassung--Rechnergesteuerte Computer-Verfahren wurden entwickelt, um den Herzschlag von V6geln aus FermeB-EKG-Daten abzuteiten. Ein Minicomputer vom Typ P D P - 8 / L mit Realzeituhr und einfacherxy-Anzeige wird zu diesem Zweck verwendet. Es wurden automatische Korrelationsverfahren entwickelt, die Signale verarbeiten, die durch Elektromyogramm- und Stromst6rungen iiberlagert werden. Diese Verfahren zeigen zuverl~issige Ergebnisse unter Umst~inden, wo konventionellere Verfahren unzuverlfissig sind.
Medical & Biological Engineering & Computing
November 1977
665
