1

2 2

 

'vvw‘vw u

    

      

flu. 2,: .

   

        

    

.2?
.22 .2 -
. ,2 .. .
.2. .1.’
IO u .u c2
2 ..... 2.
I. .6 .II
.I 224......
I. . .. .
(a 2... .r.
. P2112. , .
u 2..-”. .U .l 5». $9.
2.2.. 2...:
.. . .

wi-fi:

2&2...

l.
2 .o I.

. I .r
.. I. v . 2- l . . a...
2 2 2. .ar.
20 . .5...- h...‘ In? Lv\l2
h... D‘dnlhh‘.'clfz

gsumhy. b!"

Mill's}.

O

      

 

4...... .

 

L'\.|.

 

“\I

  

 

   

 

 

 

       

 

 

. .v . 2 . I . . . 2. .51.; it! I
.. . 2 2 I.u2|o.I.(.c»....2|.I..r.ru¢! ell
.. a2 2 I . I .92.?!2.n;r:).rnf.-I . ‘0...
. 2 . . ..I .2 2 .L 3.2.0.31...— -
.2 . . 3 . . . . .. 23.253.31.122» .
.. 2 2 2 2 2. . . 2; 29‘ .22.:
. y . 2 . 2 . ..
2. 2 . . . . .
. 2 2 2
L. 2 .22 . .
. 9 . . 2222 ..:
. 22 5 , I. 2 3.1:...) ..
I I 2 II ‘ a .i 2 ‘r
2 2 . . . xuwi.unn.hn.._.tfb§x'lnt
. :‘altitlhtl‘a' 2 I.
2 2 2 , ii,l.lip'2 l.‘
I 2 ‘C V .Y.
. . . 2 . . . 2.
2 2 2 L
. . r I. 2 t
. . . O .2
u 2. l2 ‘ .v- -
.L' 2 2 . v 2. o: .
2 . . 2.21.. 2.
. 2 2 . . . 2 .
. . ...- . .. 2. .
a . I . . .. . f2 .
22v 2.: 02... I2m.l.2 . .. ..
a , 2 o . ’ . . 2 . 67.. .
. . .15. i}. f ‘ II
c . ., V. I II vi I ’p92‘nlts ‘9‘. Z 2‘ .r? A ‘
. . 2 f n ‘f‘l‘linuu ‘4-.I'II.-\u'£‘vl-’ . I
A . l . . 5.} Mr... 2. 5.5.1.23...” .1... l
t . .., ..v (I... I .212. {but :3 113'!)-
2; . t . I I 32.24.120.532 fun”... infgl12‘
. . . . :2 .0. utb>i3 .Ilylf‘ii‘iu.‘ v|
o I . 2 .2 1255.13.15. .\I . ring. I...
2 2 5. I .1 . Ru 2 ‘.U§5§n..u.3b.:~.7.tr “to .8. .Ir".‘
. 2 I 2 . 2 .Il . 10.31231! it] fits Ill
1 II“. .221 u .. 2| .4 I I a. urn...rhl...r-~u.uu. [1..luaiitit .VA
. -c :1, (2.1.3.. ‘ . ‘I 2.
. .. 5. 355...... :2 -flauinflfi 2
I | o H: 21.. ... .2 01-2. hl,ol.....t| . .2 6| ..
h. 22.2. 2 2 2 ..v..l;... 3 .v. . on. 8....IV, .1‘2. Ilbx‘ w..l)lll 2. i '..l.2uv2..lfl >nh 3". 2 ..22. ....p:vI II . 2 I. 2 .... r . . . » 22 12.7 . ’2.an462..'.
. . u. .2 3%....» b... v0 u.‘ 2 . . .v . . .. {)2 22 L. .a. .. . . . 1.. h. Pb! in! 4...”. . . 2 .274. .nilnwnol. £2US‘D,,.:.‘§$. l «2...... L :0 «i! . .11.. 3.9,. 2
l ’I I 2' n 2'22"; ’2‘ ‘2}! 2 I'D i. '1

 

 

 

«MEN‘S

2 5/ 2.2:; 7 M

v

E UNIVERSITY LIBRARIES

2.2;.- 22222222222222
Michigan State
Universal—L

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

This is to certify that the

thesis entitled

SIGNAL PROCESSING FOR A mN-CONTACTING
HEARTBEAT DETECTOR AND ESTIMATOR

presented by

Shawn David Hunt

has been accepted towards fulfillment
of the requirements for

 

 

M. S . degree in Electrical Engineering

Major professor

Date—MW

I

0-7639 MS U is an Affirmative Action/Equal Opportunity Institution

._.. v____-_-___~.—--v-_._fi .v - ,7 7, ..__ _ . . fi .i, W. 7. i _ .7—

PLACE IN RETURN BOX to remove this checkout from your record.

TO AVOID FINES return on or before date due.

r_——————=—-————
PATQESDME DATE DUE DATE DUE

 

 

 

 

J

SE2 2 o 2000

 

 

 

 

 

 

.T

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

MSU Is An Affirmative Aetlon/Equel Opportunity Institution

 

 

 

SIGNAL PROCESSING FOR A NON-CONTACTING
HEARTBEAT DETECTOR AND ESTIMATOR

By

Shawn David Hunt

ATHESIS

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

MASTER OF SCIENCE

Department of Electrical Engineering and System Science

1989

.1”

(,7 a

q

7. (3/? '7 ,
VF

ABSTRACT

SIGNAL PROCESSING FOR A NON-CONTACTING
HEARTBEAT DETECTOR AND ESTIMATOR

By

Shawn David Hunt

An algorithm for a non-contacting heartbeat detector is
described. This algorithm is used in a system consisting of a
microwave transceiver, an analog section and a digital section.
The algorithm is responsible for detection and estimation of a
heartbeat using the microwave return signal. It must be able to
distinguish heartbeats mixed with clutter from clutter alone for
detection. After detection, the signal with heartbeats has neither
a known or constant shape and period. The algorithm used
combines non-linear filtering with partial correlation and
pattern recognition. The non-linear filter is used to produce

non-zero regions where heartbeats are likely. One of the non-

' zero regions is chosen as being most likely to be a heartbeat and

is correlated with the other non-zero regions. If there are non-
zero regions that are similar in shape and almost equidistant,

these regions are said to be heartbeats.

ACKNOWLEDGMENTS

I would like to thank Albert Roseiro and Dr. Marvin Siegel for their
support.

ii

TABLE OF CONTENTS

Introduction .........................................
Algorithms Previously Investigated .......................
New Algorithm ......................................
Pre-Filtering Algorithms ...............................
Derivatives of a discrete sequence ...................... 21
Introduction of the Derivative Filter ......................
How the Derivative Filter reshapes signals ..........
Improvements and Modifications of the Filter .........
Comparison of two variations of the Derivative Filter . . .
Use of the derivative filter in detection and estimation ........
Sampling frequency ..............................
Determination of normalization factor ...............

Determination of average number of zeros surrounding
each non-zero section ............................

15

24

29

31

35

45

45

48

55

56

Conclusions and Recommendations

Appendix A: Algorithms used in the Results section

iv

LIST OF TABLES

table 1 ..............................................

table 2 ..............................................

74

74

LIST OF FIGURES

figure 1 Heart data and it's Autocorrelation .................
figure 2 Heart data, power spectral density,

psd of absolute value of data ..................
figure 3 Heart data, zero crossing intervals ...............
figure 4 Block diagram of algorithm using

median-PSD-autocorrelation ....................
figure 5 PSD of median filters and unfiltered data ..........
figure 6 Corresponding autocorrelation of signals in figure 5
figure 7 Median filter operation ..........................
figure 8 Heart data, output of median filter ................
figure 9 Example outputs of non-linear filtering .............
figure 10 ............................................
figure 11 ............................................

figure 12 ............................................

figure 17 ............................................
figure 18 .............................................

figure 19 Heart dat, output of the first and the second variation of
derivative filter ............................. 34
figure 20 ............................................
figure 21 .............................................
figure 22 .............................................
figure 23 Output of derivative filter for 6.4 Hz sine wave inputs .
figure 24 .............................................
figure 25 Sample outputs of first and second variations of the
derivative filter .............................
figure 26 ............................................
figure 27 ............................................
figure 28 different sampling rates .......................
figure 29 quantization intervals ........................ 49

V1

11
13
14
15
16
20
24
25
25
26
29
31
31
32
33

35
39
40
41
42

44
45
46
47

figure 30 quantization intervals ........................ 50
figure 31 normalization values ..........................
figure 32 normalization values .........................
figure 33 Different normalization procedures ...............
figure 34 Original signal filtered and normalized to 100, 8, and 4 .
figure 35 .............................................
figure 36 N on-zero segments after filtering ..................
figure 37 Examples of derivative filter output ................
figure 38 Derivative filter output pointing out possible heartbeats
figure 39 Correlation output ............................
figure 40 Derivative filter output determining where to examine
correlation .................................. 63
figure 41 Peaks of correlation in windows determined by
derivative filter .............................
figure 42 Correlation output showing threshold and heartbeat
peaks .....................................
figure 43 Algorithm used for pattern recognition ...........
figure 44 Complete processing algorithm ..................
figure 45 Example of processing on human data .............
figure 46 Example of processing on human data .............
figure 47 The test signals with no noise added .............
figure 48 Signal 1 plus maximum noise for algorithm 1 .......

Vii

50
50
53
54
57
59

61
62

64

65
66
67
68
69
75
76

Introduction

This project is concerned with the remote detection of vital life signs in
human subjects. The purpose is to develop a non-contacting heartbeat
detector. The overall system consists of a microwave transceiver, an analog
section and a digital section. The transceiver generates, transmits and
receives a microwave signal. The analog section amplifies and filters the
microwave return signal. The digital board then samples this signal and
processes it to determine the heart rate. This thesis describes the different
methods tried for solving the detection and estimation problem, and analyses
one of the digital signal processing algorithms used for detection and
estimation.

This research was done under a grant from the United States Navy
which defined the goals of this project. They want a compact, portable
instrument for the detection and estimation of heart rate. This is needed for
military personnel to quickly identify in combat situations the heartbeat of
battle field casualties and estimate it's value.

A microwave signal at 10 GHz is sent out, and after detection, it is band-
pass filtered from four to fifteen Hz and digitized at a 512 Hz sampling rate.
From the return signal, the system must decide if a heartbeat is present, and if
it is, the heartrate. The algorithm works only with the sampled signal, thus
when a heartbeat is referred to, it means the digitized waveform of the return
signal containing a large chest movement due to the heart pulsation. To detect
the heartbeat, the microwave transceiver is aimed at the chest, for details see
[1]. Along with the heartbeat, breathing motion and other physiological
clutter is returned. The heartbeat, if present, must be distinguished from this

noise. This problem is interesting and difficult because of the number of

1

2

variables, how the variables change, and the amount they change. This
problem is not new and has been attacked before [2,3].

The returned microwave signal gives information about movement of
the chest wall. The motion of the heart is attenuated and distorted at the chest
surface. Also, the received signal will be different because of differing
physiology. Not only will heartbeats from person to person be different, but
the shape of the waveform from heartbeat to heartbeat of the same person
will be different. The heart is also naturally pseudo-periodic. It may remain
almost constant for some time or it may change period dramatically in a very
short time. The algorithm then must be able to distinguish clutter from the
heartbeat, make a decision of whether a heartbeat is present, and be able to
follow a heartbeat with changing shape, amplitude and period. The qualities
described make the detection and estimation inherently difficult, but the
heartbeat also has features that should help, and can be taken into account
when developing an algorithm. A feature that has been fundamental in the
algorithm implemented, uses the fact that the heart is impulsive in nature.
This helps in distinguishing it from periodic and pseudo-periodic signals that
are slowly varying. The problem in detection and estimation is that a clean
signal with a high signal to noise ratio cannot be guaranteed. This is the main
reason for the effort put forth to develop good filters and processing.

The microwave signal that returns after reflecting off the chest has
information of the heart rate mixed in with ambient as well as physiological
clutter. In this case everything not pertaining to the heart is referred to as
clutter. Because of the clutter and variability, having an algorithm work
under normal conditions without any pre-filtering is difficult. Thus methods
such as peak detection or correlation cannot be used reliably without pre-

filtering. The first thing to be done was to find or develop a filter that can

3

reduce the level of clutter while leaving enough heartbeat information in the
signal to make the detection/estimation problem easier to solve.
The organization of the remaining sections is:

(a) A summary of the efforts of others involved in this research before
me, and my own previous efforts.

(b) A summary of the algorithms tested for non-linear filtering.

(c) A description of the derivative filter.

((1) Justification for different parameters.

(e) Detection and estimation after the non-linear pre-filtering.
Different parameter settings are justified and typical results are
given.

(f) Results to compare the final algorithm with four other algorithms.

(g) Conclusions.

The appendices contain the implemented algorithms, written in C, and

some sample four second data segments used for the result section.

Algorithms Previously Investigated

The first algorithm implemented was autoconelation of a four second
segment [4,5]. Because the heartbeats are not purely periodic, and vary in
shape, this method does not perform as well as is needed. An example of a
segment and it's autoconelation are shown below. As is apparent, if the
signal is easy to pick out visually, and is almost periodic, autoconelation will
perform well. However with a poor signal, as that shown in figure 1,the

autoconelation method does not work.

 

300

200 J

.001 0
o I.

1

'300 -‘

. Afifi

)e+6

 

 

 

 

 

 

 

J
)e+6 -l

)e+6 -

)e+o - WWW“.-
)e+6 a

q

3406 ' I r r ' I '
0 l 2 3 4

time (sec)

 

 

 

 

figure 1: Heart data and it's autoconelation.

5

W. Byrne used an adaptive filter approach [6,7]. Pat Mahoney is
presently continuing this work on an adaptive filter algorithm.

G. Hoshal implemented a number of algorithms [8]. One algorithm used
a spectral approach in the frequency domain, one used a comb filter in the
frequency domain [9,10], and another used features to form a multi-
dimensional classification space [11].

The spectral approach uses the Fast Fourier Transform of the input
signal and estimates the heartbeat frequency from that. It was found that the
breathing interfered if the signal was not filtered by a low-pass filter. It was
also found that the return signal due to the heart has many strong harmonics,
but that due to the breathing does not. Because of this, the signal is band-
passed filtered from 4 to 15 Hz. The fundamental component of the heart
signal is blocked out along with the fundamental of the breathing. The heart
signal is detected by it's harmonics. The breathing plays little part because of
it's very small harmonics. This procedure is very effective with a strong
heartbeat present, but is not good otherwise. He took the absolute value of the
signal before performing the FFT.

e+l

 

e+t
e+0d
e+H

q
e+rJ

 

 

time (sec)

 

e-I-

 

 

b)-.. u-J\

2-1

 

 

 

 

 

C )-’>n I

I
0 l 2 3 4

 

 

' V

I U r. tT‘f r. I 1 I I I ‘fi' ' '

l
6 7 8 9 10 ll 12 l3 14 IS [6
frequency (Hz)

I
5

figure 2: a) Heart data, b) power spectral density, c) psd of absolute value of
data.

7

Another approach Hoshal used was an adaptive comb filter. It used
adaptive signal filtering to select the best fit and from this determine the
heartbeat.

The final approach taken by Hoshal was to select several features as
different dimensions of a classification space. Data from many subjects was
gathered to determine the thresholds for decision regions in this space.Thus,
the input signal could be processed to determine the features and from these a
decision about the presence of a heartbeat is made. Many features were
investigated, such as the average value of the power spectral density, and zero

crossing count.

New Algorithms

Since the previous signal processing techniques were examined, a
number of additional processing methods have been investigated. A statistical
approach using correlation was tried first. This uses two features from the
autocorrelation to form a two dimensional decision space. The two features
selected were the sequential distances between the three largest peaks. The
algorithm however will work with any feature from the autoconelation.
After the two features of the signal autoconelation are selected, these
features are computed from as many heartbeat files as possible. From the
pool of these two sample features, a correlation coefficient and a regression
line was determined. When new data is acquired, it's autoconelation
features are compared with the regression line and a decision is made about
detection. Because the features selected give a heartbeat estimate, both
heartbeat detection and heartbeat estimation were be determined at the same
time. This works reasonably well with a good signal and windowing. The
windows work as follows. If it is decided to detect a heartbeat between 60 and
120 BPM (beats per minute) we can segment, or window, a section of the
autoconelation corresponding to this frequency. The corresponding
frequency for 60 to 120 BPM is 1 to 2 Hz, or a 1 Hz difference. Windows
would be placed at intervals corresponding to every 1 Hz. Separation
between the largest peaks in every window should correspond to the heart
frequency. The drawback in estimation range is immediately apparent. For
non-overlapping windows, the highest frequency estimated is twice that of
the’lowest frequency. For example, if 50 BPM is the lowest frequency, then
100 BPM is the highest. However the largest drawback for this algorithm is

8

9

the correlation technique itself. To extract the wanted features, an
autocorrelation must be done before any other processing. Practically, for
autoconelation to work, the signal must be periodic and the noise level must
be small. The heartbeat is not always periodic, and large changes in heartbeat
frequency in a short time can occur. Typically, four seconds of data sampled
at 64 Hz (256 samples) have been taken. If in these four seconds the heartbeat
remains approximately constant, and the signal is relatively noise free, the
correlation will work. These are ideal conditions however, and will not work
satisfactorily in most situations. Typical heartbeat data is shown along with
it's autoconelation in figure 1.

The difficulties with autoconelation required new algorithms to be
developed.

There is a similarity between estimating pitch in speech and estimating
the heart rate. The problem of estimating the fundamental frequency of
speech has been widely studied and the ideas that worked well in that area can
be easily tested to see if they have an application for the heartbeat estimation
problem. Four main algorithms were tested.

Methods using zero crossing count and zero crossing intervals seem
good choices [12,13] because the zero crossing count is lower in the segments
where there is a heartbeat, and the zero crossing interval is almost constant
from heartbeat to heartbeat. An algorithm to determine heartbeat frequency
by zero crossing interval was developed, and the approach had no problem
with aperiodicity. However, it was found that the differences of the intervals
of the signal containing a heartbeat and signal containing no heartbeat did not

lead to a good detection or estimation scheme, as shown in figure 3.

10

 

2e+3

2e+3 '-

Ze+3 -

2e+3 -'

2e+3 -

)Afi?

10

 

 

 

 

 

 

 

"NO v I f I '
0 l 2 3 4

time (sec)

-
1

figure 3: Heart data, zero crossing intervals.

The zero crossing count did appear to lead to a reliable detection
scheme. It was not used because it did not lead to a detection scheme and the
algorithm used now does both detection and estimation.

The cepstrum approach was also tried [14]. The Fourier transform of
the signal is taken. The square root of the power spectral density is taken and
signals that are convolved can be separated if they are separated in the
frequency domain. This approach is good for data where the noise is
convolved with the signal. This was not the case and this approach showed no
improvement over simply looking at the power spectrum.

The zero infinite clipping algorithm with autoconelation was also
investigated [14,15]. First, the mean is taken out of the signal. A threshold is

11

then determined, usually some percent of the average signal strength. All the
elements above this threshold are assigned a one and all the elements below
the threshold are assigned a zero. The signal is then autoconelated. Although
this method uses the zero crossing interval to some extent, it performed
poorly because of the use of autoconelation in a signal that is not periodic.

The last technique discussed in this section is a combination of non-
linear filtering and frequency domain pattern recognition. First the 256
samples of the signal are filtered with a non-linear filter. 768 zeros are
inserted after the data to make a segment of length 1024. An FFT is then
taken. After the power spectral density is determined,an autoconelation is
done to detect heartbeat harmonics in the frequency domain. The process is
shown in figure 4.

 

 

 

 

of PSD

 

 

 

Median 7 ' 1
Filter Signal Insert Zeros H Compute PSD H Autocorrelation

 

 

 

 

Use threshold

 

If heartbeat is present. on autoconelation
determine frequency ‘ to determine if
from autoconelation heartbeat is present

 

 

 

 

figure 4 : Block diagram of algorithm using median-PSD-autocorrelation

The non-linear filter used is the median filter. It's operation and sample
input and output are shown in the next section.
After the filtering, the 256 data elements are zero stuffed to complete

1024 elements. The longer segment is used to give a precision within i0.125

12
Hz. An FFT is performed and the power spectral density is determined.

Because the heart is pseudo periodic and impulsive in nature, it has
harmonics that are periodic in the frequency domain. The signal is band-pass
filtered from four to fifteen hertz before it is sampled, so the fundamental
frequency of the heart is almost entirely attenuated before the processing
starts. The heartbeat harmonics remaining are used to determine the
heartbeat frequency. This pre-filtering will eliminate interference with
frequencies near the heartbeat fundamental, such as breathing. If they have
very small or no harmonics above four Hz, they will not interfere with this
detection scheme. The heartbeat harmonics are evenly spaced and the spacing
is the fundamental frequency. An autoconelation is done on the power
spectral density. Looking for the highest peak in the correlation file
corresponding to a fundamental frequency of .5 Hz to 2.5 Hz, a threshold is
used to determine if a heartbeat is present. This method does not appear to
suffer as much from an aperiodic signal as a straight correlation of the time
signal. This change in period between heartbeats in the sampling time will
result in wider harmonic peaks in the power spectral density. This will in
turn lead to a less accurate estimate, but the estimate is still possible,
something not seen in correlation.

An example of the PSD of the original signal and the filtered signal is
shown in figure 5. The difference of the two large peaks in a) of figure 5
correspond to the fundamental heart frequency. Few peaks in the psd mean
fewer peaks in the autoconelation function file and choosing the conect peak
is easier. If the unfiltered signal's psd and thus it's autoconelation has more
peaks, it makes the peak corresponding to the fundamental heart frequency

much more difficult to detect as shown in figure 6.

l3

 

3e-2
3e-2 1
2.2—2 1
.3e-2 4
3e-2 .

De-Z '

 

 

 

 

2e- I
)e-I -
)e-2 -

)e-2 -
)e-2 :
)e-2 -
not) a

O

 

 

 

 

 

 

 

hr.

1 6 7 10 11 12 13 1415

la
2 3 4 5
frequency (Hz)

figure 5: PSD of median filtered (top) and unfiltered data.

 

16

14

 

3e+0

3e+0 -

Oe-i '1

Oe-l a

Oe-i '1

Oe-l '

1

 

n-0n

 

 

 

 

 

leto

 

 

 

 

‘iitfb

. T . . f
I 2 3
time (sec)

i

figure 6: Corresponding autoconelation of signals in figure 5.

 

 

Pre-Filtering Algorithms

To make the heartbeats easier to detect in the clutter, pre-filtering
algorithms were developed. Like the algorithms previously implemented, the
filtering is done on a four second segment of data. In the analog section, a
band-pass filter from 4 to 15 Hz is used, but the signal is still very cluttered
and any additional linear filtering by the digital section is of little use. Of
course the device cannot extract more information than is present in the
signal,but we want the pre-filtering to reduce the clutter while leaving
enough of the heartbeat signal so the detection and estimation can be done
easier than before. A first try in using a non-linear filter was the median
filter [l6,l7,18,19].

This filter is a special case of the L-filters. It works as shown below. In

figure 7, each element is filtered using the surrounding elements.

Original Signal Order Output
...42975... 429—249 ...477...

LE3 29—279

975-579

figure 7: Median filter operation.

The elements used in the determination fall into a window. This window
is simply the number of elements used in the ordering. The median filter will
preserve discontinuities in the signal if they are long enough, yet eliminate or

greatly attenuate short discontinuities. Thus the window is chosen according

15

16
to the width of the peaks to be eliminated.The elements in the window are

ordered in ascending order of magnitude. The element now in the center
replaces the original element. This is done for each element in the data file.
For elements on the end, the missing elements in the window can be chosen as
zeros, or the value of the end element.

The effect of the filter with a window of 3 is shown in figure 8.

 

4o

 

 

 

 

 

 

I V I U l' v

2 3
time (sec)

0-:
5.1

figure 8: Heart data, output of median filter.

This filter has the desirable feature that it keeps sharp edges of peaks
wider than some threshold, while reducing peaks narrower than the
threshold. It passes any peak wider than the threshold and it severely distorts

the waveshape as shown. As previously discussed, this filter was used along

17

with frequency domain pattem recognition. This pre-filtering did improve
detection and estimation, but the complete algorithm did not perform
adequately. After this attempt, detection and estimation based on shape was
investigated. Because the median filter did not generally make the signal
easier to detect visually, it was not used. It did show that non-linear filters
have promise, and gave reason to pursue them further.

The intent was to filter the signal so it would be easier to detect. It did
not matter if the signal shape was totally changed, as long as it emphasized the
heartbeat signal. Many different methods were tried. The first and simplest
was simply a threshold (peak detection). The next used only the direction of
change not the amplitude. An adaptive filter using the error as the output was
also implemented.

The first method set all the samples smaller than a threshold to zero.The
threshold was determined as a percent of the largest sample in the four
second segment. This method did eliminate a lot of clutter, but it also
eliminated a lot of signal information and small heartbeats. It did not
emphasize the heartbeat over the clutter already present in the signal. An
example of input and output is shown in part b of figure 9.

The next filter only used information of the direction of change from
sample to sample. This filter starts at zero and goes up one unit if the
direction of change from the first sample to the second is up, down one if the
change is down and stays at zero if there is no change. Because there is only a
change of one or zero each time, it removes information of amplitude
change. It continues going up one, down one, or staying equal, according to
the change from sample to sample. An example of this filtering is shown in

part c offigure 9.

18
The problem is that the heartbeat can fall in the same range of

frequencies, can have the same magnitude, and have no simple shape
difference from the clutter. The solution found here was to go back to the
basic assumptions about the hearLA feature of the heart is it's impulsive
nature. It also must have enough amplitude to make it distinguishable in the
clutter. We require it to be at least as large as the clutter. Since the derivative
will give information about the impulsive nature, using information from
derivatives seems like a good idea. Both the adaptive filter and derivative
filter do this.

The adaptive filter is based on a difference equation where the
coefficients are determined using the least square error method. Different
order equations were tested, the fourth order was found to perform the best.
Higher order equations did not improve the performance. An example output
is shown in figure 9d.

The final test is how well it performs. The adaptive filter worked well,
but because the derivative filter appeared to perform better, in this report,
the derivative filter was used.

In the next section this derivative filter is discussed.

19

 

 

 

 

 

 

 

 

 

1e+1
2e+i
)e+0-
3e+1 -
1e+1 -
‘WU I I I f
0 2 3
time (sec)
azheart data
1e+1 *
3e+l
ie+o-
e+i -
e+i -
WU I V I I I
o 2 3
time (sec)
bzthresholded
3e+l
eui
le+0
WU I 1 I ' I
0 2 3
time (sec)

czdiscarding amplitudes

 

 

 

20

 

.e+|
2...]
fe+1-i

fl
re+04

 

 

 

 

91-2
1

9+] ..

 

 

 

' I ' l ' I T I
O 1 2 3 4
time (sec)

dzadaptive filter

“Adm

d

figure 9: Example outputs of non-linear filters.

Derivatives of a discrete sequence

The derivatives, or the discrete equivalent of a derivative will be used
extensively. To show how the absolute value in the derivative filter relates to
the discrete derivative, the first and second derivatives of a sequence are

shown.
Start with a Power series expansion of a function f(z) about a pornt a

f(z): f(a)+ Lit (a)+ +S—l-r"(a+) £2_-!.L3f "'(a) +

(1)

(1) can be put into two other forms.

first:

z=x+h h=z-a
x=z-h a=z-h
x = a

substituting into (1)

h 0 hz " "I
f(x+h)=f(x)+ Ff (x)+ 57f (x)+§;f (x)+ . ..
(2)
or:

again substituting into (1)

h , h2 ,, h3 m
f(x-h) = f(x) - If (x) + 271' (x) - yf (x) + . . . (3)

combining these we get:

21

22
(2) - (3)

f(x+h) - f(x-h) = 2f '(x) h + 2f "'(x) %?+ . . .

- - 2
f '(x) = f(x+h)2hf(x h)_ [ f'"(x) 1:231?+ . . .]

f(x+h) - f(x-h)

= 2b + CW)

(2) + (3)
h2 . h4
f(x+h) + f(x-h) = 2f (x) + 2f "(x) -2-!~+ 2f W(x) 47+. . .

f(x+h) + f(x-h) - 2f(x)

. h2
112 -[f1V(x)4-!+...]

f "(x) =

f(x+h) + flg-h) - 2f(x) + o (112)

The error terms in our application however are actually not important
since we do not want to estimate the derivative of the function, but merely use
the estimate to make the heart signal easier to detect. Because the output
signal will be used for detection and estimation, it is important to look at the
performance of the proposed filter.

In the following sections, the derivative filter, and two variations will
be introduced. How the filter and the variations reshape the signal is
investigated. The conditions for the heartbeat peak to be higher than the

surrounding noise is derived, and the conditions for the highest peak in the

23
pre-filtered signal to be the highest point after filtering (having the peaks be

filter invariant) are determined.

Introduction of the Derivative Filter

Inputsequence
. . .x0xl x2x3x4x5x6...

Output sequence
...y0yl y2y3y4y5y6...

The derivative filter is defined as follows:
yi = Xi * [(Xi+r - Xi) + (Xi - Xi-r)]
The output sequence is uniquely determined by the input sequence.
The first filter used can also be written as:
yi = Xi * [(Xi+1 - Xi) + (Xi - Xi-1)l

= xi * [xi+1 - xi-1] - the present value is weighted by the
difference of the adjoining values

The filter in this form does not work well for 3 main reasons.

1)

xi

xi
. A Xl+l , X1+l
x1—1 Xi-i

figure 1 0

24

25
The definition used for the discrete derivative does not take into account the

middle value; the sequences shown in figure 10 have the same difference

multiplier.
2)
X 1
X1- iAxi+ 1
figure 1 l

The difference in figure 11 is zero, leading to a zero output, even though the

value and the change between it and the surrounding values are not zero.

3)
x2
x3
x1 0
\X4
x0

figure 1 2

even if x1 and x3 are both positive,as shown in figure 12, the corresponding

outputs yl and y3 will be of opposite signs.
The absolute value is introduced to help resolve these problems.

yi = Xi * [IXi+1 - Xil + IXi - Xi-rll

26
Here the middle value is multiplied by the absolute value of the

surrounding differences or the middle value is multiplied by a factor
determined by the neighboring points.
There exist five possible three point sequences. Taking only direction

of change into account.

A / e e e
<1) <3)

(2)

(4) (5)

figure 13
1) if x3<x2 & x2°>xl
lx3-x2l + lx2-xll = x2 - x3 + x2 - x1= 2x2 - x3 - x1

2) if x3>x2 & x2>xl
lx3-x2|+lx2-x1|= x3 - x2 + x2 - x1= x3 - x1

3) if x3 = x2 & x2 = x1
lx3-x21 + lx2-xll = x2-x3+x1-x2 = 0

4) if x3>x2 & x2<x1
lx3-x2|+lx2-xll=x3 -x2+x1-x2=x3+x1-2x2

5) if x3 < x2 & x2 < x1
lx3-x2l + lx2-xll = x2-x3+xl-x2 = xl-x3

Cases

27

A relative maximum, is weighted by the absolute value of the
second derivative.

Strictly increasing.
This is the same as the original derivative formula.
yr = X*X' = XiIXi+1 - Xi-rl

All values are equal, so the derivative, and output of filter is
zero.

A relative minimum point.

yi = X*X" = Xi[(Xi+1 - Xi)+(Xi-l - Xi)l = Xi[(Xi+1 + Xi—l - 2Xi)l
Using half step points (denoted XiiI/Z). the second derivative can
be calculated as

Xi' = [(Xi+r/2 - Xi) + (Xi - Xi-1/2)l = Xi+1/2 - Xi-rrz
Xi+1/2' = [(Xi+1 - mm) + (Xi+1/2 - Xi)] = Xi+1 - Xi
Xi-rrz' = [(Xi - Xi-1/2) + (Xi-r/z - Xi-1)] = Xi - Xi—l

Xi" = [(X'i+r/2 - X'i) + (X'i - X'i-1/2)l
= [(Xi+1 - Xi) - (Xi+r/2 - Xi-1/2)l +
[(Xi+1/2 - Xi-r/z) - (Xi - Xi-1)l
= Xi+r + Xi-r - 2Xi

Neglecting the h2 factor, this is the same as the second
derivative shown in the previous section. X2 is weighted by the
second derivative.

When strictly decreasing, we have a negative derivative, this is
the negative of the original formula.
yi = -X*X' = -Xi[Xi+1 - Xi-rl = XilXi-r - Xi+1l

28

From these five cases it is seen that when the samples are strictly
increasing or strictly decreasing, the output is X, multiplied by the absolute
value of the normal discrete derivative. When it is a relative minimum or
maximum, where the first derivative becomes very small, or goes to zero,
the value is multiplied by a factor of the second derivative. The sample will
not be multiplied by zero even when the first or the second derivative goes to
zero. The only case of a zero will be when three consecutive points are equal,
when both derivatives are zero.

29

How the Modified Derivative Filter Reshapes Signals

Two of the five possrhle three point sequences are the inverse of the

other, so only three need to be investigated. Because three consecutive equal

points result in a zero output, this case will not be investigated. To be

determined is what kind of signal, in terms of relative amplitudes and relative

derivatives, is needed for the heartbeat detection and estimation scheme to

work. This criterion will also lead to choosing the best variation of the filter.

Because the intervals between samples is constant, the relative amplitude is

directly related to the derivative.

To get a intuitive feel for the filter, here are some sample of how it

reshapes three common signals.

input

l [.1

r‘Iril m
U U

 

Wit MUU

output
comments:
the filter will severely distorts
not pass regions where

constant levels derivative is small

figure 14

/\/\
V
/\/\
V

constant amplitude
change

30

Smooth shapes and constants will not be passed, those with constant change
will not be reshaped.

31
Improvements and Modifications of the Filter

We want to emphasize the heartbeat signal. In our signals, the largest
derivatives are where the signal is the smallest. Where the signal is large the
derivative is small. This leads to output signals that look like:

before filter after filter

figure 15

This leads to a trade off. To get a good approximation of the signal, a
high sampling rate must be used. Thus the points with the largest amplitude
will have small derivatives. If a low sampling rate is used, the derivative will

not be small, but the peak might be missed as shown in figure 16.

figure 1 6

32
It is possible to improve the situation in two ways. One uses a clever data
reduction scheme, the other uses a new method to calculate the derivative.

The data reduction scheme could be as illustrated below:

figure 1 7

This method attempts to approximate the signal with a piecewise linear
funtion. It would help in the problem of peak reduction, but it is unattractive
because the spacing between the samples may not be equal.

A new method to calculate the derivative is as follows: Use only the
extreme points to calculate the derivative, then multiply all the samples

between these extreme points by the same derivative.

_ yi=xi*a i<i<12

yi=xi*b11<i<18

 

figure 1 8

This method retains much of the waveshape, keeps the spacing equal and
eliminates the dips encountered before. Keeping the spacing of the samples
uniform is necessary for later processing, so only the derivative filter using
the absolute value of the differences of the 2 surrounding points, and the
derivative filter using extreme points to calculate the derivative will be

compared. These will be designated the first and second variation,

respectively. An example of both filtering methods is shown in figure 19.

34

 

 

 

6*2-

 

e+21

 

a)nat’3
6+1

ee]:
2 11

e+0-~,

 

 

 

 

 

e+i-

 

b)...

e+1

 

 

 

 

 

 

 

WU| r r f I
0 l 2

time (sec)
figure 19: a) Heart data, b) first variation, c) second variation filter outputs.

35
Comparison of two variations of the Derivative Filter

This filter is used to point to the heartbeats in the original unfiltered
segment, so we would like the peaks to be filter invariant.

The first variation: using the absolute value of the surrounding points

Let xi.2 and xi be fixed and let xi.1 move on a vertical line located at -h.
This constraint applies in this system because the interval between samples is
uniform.
Let xi be bounded below and above by x24 and xi . This section will

determine the amplitude constraint on x“ so that xi will be remain higher

after the filter.
L—a
b
’2h "h

0

 

figure 20

If xi.2 (at -2h), xi (at 0), and xi.1 (at -h),all fall on a line, a triangle is
formed, and it has been shown this will pass with only a scaling factor. If xi-1

does not fall on the line formed by xi-2 and X, then

Ixil=a a-b=v
v
lxi-1l=c y=fit+a

lxml = b

For the first variation:

yi = xi * [lxm - xil + lxi - x141] - to find the positive peaks, find the
zeros of the first derivative and
negatives of the second derivative.
. (ym - yi) + (Yr ' Yi-r)
YI = 211

because this derivative never, or almost never, goes to zero, will
define critical points as points where (ym - yr) and (yi - y“) have opposite
signs.

to _ (Yttfiji-l -ZYI)
Y! -' 112

when yi" is negative there is a local maximum and when yi" is positive
there is a local minimum. We would like these local maximums and .
minimums to be the same samples as before, only the relative heights of the
samples changed.

' ' + " ' .. M
._ 9111-1 - Yr) + (Yi- Yi-r) .. _ W)
y‘ — 2h Yi " hz

 

considering only positive peak points,

Yi' =\f([xi+l * (IXi+2 - Xi+1| + lxi+1- Xil)l - [Xi-r * (lxi - Xi-rl + IXi-i - Xi-2')],2h)

37

case 1
Xi > 11-1 > 11-2
Xi > Xr+l > Xr+2

this extreme point in the x] series will remain the extreme point in the
y] series, if a and b are of different signs.

yi+1 - Yr = [Xi+r(lxi+2 - Xi+r| + IXi+1 - Xi|)] - [Xi(lxi+1 - Xil + lxi - Iii-11)]
= [Xi+r((xi+r - Xm) + (Xi - Xi+1))] - [Xi((Xi - Xi+r) + (Xi - Xi-r))]
= xi+i(xi - xi+2) - Xi(2Xi - Xi+r - Xi-l)
= 2Xi+rxi - inxi - xi+lxi+2 + XiXi-l

Yr - Yr-r = [Xi('Xi+1 - xii + lxi - Iii-11)] - [Xi-1(lxi - xi-ll + lXi-r - xi-2')l
= [Xi((Xi - Xi+1)+ (Xi - Xi-r))] - [Xi-1((Xi - xii) + (Xi-1 - Xi-2))]
= Xi(2xi - Xi+r - Xi-r) - Xi-1(Xi - M2)
= 2Xixi - 2Xixi.r - Xi+1xi+r + Xi-1Xi-2

to get:
(Yi+1 - Yr) < 0
(Yr - Yi-1)> 0

will need:
2Xi+rxi - inxi ' Xi+rxi+2 + XiXi-r < 0
2xixi - inxi-1 - xixm + xi.1xi.2 > 0

2Xi+rXi - 2XiXi - xi+lxi+2 + XiXi.r < 2XiXi - 2XiXi-r - Xi+rXi+1 + Xi-1Xi-2

3Xi+rXi + 3XiXi-r - Xi+rXi+2 - Xi-lXi-Z < 4XiXi (5)

assumed: xi > xi“ > Xi+2
Xi > Xi-r > Xi-2
introducing two more constraints,
xi+1 = xi-l
xi+2 = xi-2

38

(5) becomes,
6Xi+1xi - in+rxi+2 < 4Xi2
3Xi+1xi - Xi+1xi+2 < 2Xi2

 

3C3 - Cb - 232 < 0 01‘ 3xi-1xi - xi-1xi-2 < 2Xi2
c(3a - b) - 2a2 < 0 xi.1(3xi - x”) < 2x22

262 .221;
C< 3a - b Xi-l< 3xi-xi-2

-2
So x“ must be less than 3%; for the peak to remain in the same
I ' 1-

place. Of course, the same can be done for xi” the first sample after the peak.
On the negative peaks the inequality is reversed.

It becomes apparent from this equation how nonlinear the process is.
Suppose xi: 10 and xi-2 = 5. xi.1 must be less than 8. If xi = 5 and xi.2 = 0, the
same relative difference, xi.1 must be less than 3.33. Thus, the same
difference is changed in amplitude, the relative difference of xi and xi-1 does

not stay the same.

39
50 . 500 10. 110.
1.1 / 1.5 ’ 1.045
40 L 400 o

L 100

 

mu1t1p1es produce the same relative differences
v same xi-i values do not
100
1.1
80
figure 21

In the leftmost figure above, if xi-2 = 40 and X, = 50, xi-1 must be
less than 1.1 of xi. Any multiple of these two numbers will have the
same requirement for x“. The same relative differences, such as 0 to
1 0 and 100 to 1 10, do not have the same requirements for x“. In
general, xi-1 must be about on the line formed by xi and xi-2 or below to get
good results. This is very important. With a sampling rate to assure the peak
is recorded, this will not be the case in general.Looking at 41 segments of
typical heart data both filters were used and the results are as follows. Using
only three points to determine the derivative, 78 heartbeat peaks did not
change position, and 41 did. Using the extreme points to calculate the
derivative, 93 heartbeat peaks did not change position and 26 did.

The second variation of the filter: using extreme points to calculate
derivative.

xf

value
of

“N
slope 1325K 33's 23'

f-b
figure 22

 

AT =—?L:—§b' the slope from x1, to Xf
All horizontal distances between the points are the same. Between two
extreme points, all the slopes are either all positive or all negative, and AT is
just the average of all the slopes between these extreme points.
xi - xi-1 _ Axi
h — h

AXf+ AXf-1+. . . + Axb+1 _ Xf- xb
f-b ' f-b

The first variation uses lei+1l + leil for a multiplier, which for high

sampling rates and sinusoidal signals, is smaller than AT giving the second
filter better performance. Also,because Xf is higher by definition than the

41
other samples, and all the points between xb and Xf are multiplied by the same

derivative value, Xf will be the highest point after the multiplication.
The Axf and Axb+1 for the sinusoidal signal are smaller than the other

differences and using lei+1l+ leil for the multiplier can make some inner

points larger than the extreme points. This is shown in figure 23.

 

10

 

 

 

first variation

 

 

 

 

second variation

figure 23zoutput of derivative filter for 6.4 Hz sine wave input

Our main purpose is to increase high derivative high amplitude signals
over low derivative high amplitude or high derivative low amplitude signals.
Using the first variation we would need to sample at least as low as shown

below to get as good results as the second variation.

42

figure 24

Because we have a band-pass pro-filter, we can set it at the highest
heart frequency component so no clutter is present of higher frequency, or
has a higher derivative than the heart. There is a trade off in the filter: use a
high sampling rate to assure getting the highest point and then have smaller
derivatives; or use a smaller sampling rate to keep the large changes between
samples and maybe miss the highest point. For square waves or triangle
waves, both variations have the same output, so the higher sampling rate that
can be used with the second variation only assures the peak is obtained. In the
heart signal measurement,we deal with more sinusoidal shapes. In these
waveforms the highest slope is where the signal is smallest; we would like the
slope large for large amplitudes. To make the signal large, the filter is non-
linear. In any sinusoidal signal, with a high enough sampling rate,

XfIAT' > Xf [lAXf+1l + lAXfI]
as was shown in figure 23. This effect, using the first filter, becomes

intolerable at high sampling rates. As seen before, there can be double peaks
produced. The second filter will eliminate double peaks because, by

43
definition, the extreme point is higher than the others. In the detection and
estimation algorithm the pre-filter is used to determine places where a
heartbeat is likely to be. So double peaks are tolerable if this algorithm is
better at distinguishing heartbeats from clutter. As shown, however, this is

not the case and the second variation of the filter performs best.

44

a)

 

2e+i
1e+i -
)e+0 ~

e+i- i

’1“.

 

 

 

 

 

 

 

 

 

b)

 

2e+ I

1e+1-1

)e+0 '-

 

e+lq

 

 

 

lmb ' I ' I ' I
0 l 2 3
time (sec)

figure 25: Sample outputs of a) first and b) second variations of the
derivative filter

45

Sampling frequency

We need to determine an adequate sampling frequency for the
derivative filter. We have an 8 bit, or 256 level analog to digital converter.
We use at least fifty percent, 128 levels, of these levels for the heartbeat
peaks. Our analog filter cuts off at 15 Hz.

Using a 15Hz sine wave for our analysis, we will determine the
sampling rate needed to discretize a sample within a pre-specified percentage
of the analog peak. From our algorithm, we want to locate the heartbeat
peaks. This means no matter what the phase of the signal is with respect to the
sampling (0-21t) at t=0, we want to be assured of at least having one sample
within 90% of the highest point in the signal.

 

 

192 ?\
128 ) t
128
64 V
T
2
figure 26

for 15Hz signal and a peak to peak sinusoidal signal of 128,

46
to changes 1800/31-6sec => 5400:-

0.9(192) = 64 cos(¢) + 128 => 45.570

{KY/Y 90%line
I34 180 225

figure 27

for Nyquist criterion:
SR(samp1ing rate) > 130 sec => 30 Hz SR

for 90% criterion:

45.57o*2 _

——0—- sec => 59.24 Hz SR
5 IOUsec

If we use all 256 levels, the sampling rate for 90% criterion is 73.23
Hz. The 64 Hz sampling rate used is adequate, as shown.

47

 

2+2
3+2-

3+2 .1 J P I L

 

2+o-, ' 2
H2

H2 1

\A.” ,

a) 256 Hz sampling rate

 

 

2+2
3+2-
2+2 " 1 1

I
3’0 " 1

 

 

 

 

 

 

 

 

b)64stamplingratc

 

+2

 

 

 

 

 

+2-

 

 

 

 

 

«o0 . . . . .
0 1 2 3

time (sec)
c) 32 Hz sampling rate

figure 28: different sampling rates

 

48

Determination of normalizing factor

To make the filtered data easier to process integer normalization is
done. This division factor, 5, is chosen so that some inter-heartbeat clutter is
eliminated. Having some of the data points equal to zero, will make the
detection and estimation easier, only having to look at a finite group of non-
zero areas in each segment of data and determine if these are heartbeats or
noise. Thus we want as many points as possible to be zero after filtering,
without eliminating any heartbeat peaks. The trade off will be a greater
number of zeros vs. reduction in precision. A method for determining this
normalization and the normalization factor must be determined.

First, a method for calculating the normalized output must be
obtained.
In the presentation that follows, 3 is the normalizing factor. The x's

and arithmetic is all integer, which truncates all fractions.
1) The fast method used'is as follows.

To keep as much precision as possible, x, is rounded off to the
nearest 8.

Xi+ sgn(xi)*%

 

S

where sgn(x) = 1 if xi 2 0
-1 if Xi < 0

49
With this method, all xi S 13/21 => go to zero. In other words, a factor s
is chosen so that all samples less than half s will be zero and precision will be
reduced by this factor also.

All xi falling within 1 s/2 of some (ks) will be rounded off to the

 

sarnevalue.
xi=skfor sk-s/2 < xi 5 sk+s/2 k=0,1,2,3...
0
-2s -11/2s -s -s/2 s/2 s 11/25 25
l l I I I l | l
-2 -1 O l 2

figure 29 : quantization intervals

2) The second method is simply integer division by s.

...’.‘_i
3

There will be reduction in resolution or precision, but if the s here is the
same as the s in the first quantization method, we will have the same
quantization intervals except around zero, where it will be twice as long. We

need to decide how much loss of resolution can be tolerated.
All xi falling within is/2 of some (ks) will have the same value.

xi=skfor sk<xiSSk+s k=0,1,2,3...

50

Negative values will be the mirror image.

 

0
-2s -11/2s -s -s/2 l s/2 s 11/23 25
1 1 1 1 1 1 1 1
W
0

figure 30 : quantization intervals

The trade off is precision vs. eliminating clutter. Examining the

precision of each one:

figure 31 : normalization values
the first method has a maximum quantizing error of .53

S 28

21111214111 11:
W

-l l
O

 

 

 

 

 

 

 

figure 32 : normalization values

51
the second method has a maximum quantizing error of s

The first method has a quantizing error of $.55 and the second an error

of s. The second method however, makes everything from -s to 3 go to zero,
and the first method makes everything from -s/2 to s/2 go to zero. If the
second method is used with s/Z as the normalizing factor, it will have a

quantizing error of s/2 and will make everything from -s/2 to s/2 go to zero.

To get the mean square quantizing error of each [20] (MSQE)

L
MSQE= 2(MSQEi)Pi L is #of levels
i=1
Pi is probability of sample in ilh level

MSQEi= ] (x - 392 p(xli) dx
11

p(xli) = E‘g-i’i Pd(x) - probability density of voltage Bi

quantizing value

MSQE= ] (x - 1302 Pd(x) dx

11

.m—

i:

Integral evaluation for the first case:
Using 16 as the normilizing factor, the interval around zero will be
from -8 to 8. All the intervals have the same MSQEi.

8
3
MSQEi= £(x - 0)2 dx = x? =lg—Zi= 341.33

52
Integral evaluation for the second case:
Using 16 as normalizing factor, the interval around zero is from --16 to
16, length 32. The other intervals will be of lentgth 16. The MSQEi
shown here is the average of all the MSQEi's.

32 16
1023 2
MSQEl=l-0_251£(x - 16)2dx + T025_1-L(x)2 dx

_1023[g 2 2 16 ]3_2 2 [x_3_ -16
‘1025 3’ x + x16+1025 3 16
=1367.997

Using 8 as the nonnilizing factor, to get same zero region as the first
case,

16 8
2046 2
woes—2050 -g (x - 8)de + 50—50"; (x)2 dx

= 170.667
A Thus the second method has a smaller MSQE for the same zero area. In
s
x+sgn<x>*(§)
s

 

figure 33, the signal was normalized to :8, using x = on top,
and x =% on the bottom. The second normalization method, x =§ has been

chosen for normalization.

53

e+1

 

6+0-

 

 

«‘1

e+1

 

 

6+0-

 

 

 

 

"+01? 1I I I I I I
0 l 2 3 4
time (sec)

 

figure 33: Different normalization procedures.

Determining the normalization factor, 3 in the equations above, is a
matter of trial and error. The method used was to use many heartbeat
segments and determine the largest factor that does not eliminate any
heartbeat. Figure 34 shows the output of the filter for different
normalization factors. This is a very good signal, so a small normalization
factor, such as 4 works very well. In general there is a large difference
between heartbeat amplitudes and a larger factor is needed so none of the

heartbeats is normalized to zero. The normalization factor chosen was 8.

54

 

 

-an

100

 

 

 

-rnn

 

 

24

 

 

 

 

 

40L}. 1 r I ’ I I I
2 3 4
time (sec)

-1
.1

0

figure 34 : Original signal filtered and normalized using 5 of 100,8,and 4.

55

Determination of average number of zeros surrounding

each Nonzero Section

First the signal is filtered. Then the signal is normalized using the
integer division or normalization of the previous section. The normalization
will produce non-zero groups surrounded by zeros. We need to determine
how many zeros in a row will separate non-zero regions.What is desired is to
choose a number of consecutive zeros so that a single heartbeart would not be
split and consecutive heartbeats would not be joined. One or two zeros is too
few, as it was found that one or two consecutive zeros can occur in the
waveform of the same heartbeat. On the other hand there tend to be non-zero
regions close to the heartbeat waveform and ten or twenty zeros is too many.
How many zeros in a row should distinguish two non-zero groups, each
possibly a heartbeat, needed to be determined. A statistical analysis of
heartbeat data was used.

Using 84 segments with heartbeats and 33 with only noise, the average
number of consecutive zeros for heartbeat segments was 6.7 and the average
number of consecutive zeros for noise only was 4.9. Goodness of fit analysis
was performed to determine an adequate model. However, this did not lead to
any useful results. Detection is not performed at this level so the noise did
not come into consideration. It was found that for a signal with a reasonable
signal level,where the heartbeat could be detected, the number of consecutive
zeros in a heartbeat did not exceed four and the nmnber of consecutive zeros
between heartbeats was larger than this. Thus, five consecutive zeros was

chosen as the threshold for two non-zero segments to be considered distinct.

Detection and Estimation Algorithm

A good algorithm for determining if two waveshapes at different times
are close enough in some sense to be called consecutive heartbeats is needed.
When a observer looks at a collection of data, typically four seconds of data,
he is able to determine if a heartbeat is present, and if it is, where each
heartbeat is located, even when the signal to noise ratio is too low for
computer algorithms such as peak detection or correlation to work. The
observer uses pattern recognition for his detection and estimation system.
Without elaborate programming it would be very difficult for a computer to
take four seconds of data and do a similar procedure as the human. This is
because information easily interpreted from a picture is not available to the
computer. Information such as knowing if the heartbeats are there, where
they are, what is their period, what is their shape, and their amplitude. It is
also probable that the period, the shape and amplitude will change between
heartbeats. An observer looking at a graph can easily adjust to these changes,
but it is much harder for a computer working only with numbers. Following
this reasoning, the following algorithm was developed. Unless a good model
for the heart signal is determined, selecting an optimum algorithm cannot be
done. The frame of reference for judging algorithms, is how well it
performs compared to human recognition.

One basic premise is that if a human can look at a four second segment
of data and pick out the heartbeats, if present, an algorithm which can do the
same must exist. We may not be able to determine if a given algorithm is
optimum, but since we are using human recognition as the criterion,if it can

work as well or better than human recognition, it is an acceptable algorithm.

56

57
We want to have a process that can detect the repetition of similar shapes in
the segment. We cannot use typical correlation of the data because the
heartbeats are not periodic. We also cannot use a matched filter because the
shape of the heartbeat is not known and is continually changing. Because it is
a small region that is repeating, if we could take this region and correlate it
with the whole segment, as shown in figure 35, we would get rid of the

problem associated with non-periodicity.

1) I choose a sub-segment

    

 

2)

make a segment of this sub-segment and zeros

 

correlate the segments

 

 

 

 

figure 35

58

The problem is finding a subregion to correlate with the whole
segment. Sampling at 64 Hz we have 256 samples in four seconds. A typical
heartbeat is 8 samples. This means there are 248 regions that can be chosen as
a heartbeat. It is too time consuming to use each of the 248 sub segments and
correlate with the whole segment and see which has good correlation at
regular intervals in the segment. This is where the derivative filter plays a
role.The normalization factors of the filters can be set so few heartbeats will
be set to zero by normalization, and will produce zero regions. However, this
will not set all the noise regions to zero. So the filtering and normalization by
itself is not a good criterion for deciding if a heartbeat is present. However,
it will point to non-zero regions that could be chosen for the correlation.
Instead of having the algorithm look at the entire data segment to decide if
there are repeating pseudo periodic regions that are similar in shape, it will
only have to compare a smaller number of intervals produced by the non-
linear filter. Thus a combined system using the derivative filter to segment
the data into possible heartbeats, and a pattern recognition algorithm to get
rid of false detections will be a good detection and estimation scheme.

The output array from the filter contains zero elements separating

non-zero regions as shown in figure 36.

59

 

arrows point to these nonzero regions

2- .1 ‘1

 

J

 

 

«.100 . t . . . I .
0 1 2 3

time (sec)

figure 36 : non-zero segments after filtering

In this algorithm there must be five consecutive zeros to consider

regions distinct and different. For segments containing only heartbeat pulses,
typically the number of non-zero regions in a 4 second segment varies from
three or four,to fifteen or twenty. For noise signals, the number of non-zero

regions typically varies from one to more than twenty. Sample outputs are

shown in figure 37.

 

60

 

2O

10'I

 

-lOn

 

 

 

 

 

 

«at? I . , . I . .
0 1 2 3 4
time (sec)

figure 37 : examples of derivative filter output

The objective is determine if these non-zero regions are close enough
in shape, and are pseudo periodic in order to consider them heartbeats. One
way would be to use each sub-region to make a new segment as shown in
figure 35. The correlation of these regions would give an indication of the
period between regions and their similarity in shape. This procedure works
well, but it is also too time consuming because a new segment then
correlation is performed for each sub-region. The procedure chosen selects
one of the regions according to amplitude and derivative as being the most
likely to be a heartbeat and uses this region for correlation. Because the

derivative filter re-shapes the signal, and because the heartbeat waveform

61

usually extends beyond the non-zero regions determined by the filter, the

correlation is done using the original unfiltered signal.
40

 

20"
0-
.1

”20‘
2-40“

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

(108
o l 2 3 4
time (sec)

figure 382derivative filter output pointing to possible heartbeats

A 64 Hz sampling rate is high enough to capture the peaks and high
enough for the derivative filter, but is inadequate for correlation. Typically
there are three of four samples per peak and this does not give enough
discrimination between shapes. When doing the shape discrimination visually
using a graph, the points are interpolated and 64 Hz is an adequate sampling
rate. To get finer definition in the correlation, a higher sampling rate, 512
Hz, is used. Because because a 64 Hz sampling rate is high enough for good
peak definition, when doing the correlation on 4 seconds of data using the

512 Hz array (2048 samples), the region is slid over eight samples instead of

62
one after each step in the correlation.There are then 256 steps in the partial
correlation.2048/8=256.

An example of a 256 element output correlation file is shown in figure
39.

 

Se+4

4e+4 _,
2e+4 -
)e+o ‘
2w 1

1e+4 -1

 

 

0 1 2 3 4
time (sec)

figure 39 : correlation output

Using the region boundaries found from the derivative filter, the
highest peaks are found in the corresponding regions of the correlation file.
If there is no value greater than zero in the region, no value is taken for this

region. The final output file is all zero except at these highest points.

63

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

out) , - -- —--- -— --
o 1 2 3 4
time (sec)
figure 40: derivative filter output determining where to examine correlation
These points are then normalized, the largest becoming 100, using
integers, so that points smaller than 1/100 of the peak go to zero. An example

is shown in figure 41.

64

 

1e+2

1e+2-

8e+1-1

6e+1-

 

4e+1-I

26+l a

 

 

 

0400 T A 1
o 1 2 3
time (sec)
figure 41 : peaks of correlation in windows determined by derivative filter

One of the peaks in the correlation is the autocorrelation of the region
selected. This value is taken as a reference for determining likeness in shape.
Looking at i some threshold about this autoconelation value will give the
sub-regions that are similar in shape to the sub-region selected. If there are at
least two peaks within this threshold region, the algorithm will continue. If
there is only one peak, the sub-region selected is said to be dissimilar from
the other sub-regions and the algorithm proceeds to the next 4 second
segment. If there are only two peaks within the threshold region, the distance
or interval of these peaks is determined. If there are peaks separated by this
interval through the 4 second segment, even if the peaks are outside the
threshold region, the segment is said to contain heartbeats . If there are more
than two peaks in the threshold region, each inter-peak interval is
determined. Starting with the shortest interval, proceeding to the longer
ones, if there are peaks in the 4 second segment separated by this interval, the

segment is again said to contain heartbeats.

1e+2

 

“”2 ‘ threShold

8e+1-

6e+1 '1 ¢

/ other heartbeat peaks
4e+1 - J

2e+1-

 

 

 

 

 

 

 

0400 . r I |
0 1 2 3

time (sec)

figure 42 : correlation output showing threshold and heartbeat peaks

This file of peaks is passed to the final pattern recognition scheme, and
two features are used to determine if heartbeats are present, the amplitude of
the output peaks of the correlation and and their periodicity. The algorithm
followed is shown in figure 43.

66

 

Ecate peak positions I

J

determine which peaks are
within +the threshold of the
autocorrelation value (tpeaks)

 

 

 

 

 

 

 
   
   
 
  

obtain tpeak
intervals

  
   
 
 
 

intervals
almost equal

yes 0

 
     
 
    
   

  
   

  

lo the tpeaks exten

through the segment enough tpeaks

in interval

 

  
   
     
 

 

  
   

no
yes
reduce interval
check to see if peaks
are there but out of
threshold range
yes no
estimate
heartrate estimate no detection estimate
heartrate heartrate

 

 

 

 

 

figure 43 : algorithm used for pattern recognition

 

67

A flow chart of the complete algorithm is shown in figure 44.

 

ISample analog signal for 4 seconds. Get 2048 element
array.

 

 

 

Pe-sample to get 256 element array. 1

4 -

Pass small array though derivative filter and
eparate signal into segments.

1

hoose one segment as being most likely to be heartbeat
nd correlate it with the others using large array.

 

 

 

se output of correlation to determine if heatbeats are
resent. Determine heartrate if heartbeats present.

 

 

 

 

 

figure 44 : complete processing algorithm

Samples of complete processing are shown in figure 45 and figure

46. The top figure (a) is the original signal, (b) is the output of the derivative

filter, (c) the correlation output, and (d) the normalized peaks of the

correlation in the intervals determined by the derivative filter.

 

 

 

b) -n

 

 

5e+4
4e+4-
2e+4~
)e+0 a
2e+4 -
1e+4 .
c)..i‘

 

 

e+2
9+2 .
3e+1 -
3e+1 -1

1e”-1

2e+1-

 

 

61.... '

 

 

 

 

"NU |

1 i if

time (sec)
figure 45 : example of processing on human data

 

69

 

300

200 -

100 -1 n
q

- 100 j

-200 -
1

 

 

 

 

 

'300 -

a). arm .
20

 

 

 

 

 

 

 

b)-on
3‘6

2+6 -
3+6 :1 n

3+6 ~

 

 

2+0 -

3+6 - U

3+6 '1

C)-\¢‘ q
e+2

 

 

 

 

 

 

 

e+2'I
e+l -1

e+1-

 

6+1“

3+] -1 L
6)....) l
’7

 

 

 

QUU l I j I

time (sec)
figure 46 : example of processing on human data

RESULTS

This section compares this algorithm with new thresholds (l), and old
thresholds (2), with three others by having each classify heart data and noise
(for enumeration and explanation of the thresholds see the programs in
appendix A). One algorithm was developed by Albert Roseiro (3), the next is
simply autoconelation with a threshold (4), the last passes the data through a
median filter then uses the power spectral density for detection and
estimation (5). The results are shown below.

Eighty-four four second segments of heartbeat data and thirty five
segments of noise were classified by the five algorithms. In the table below,
correct detections means that the segments were classified as noise if noise,
and as having a heartbeat if a heartbeat was present. Correct estimations
means of the heartbeat segments correctly classified, the percentage of
correct estimations was calculated. As can be seen, there is a trade-off
between the percentage of correct detections and the percentage of correct
estimations. In general, the higher the correct detections, the lower the
correct estimations. This is expected because the lower number of correct
detections when a heartbeat is present means the algorithm perfonns

estimation only on relatively good segments.

70

71

correct detections (heart) correct estimations correct detections

algorithm (of those correctly detected) (noise)
1 43/84 (51.2%) 27/43 (62.8%) 30/35 (85.7%)
2 25/84 (29.8%) 18/25 (72%) 32./35 (91.4%)
3 20/84 (23.8%) 14/20 (70%) 32/35 (91.4%)
4 27/84 (32.1%) ‘ 11/27 (40.7%) 27/35 (77.1%)
5 65/84 (77.4%) 29/65 (44.6%) 22/35 (62.9%)

Deterministic signals were synthesized to compare the algorithms in a
easily controlled manner. Four different signal structures were used. The
first had one cycle of a sine wave repeated four times at equal intervals as
shown in figure 47. In actual human measurements the heartbeat pulse was
found to contain approximately 80 samples at 512 Hz sampling rate. The sine
wave was chosen to be 6.4 Hz to give 80 points for one cycle of the sine wave.
The peak of the sine wave pulse was 1000. The second structure placed the
four sine wave pulses at different intervals. One was placed to start at sample
2, the next at sample 600 the third at sample 1039 and the last at sample 1551
also shown in figure 47. The third structure was the same as number two,
except the sine wave pulses now had different frequencies.Gaussian noise of
known variance was the added to the first three signal segments. The variance
of the noise was increased until the algorithm could no longer give a conect
estimate. The last segment used was actual heartbeat data with a good signal
to noise ratio. The noise added was data from the microwave unit when this
was pointed at inanimate objects. This noise was then scaled to different

levels and added to the heartbeat data before the processing took place.

72

From the results in table 1 and 2, none of the four algorithms had
problems determining the pulse period of the first signal before noise was
added. After the noise was added, correlation and the median filter/power
spectral density did very well. Algorithm 1 is close but came in third.
Correlation does very well because the signal is perfectly periodic with
exactly the same signal repeating. The Median filter/PSD algorithm does well
because the median filter is very good at filtering gaussian noise, and because
the sine wave pulses are placed at equal intervals.

To test these assertions the second segment with unequal intervals
between sine pulses was used. As is seen it has no effect on the algorithm 1,
but has a large effect on the Median/PSD and correlation algorithms. The
correlation algorithm could not give the conect estimate even before noise
was added, so the noise level given is the largest noise variance possible for
the correlation algorithm to give the same estimate it gave before noise was
added.

To try and put algorithm 1 at a disadvantage, the third segment was
synthesized. The sine wave pulses now do not have the same period. Because
algorithm 1 uses a partial correlation, it is very sensitive to zero crossing
intervals. The algorithm was developed after noticing that even though the
period and amplitude of the heartbeat changes dramatically, the zero crossing
interval stays almost constant in actual measurements.

The last segment is an attempt to have a realistic comparison between
the algorithms. Roseiro's algorithm had done poorly with the previous signal
segments. The reason for it's usefulness is demonstrated here. It was, along
with algorithm 1, developed for the noise and heartbeat signals from the

microwave unit.

73

The results for the last signal segment are as expected. The results for
the other signal sections were very encouraging, indicating that the algorithm
presented may be used for a variety of signal shapes. This implies the
algorithm could work in many different environments, even if the signal
used for processing did not come from the microwave unit. The table below
shows the maximum variance of the noise that could be added to the
synthesized signal while the algorithm was still able to correctly estimate the
heartrate. The second table shows the values in the first table normalized to

one.

 

74

 

 

 

 

algorithm

 

 

 

 

 

 

 

 

 

 

 

 

algorithm

 

 

Table 1
segment
non-periodic heart
periodic non-periodic varying DUlSGS data
210,000 210,000 140,000 0.21
220,000 5,000 5,000 0.0
40,000 30,000 70,000 0.14
230,000 1 10,000 130,000 0.003
Table 2
segment
non-periodic heart
periodic non-periodic varying pulses data
0.9] 1.0 1.0 1.0
0.96 0.01 0.04 0
0.17 0.14 0.5 0.67
1.0 0.52 0.93 0.01

 

 

 

 

 

Algorithm 1 : Derivative Filter -Pattern Recognition.
Algorithm 3: Roseiro's Algorithm.

Algorithm 4: Autocorrelaton.

Algorithm 5: Median Filter-Power Spectral Density.

 

75

Below the four deterministic signals before noise was added are
shown.

1000

 

 

 

 

 

 

 

-1000 I I '7 v

 

 

 

1000
separation

corresponding to
53 BPM 70 BPM 60 BPM

 

 

0

 

 

-1000 . . . .

 

1000

 

01
is
I
N
\i
I
N

5.6 HZ

 

 

0

 

 

 

 

 

 

 

I

60

 

40‘
20"
o.
-20-
-4o~
-1060‘ r r r I .

 

 

 

0 1 2 3
time (sec)

figure 47 : the four synthesized signals with no noise added

76

 

 

 

2000
It
1000-
O -
.
1000~
W . 1 . . . A
1 2 3
thne(sec)

figure 48: segment 1 plus maximum noise for algorithm 1

CONCLUSIONS AND RECOMMENDATIONS

The algorithm described here is an attempt to have the computer
simulate the process humans go through when looking for heartbeats in the
microwave return signal. Because the return signal with heartbeats has
neither a known or constant shape and period, normal radar or similar
detection and estimation schemes will not work. This algorithm combines
non-linear filtering with partial correlation and pattern recognition. A four
second, 512 Hz sampling rate, 2048 element array is read in from the analog
board and re-sampled to 256 samples. This small array is passed through the
derivative filter to produce non-zero regions that give an indication of where
the heartbeats are. To determine if there are heartbeats within the non-zero
regions the following procedure is used. One of the regions is chosen as being
the most likely to be a heartbeat, based on derivatives and amplitude. Using
the large array, this region is correlated with the other regions and the output
is used as an indication of similarity in shape. If there are non-zero regions
that are similar in shape and are almost equidistant, these regions are said to
be heartbeats.

As seen from test results, simple correlation will not work well except
with a periodic heartbeat signal with a high signal to noise ratio. To help the
signal to noise ratio, a pre-filtering to accentuate the hearts impulsive
characteristic was developed. The idea for using integers to do the
calculations and nonnalizations, thus producing the zero regions, came from
Albert Roseiro. There is a large improvement using this pre-filter, and
judging from our experience trying many other pre-filters, not much

improvement can be expected in this area.

77

78

To be able to compare shapes adequately, the 512 Hz sampling rate had
to be used.In this algorithm one region is chosen and correlated with the
others to get an estimate of likeness of shape. The disadvantage of doing this
is the possibility of choosing a non-zero segment that is not a heartbeat when
heartbeats are present. Thus the major improvement in this area is a better
technique for selecting the region for correlation. Another possibility is to
use a lower sampling rate and having another method for measuring shape.
This could be done by finding an equation that interpolates the points of the
region, then compare equations of the non-zero regions. This may give a
better estimate of the shape, and it also eliminates the need to select one
region for use in correlation.

Only two features, the correlation output and periodicity are available
for the pattem recognition following the correlation. The algorithm would
be improved by finding more features. Many were tried, such as zero
crossing intervals, but no additional features were found. The next
improvement should be finding a better algorithm to determine if the output
correlation file contains heartbeats.

Another possibility is to combine this algorithm with the one
developed by Albert Roseiro.

As seen in the results presented,the algorithm developed in this report
works well under many conditions, sometimes even outperforming visual

inspection.

1)

2)

3)

4)

5)

6)

7)

8)
9)

REFERENCES

J .C. Lin, J. Kiemicki, M. Kiemiki and RB. Wollshaeger, "Microwave
apexcardiography," IEEE Trans. Microwave Theory Tech, vol.
MTT-27.PP.618-620, June 1979

G. Hoshal, M. Siegel and R. Zapp, "A Microwave heart monitor and
Life Detection System," Proceedings - Sixth annual Conference
IEEE Engineering in Medicine and Biology Society, Los Angeles,
pp. 331-333, September 1982

MA. Popovic, K. Chan and J .C. Lin, "Microprocessor-based
Noncontact Heartrate/Respiration Monitor," ibid., pp. 754-755

A.G. Favret and AF Caputo, "Evaluation of autoconelation
Techniques for detection of fetal electrocardiograms," IEEE Trans.
on Biomed. vol. BME-13, pp. 37-43, 1966

G. Hoshal, S. Ivkovich, M. Siegel, and R. Zapp, "Remote
Noncontacting Heart Rate Measurement Using Microwave Signals and
Digital Processing," Proceedings - Computers in Cardiology, Salt

Lake City, Utah, pp. 409-412 ,September 1984

W. Byme, P. Flynn, R. Zapp, and M. Siegel, "Adaptive Filter
Processing in Microwave Remote Heart Monitors," IEEE Trans. on
Biomed. vol. BME-33, no.7, pp. 717-722, July 1986

W. Byme, M. Siegel, "Adaptive Filter Signal Processing in
Microwave Heart Monitors,"

Personal communication with Greg Hoshal , Michigan State University

G. Hoshal, M. Siegel, "‘ "Advances on Identification of

Microwave Cardiograms for Remote Heart Rate Estimation",
Technical Report,MSU-ENGR-86-003, February 1986.

10)

11)

A. Nehorai, and B. Porat, " Adaptive Comb filtering for Harmonrc
Signal Enhancement," IEEE Trans. on Acoustics, Speech, and Signal
Processing, vol. ASSP-34 ,no. 5, October 1986, pp.1124-1138

L. Rabiner, S. Levinson, A. Rosenberg, J. Wilson, "Speaker-
Independent Recognition of Isolated Words Using Clustering
Techniques," IEEE Trans. on Acoustics, Speech, and Signal
Processing, vol. ASSP-27,no.4, August 1975, pp. 336-349

79

12)

13)

14)

15)

16)

17)

18)

19)

20)

80

R. Niederjohn, " A Mathematical Formulation and Comparison of
Zero- Crossing Techniques which have been Applied to Automatic
Speech Recognition," IEEE Trans. on Acoustics, Speech, and Signal
Processing, vol. ASSP-23 ,no. 4, August 1975, pp. 373-379

S. Soyegh, Y. Kok, and J. Hong, " An Algorithm to find two-
dirnensional Signals with Specified Zero Crossings," IEEE Trans. on
Acoustics, Speech, and Signal Processing, vol. ASSP-35 ,no. 1,
January 1987, pp.107-111

L. Rabiner, M. Cheng, A. Rosenberg, C. McGonegal, "A

Comparative Performance Study of Several Pitch Detection
Algorithms," IEEE Trans. on Acoustics, Speech, and Signal
Processing, vol. ASSP-24,no.5, October 1976, pp. 399-418

M. Sondhi, "New Methods of Pitch Extraction," IEEE Trans. Audio
Electroacoust., vol. AU-l6, pp. 262-266, June 1968

L. Rabiner, M. Sambur, C. Schmidt, " Applications of a non-linear
Smoothing Algorithm to Speech Processing" IEEE Trans. on
Acoustics, Speech, and Signal Processing, vol. ASSP-23 ,no.6,
December 1975, pp. 552-557

G. Arce, and M. McLaughlin, "Theoretical Analysis of the
Max/Median Filter," IEEE Trans. on Acoustics, Speech, and Signal
Processing, vol. ASSP-35 ,no. 1, January 1987, pp. 60-69

J. Astola, P. Heinonen, and Y. Neuvo, " On Root Structured of Median
and Med-type Filters," IEEE Trans. on Acoustics, Speech, and Signal
Processing, vol. ASSP-35 ,no. 8, August, pp. 1199-1201

V. Rao, and K. Rao " A New Algorithm for Real-Tirne Median
Filtering," IEEE Trans. on Acoustics, Speech, and Signal Processing,
vol. ASSP-34 ,no. 6, December 1986. PP.1674-1675

R- Gagliardi. Iatrodnctioutcflommunicanoufliugincering Wiley-
Interscience, New York, 1978

Appendix A

This is the main program. It accepts a 2048 point data
segment, a 4 second segment sampled at 512 Hz. This program
first removes the mean then sub-samples the 2048 point segment
to 256 points, corresponding to a 64 Hz sampling rate. This
segment then calls the processing subroutines for heartrate
detection and estimation in order.

#include <math.h>
#include <stdio.h>

#define St 512

#define seglen 2048

#def‘ine reduc 8

#define threshl 35 /*threshold for acceptance 01‘ \(llUCS

near autocorr*/
#define varthresh 9O /*thresh for variance*

#definc thresh3 10 /*threshold for accepting .rnml:
correlations*/

#define threshS 0.4 /*threshold for window in hiding
subroutine*/

#define thresh6 1.5 /*threshold for ends in hiding
subroutine‘v

#define offend 20 /*threshold for missing subroutine-*r

long data[2050],templ,der[260],w[2050],dat[260]: long llle'llL

int scgmentl,redseg;

int ibeg[260],iend[260],ilargc,numscg,autocorrsegmenu): 3:11'11‘11'02
FILE ‘ifpl,*fopen(),*ofp;

main()

float sr2,fheartrate;
int 5,563 l ,seg2,i,ns,sr l ,c,bit,j,i l,k,store.sr3.hea r t ra tes;
char fil_name[20],fils_name[20],store_l‘il[20].d.d I=""

printf("data file to process:");

scanf("%s",fils_na me);

if((ifpl=fopen(fils_name,"r"))==NULL) exiti l 1:

l‘scanf (if p l ,"%d,%f ,%d,%c%c %c,%c%c
%c\n",&ns,&sr2,&c,&d,&d,&d,&d,&d,&d1;

srlsns/seglen;

printf("\nnumber of samples=%d number of data bloets="ud'.n :.51-1 1;

printf("\nEnter starting and number of segments to he is n12" ;

scanf("%d,%d",&seg l ,&sch);

81

82

sr3=(int)sr2;
sr3=sr3/reduc;
redsegsseglen/reduc;
segmentl-se32*(redseg);
588681 *seglen; f or( 520; j<s;++ j){
fscanf(ifpl,"%ld\n”,&data[j]); }
for(s-0;s<segz;++s){ f or( jaO; j<seglen;++ j){
fscanf(ifpl,"%ld\n”,&data[j]); } printf("\nsegmentatt
%d",s+l); mean-0; for(i-0;i<seglen;++i){
meanamean+data[i]; } mean=mean ”seglen;
f or(i-0;i<seglen;++i){
data[i]-data[i]-mean; }
f or( j-O,k=0; j<redseg;++j,k=k+reduc){
dat[j]=data[k]; ) derivative(dat,der);
filter(dat,der); segments-segment(dat,dcr);
if(segments>49){goto out;} correlate(dat,der.data):
heartrates=pat(der,autocorr,thresh1); ifiheartrates<19 ||
heartrates>72){
f heartrate-0.0; } else{
fheartrate=(3840.0)/(float)heartrates: }
printf("\nheartratc=%d",(int)fheartrate1;
outz;
}
l‘close(ifp l );
}

83

This subroutine accepst a down-sampled version of :he 2018
length to 256 elements. The output is the derivative filtered
signal, normalized to +8.

filter(in,out)

long in[260],out[260];

{

int i,j,k;

long peak,ten=lO,hpeak,f1,f2;

for(i=- l ;i<redseg- l ;++i){
fl-in[i+l]-in[i];
f2-in[i]-in[i- l ];
w[i].f1-f2; ) w[0]=0; w[redseg-l]=0;

for(i-0;i<redseg;++i){
out[i]=in[i]*labs(out[i]);
w[i]-in[i]"labs(w[i]); } for(i=0;i<redscg-2;++i)1
if(w[i] =8 0 && w[i+l] == 0 && w[1+2] == 0'“,
out[i+l]=0;

} } peak-=0; for(i=0;i<redscg;++i){
if(peak<labs(out[i])){
peak=labs(out[i]);
} ) peak=(peak+4)/8; if(peak== ){pezil.=1:;
for(i=-0;i<redseg;++i){
out[i]=out[i]/peak; } for(i=l;i<redseg-l;++i‘1{
if(outli] is 0 && out[i-l] ==0 && outli+ll == 01-:
out[i]=0;
} }

}

derivative(in,out)

long in[260],out[260];

{

int i,j,k,oldi,ncwi;

long diff,oldpeak,newpeak;

newi=0;
newpeak=0;

for(i=0;i<redseg;++i){ out[i]=0;

)

for(i=l;i<redseg-l;++i){ if(in[i]>=in[i-l] && in[i]: :3:1.1i- l1,
Oldi=ncwi;
newi=i;

diff=newi-oldi;
oldpeakzncwpeak;

}

85

newpeak=in[i];

jeva(oldpeak,newpeak,diff,i,oldi1; } else
if(in[i]<-in[i-l] && in[i]<=in[i+l]){

oldianewi;

ncwi=i;

diff-newi-oldi;

oldpeaksnewpeak;

newpeak=in[i];

jeva(oldpeak,newpeak,diff,i.oldi); } else if: iniil==0
&& in[i-1]-=0 && in[i+l]==0){

oldi-newi;

ncwiai;

diff=newi-oldi;

oldpeaksnewpeak;

newpeakainfi];

jeva(oldpeak,newpeak,dii‘i‘,i.oldi‘1: :-

for(i=0;i<redseg-l;++i){ if(out[i] != out[i+l]){

}
}

out[i+l]=(out[i]+out[i+l]+l)/2; }

jeva(oldpeak,newpeak,diff,i 1 ,i2)
long oldpeak,newpeak,diff;

int

{

int

ii;

long slope,magdiff,hdiff;

magdiffalabs(newpeak-oldpeak);
hdiff=(diff+l)/2;
slopc=(magdiff+hdiff)/diff;
for(i=12;i<il;++i){ der[i]=510pe;

}
l

86

This subroutine is used in the nornmliznti...~. DIV-c":

a I if the input is positive, a -1 if it is negaizxe.

sign(in)
long in;

if(in<0){return(-l );}
else{return( l );}

( 1

1'1.“ 'll'll

88

This subroutine inputs the derivative filter output 1111(1
determines the regions of possible hear-theory 3111! 31-0-05 :5 the
region most likely to be a heartbeat.

int segment(dat,der)
long dat[260],der[260];
{

int i,j,k,iderl,ider2,ider3,corrlen,repeats;
long peak,peak2,peak3,datpeakl,datpeak2,datpcak3:
long diffl,diff2,diff3,avpeakl,avpeak2;

for(i=0;i<redseg;++i){
ibeg[i]-0;
iend[i]-O;

}

/* determine segment boundaries */

peak=0;
for(i-3,numseg-0;i<redseg;++i){
if(der[i-3]-s0 && der[i-2]-=O
&& der[i-l]==0 && ClCl‘[l] . n
ibeg[numseg]=i;
for(;i<redscg-3;++i){
if(der[i]!=0 && dcr[i+3]==0

&& der[i+2|==0 &.’\'- rli- '1. .7,
iend[numseg]=i;
++numscg;
goto out;
1
l
}
out:;
}

il‘(numseg>49){goto end;}

/"' determine segment used in correlation *.

peakao;

pcak2=0;

peak3=0;

for(i=0;i<numseg;++i){
corrlcn-iend[i]-ib¢8[il;
if(corrlcn>5 && corrlen<l7){

for(k-ibeg[i];k<-iend[i];++k){

if(peak<der[k])(
pea k=der[k];
iderlsi;

)

il'(pcak>dcr[k] && pcak2-:der[k| :

89

peak2=der[k];
iderZ-i;

)
if(peak>der[k] && peak2>dcr[kl .4 4' or"?- s’cr'T-‘ll

peak3=der[k];
ider3=i;
}
}
l
}
datpeak3=0;
datpeaklao;
datpeastO;
if(peak>0){
for(k=ibeg[iderl];k<=iend[iderl];++k)1
if(datpcakl<dat[k]){
datpeakl=dat[k];
}
)
}
if(peak>0){
for(k=ibeg[ider2];k<=iend[ider2];++k){
if(datpeakl>dat[k] && datpeak2<dat[k]>f
datpeak2=dat[k];
}
}
}
if(peak2>0){
for(k=ibcg[idcr3];k<=icnd[ider3];++kll
if(datpeakl>dat[k] && datpcakbdml»i , ~21 ;‘- ;ll;1t[k]){
datpeak3=dat[k];
}
}
}
ilarges-l;

avpeakl=(peak+peak2+peak3)/3;
avpcak2=(datpeakl+datpeak2+datpeak31 3:
if(avpeakl ==0 ){avpeakl=l;}

if(avpeakZ ==O ){avpeak2=l;}

diffl=(peak/avpcak l )+(datpcak l/avpea k: 1:
diff2=(pcak2/avpeakl)+(datpeak2/avpcak21;
diff3-(peak3/avpeak l )+(datpeak3/avpca R: 1:

if(diffl>=diff2 && diffl>=diff3){ilargc=idc1' l2;
else if(diff2>diffl && diff2>diff3){ilargc=ide1'1}
clsc if(diff3>diffl && diff3>diff2){ilarge=ide:‘E;2

if(ilarge<0){

goto end;

}

end: return(numseg);

}

9O

91

This subroutine inputs the original 2048 element segment and
correlates the segment choosen by the segment subrrmtinc with
all the other possible heartbeat regions. The liczrk 61' 1111
correlation in each possible segment is 11:: «nuns-a

correlate(dat,der,data)

long dat[260],der[260],data[2050];
1

int i,j,k,s,corrlen,corrlen2;
long peak,hpeak;

float fpeak;

peakao;

corrlen-iend[ilarge]"reduc-ibeg[ilarge]*rcd u 04-1'11'.‘ -
corrlen2=corrlen/(2*reduc);

for(i=0;i<corrlen2;++i){

dat[i]=Ol;

}

for(i=0,j=ibeg[ilarge]*reduc;i<corrlcn ;+ +1. + +1 1f
Wlil=data[i];

for(s=0,i=0;s<seglen-corrlen;++i,s+=8){
for(j-0,temp l =O;j<corrlen;++j){
temp l +-(data[j+s]*wll]);

}
da t[i+corrlen2]=temp l;
}
for(i=rcdseg-corrlen2;i<rcdscg;++ i ){
dat[i]=0|;
}
for(i=0;i<rcdscg;++i){
der[i]=0;
)
for(i=0;i<numscg;++i){
peak--100000;
for(kaibeg[i];k<iend[i];++k){
if(peak<dat[k]){
pea k=dat[k];
j=|<;
}
)
if(peak>0){
der[jl=peak;
if(i==ilarge){autocorr=j;}
}

fpeak=(float)der[autocorr];
fpcak=(fpcak)/100.0;

92

if(fpeak==0.0){fpeak=1.0;}
for(i=0;i<redscg;++i){

der[i]-(long)((float)der[i]/fpeak);
}

end:;

}

93

This subroutine inputs the peaks of the autocorrelation and
uses the peak amplitude and their separation to deternnnc il‘
heartbeats are present.

int pat(der,autocorr,peakthresh)
long der[260];
int autocorr,peakthresh;
{
int i,j,k,corr[50],diff[50],ihigh[50],numpeak.numbe:tt.~:
int avdiff,numdiff,smalldiff,ipeak,imedinn.threshltargedil‘l‘;
int oldipeak,sumbeats;
long x,var,peak;

if(heartrate-=O){heartrate=64;}
for(i=0;i<redseg;++i){
if(der[i]!=0){
printf("\ni=%d.dcr[i]=9’old".i.tler[ill:
}
}
beginz;
for(i=0,numpcak=0;i<rcdscg;++i){
if(dcr[i]>=(l00-pcakthresh) && dcr[i]<=! I'l'IHDCJl ' us 1‘: ;
ihigh[numpeak]=i;
++numpeak;

}

printf("\nnumpeak=%d”,numpeak);
if(numpcak>40 || numpcak<2){
heartrate=0;
goto end;

numdiffanumpeak-l;
l‘or(i=0,avdiff=0,smalldil‘f=lOOO.largtnlil‘I‘=-I- +73 :- i ii '- : ~ - '< ,;
diff[i]=ihigh[i+l]-ihigh[i];
avdiff+=diff[i];
if(smalldiff>diff[i]){smalldiff=dil‘l‘[i];}
if(largedil‘f<diff[i]){largediff=dil‘l‘[i];}

printf("\ndiff[i]=%d ihigh[i]=%d".dil‘f[i].ihighliIV

}
imedian=dif|‘[numdil‘I‘/2];
avdifl‘ /= (numdifl‘);
printf("\navdiff=°/ud smalldil‘l‘=”/ud".avdil‘l‘.sin:t-!
for(i=0,var=0;i<numdil'f;++i){
var+=((diff[i]-avdiff)*(dil‘l‘[i]—avdiI'l'H:

var/=numdiff;

94

printf("\nvar=%ld",var);
printf("\nimedian=%d",imedian);
if(var<varthresh){
goto missing;

)
else{
if(avdiff<30){
peakthresh-s3;
if(peakthresh<=0){
heartrate=0;
goto end;
}
goto begin;
}
goto hiding;
}
missingz;

printf("\nmissing");

if( (ihigh[0]-avdiff-ofl‘end)<0 &&
(ihigh[numdiff]+avdil’f-l-ol‘l‘end)>256 l{
heartratc=avdiff;
goto end;

}

else{
goto hiding;

)

hiding:;
printf("\nhiding");

i=autocorr-heartrate;
thresh4=heartrate‘threshS;
printf("\nautocorr=%d",autocorr);
numbeatSIO;
sumbeats=0;
oldipeak-autocorr;
while(i>(thresh4*thresh6)){
printf("\ni-%d heartrate=%d thresh4=%d".i,hea rtl‘11tc.th reshd)
for(j=i-thresh4,peak=0;j<i+thresh4;++j){
if(peak<der[j]){
pea k=der[j];
ipeakaj;
}

printf("\nl peak=%d",peak);
if(peak<thresh3){
goto down;

96
}
i=ipeak;
i-=heartrate;
sumbeats+=oldipeak-ipeak;
++numbeats;
oldipeakaipeak;
}
i-autocorr+heartrate;
oldipeak-autocorr;
while(i<255-(thresh4*thresh6)){

printf("\ni=%d heartrate=%d thresh4=%d",i.hea r: witty: 3 v -

for(j=i-thresh4,peak=0;j<i+thresh4;++_i l:
if(peak<der[j]){
peak=der[j];
ipeakaj;
}
}
printf("\nz peak=%d",peak);
if(peak<thresh3){
goto down;
}

i=ipeak;

i+=heartrate;
sumbeats+=ipeak-oldipeak:
++numbeats;
oldipeak=ipea k;

if(numbeats!=0){
heartratezsumbeats/nu mbea ts;

}
endz;
return(heartrate);
down;
printf("\ndown");
if(avdifl‘ != heartrate && nvdil'l' != ("M

hcartratc=avdifl‘;
goto hiding;

}

if(smalldiff != heartrate && smalldil‘l‘ != O){
heartrate=smalldiff;
avdiff=heartrate;
goto hiding;

}
if(largediff la heartrate && largedil‘l‘ != OH
heartrate=largedifl‘;
avdiff=hcartrate;
smalldiff=heartrate;
goto hiding;
}

else{

97

heartrate=0;
goto end;

98

This program inputs a segment of ASCII up to 512 samples and performs
a median filtering then a FFT on the data. The power specual density

is computed. The interpeak distances of the PSD l‘ilest re determined.

The interpeak distance ocurring with the most frequency is determined
to be the heartbeat period. The mainprogra m below in outs the ASCII
characters and calls the subroutines.

#include <stdio.h>
#include <math.h>

/* median filtering and power spectral density *9

float data[512],templ,mean;

float w[512],YYl[2][2048],YY2[2][2048];

float tp,psdmax,fest,startf,endf;

int segment],sr,seglen,window;

char fil_name[20],fils__name[20],store__l‘iI[20].d.cl i=“".

main()

{
calculate_wind0W();
media n();

psdl);

}

99

This subroutine performs the median filtering.

median()

{

int ‘fp,s,segl,seg2,i,ns,c.bit,one,j,i l,k,store:

double sqrt();

float xy,xi,x2,yi,y2,dat[256],rho,EY,EX,sigmax.'sigtn:t)‘.ft.l‘..'nax.min:
FILE *ifp,*ifpI,*fopen(),*ofp;

printf("data file to median filter:");

scanf("%s",fils_name);

if((ifpl=fopen(fils__name,"r"))==-NULL)
exit(l);

printf("\nEnter l to store filtered signal"):
scanf("%d”,store);
if(store == l){
printf(”Name of storing file2"):
scanf("°/os",store__fil);
if((ofp=fopcn(store__fil,"w"))==NU LL l
exit(l);

}
fscanf(ifpl,"%d,%d,%d,%c%c %c.”/oc‘?/nc ”mean".& ns.&sr.emtf..\'(1.8 d.&d.&d,&d,&d);
printf("\nEnter segment lengthz");
scanf("%d”,&seglen);
srsns/seglen;
printf("\nsegment Iength:%d",seglcn);
printf("\nnumber of samplcs=%d number of data blacks-”ltlflnxsr);
printf("\nEnter starting and number of segments In; I‘ .‘"l
scanf("%d,%d",&scgI,&scg2);
printf("\nEnter windowz");
scanf("°/od",&window);
segmentl=se32*(seglen);
printf("\nsegment l =%d",segmentl );
if(store =-- l){
fprintf(ofp,"%d,64,0,%cT %c,%cR ‘Voc\.n".segment l .d l.tl I .d 1. ll :
}

s=segl*scglen;
for(j=0;j<s;++j){
fsca nf(ifp l .""/uf\n".&data[j] l2
}
l'or(s=0;s<scg2;++s){
for(j=0;j<scglcn;++j){
fsca nf(ifp l ,"%f\n",&da ta[j] l:
}
mean=0.0;
for(i=0;i<seglen;++i){
mean=mean+data[i]:

}

100

mean=mean/seglen;
for(i-O;i<seglen;++i){

data[i]adata[i]-mean;
}

for(i=--window;i<0;++i){
data[i]=data[0];

for(i=seglen;i<(seglen+window);++i ){
data[i]=data[seglen-I];
}
for(i-O;i<seglen;++i){
for(j-O;j<2*window+ l ;++j){
dat[j]=data[j+i-windOWI;
I

for(j=0;j<2"‘window;++j){
for(k=j+ l ;I<<2"window+ l ;++k ){
if(dat[k] < dat[j])(
temp l =dat[j];
dat[j]=da t[k];
dat[k]=tcmp I;

}
} }
YYl[l][i+l]=dat[window];
}
}
felose(ifpl);
if(store == l){
felose(ofp);
)

}

101

This subroutine calculates a hanning window used in 'hc FF'I‘

calculate_window()
{
double cos(),zarg;
float sk,sc,ti;
register int i,nsoff,iend;
sk - 0.5;
nsoff - seglen / 2.0;
i = 0;
iend . seglen / 2.0;
whilc(i < iend){
zarg . tp ‘ i / seglen;
sea sk + sk * cos(zarg);
w[nsoff+ l + i] a sc;
w[nsoff - i] - sc;
i++;
}
w[l] = 0.0;

102

This suproutine caclulates the PSD

psd()
{
tp = 6.2831853;
windows();
f f t();
magnitude();
}
windows()
I
register int i;
i =- 0;
while(i < seglen){
i++;

YYl[l][i] = YYl[l][i] * w[i];
YYI[2][i] = 0.0;
YY2[I][i] = 0.0;
YY2[2][i] = 0.0;

103

This suproutine caclulates the FFT

fft()
{
int p,k,im,pl l,p22,p33,id,ic,inh,ia,ib;
int il;
float bi,br,ti,zi,zr;
double zsr,q,log2,log(),sin(),cos(),two;
two=2.00; '
log2=log(two);
zsr - seglen;
q . log(zsr);
zsr a q/IogZ;
ti = tp / seglen;

iI= zsr;
ia= seglen / 2;
ib - I;
inh=ia;
p = 0;
while(p < il){
p++;
ic -- O;
id=ia;
k = 0;
while(k < ib){
k++;
q a ti * ic;
zr = cos(q);
zi = - sin(q);
im=ic;
while(im < id){
im++;
pl I=im + id;
p22=im + ic:
p33=im + inh;
br= zr * YYI[|][pl l] - 4i “t ‘: I-lepllf:
bi = zr * YYI[2][pI I] + at [Ill-3! l !:
YY2[I][im] = \’Yl[l][p22j + I»:
YY2[2][im] = YYI[2][p22] + bi;
YY2[l][p33] = YYl[l][p.‘.2 - 1‘3:
YY2[2][D33] = YYI[2][p.2:! - I. i;
}
ic=id;
id=id + ia:
}
ia=ia / 2;
ib = 2 * ib;
im = 0;

while(im < seglen){

104

im++;
YYI[I][im] = YY2[I][im]:
YYI[2][im] == YY2[2][im]:

105

This suproutine is used by the PSD subroutine

magnitude()

{

register int i;
float sk;
FILE ‘of p;

ska 2.0 / (sr " seglen);
i :- 0;

printf("\nName of storing psd filez"):

scanf("%s",store_fil);

if((ofp-fopen(store_fil,"w"))==NIJ L L )
exit(l);

fprintf(0fp,"%d,64,0,°/0CT %C,%CR ‘VOC\n",scgmcnt I .tl f.d l .d .d l );

while(i < seglen){
i++;

YY2[l][i] - sk * (YYI[I][i] * YYI[I][i] + x‘s-1321.3 mum);
fprintf(ofp,""/of\n",YY2[l][i]);

}

felose(ofp);

106

This suproutine performs the peak search of the PSI) .‘lll es

pcak_search()

{
register int i,imax,il,i2;
float psdmax,fest;
il - startf "‘ seglen / sr;
i2 =- endf "' seglen / sr;
printf("\nstart,finish= %d,%d\n",i l,i2);
psdmax . 0.0;
imam-0;
i il-I;
while(i<i2){
i++;
if(YY2[l][i] > psdmaxll
psdmax = YY2[l][i];
imax a i-l;
}
}

fest = (float)imax * (float)sr / (float)seglen;

}

107

This subroutine performs a correlation on an input ASCII segment
and compares the first peak after zero to a threshold. If the peak
is larger than the threshold, the distance from zero to the first
peak is taken as the hearbeat period.

#include <stdio.h>

float dat[260],r[260],threshold;

int sr,ns,lmin,lmax,lrmax,lag,nseg,segnum,numseg,locat.e;

int wsize,cwsize,bmin,bmax,imax,rect,izr,ss,f f ,sum:

double sqrt(),var,fabs(),rvar,bvar;

float hrest,bavg,ravg;

char fil__name[20],d,store_fil[20],d12"";

float suml,sum2,sum3,mean,varl,var2,s,rmax,mean2.taper.scale;
FILE *ifp,*fopen(),*ofp;

main()

I
int i,j,k;

bavg=0.0;
bvar=0.0;
war-0.0;
ravg=0.0;
i=0;
izr=0;
wsize=256;
printf("file to correlate:");
scanf("%s",fil_namc);
if((ifpsfopen(fil_name,"r"))==NULL)
exit(l);
fscanf(ifp,"%d,%d,%d,%c%c %c,%c%c %c\n".&ns.&sr.&e.&d.&tl..&*d.&d,&d,&d);
printf("# of samples=°/od # of 256 blocks=%d'\n".ns.ns 'ij;
printf("\nnumber of blocks:");
sca nf("%d,%d”,&nseg);
sum=nseg*256;
printf("Full wave rectify?(I=yes,0=no):");
scanf(”%d",&rect);
printf("rect=%d\n",rect);
lmin=0;
lmax=256;
segnum=0;
mean=0.0;
locat=0;
while(segnum<nseg){
printf("\nsegment#=%d",segnum);
i=0;

108

while (i<256){
i=i+l;
fscanf(ifp,"%f\n",&dat[i]);

}
corr();
i=0;
rmax-O;
i=lmin+l;
segnum-segnum+l;

}
felose(ifp);

109

This subroutine performs the autocorrelation of the input file and
compares it to a threshold.

corr()
{
int i,j,k;
float max,dmax,databs;
i=0;
mean-0.0;
while(i<wsizc){
++i;
mean+=dat[i];
}
mean-mean/(float)wsize;
i=0;
suml=0.0;
sum2=0.0;
mean2=0.0;
while(i<wsizc){
i=i+l;
dat[i]=dat[i]-mea n;
databs=fabs(dat[i]);
if(dmax<databs){
dmax=databs;
}
if(rect>izr){
dat[i]=databs;
mean2=mean2+dat[i];
}
}
mean2=mean2/(float)wsize;
if(rect>izr){
i=0;
while(i<wsize){
i=i+l;
dat[i]=dat[i]-mean2:
}
}
cwsize=wsize-lmin;
i=0;
Iag=lmin;
while(lag<lmax+ I ){
sum3=0;
i=0;

j=wsize-lag;
while(i<j){

110

i=i+l;
sum3=sun13+dat[i]*da t[i+lag}:
} .
S=J;
s=s"'s;’
r[Iag+I]=sum3;
lag=lag+l;
}

max=0;
for(i-30;i<-lmax;++i){
if(max<r[i]){
max=r[i];
imax=i;
}
}
printf("\nmax=%f dmax=%f",max.dmax);
thresholdamax/dmax;

if(threshold> I 600){
printf(”\nheartrate=%d",3840/ima x );

else{
printf("\nheartratc=0");

}

"’lllllllllllllllllllf