# Automatic Detection of Signalling Behaviour from Assistance Dogs as they
Forecast the Onset of Epileptic Seizures in Humans

Hitesh Raju1, Ankit Sharma1, Aoife Smeaton2, Alan F. Smeaton1,3 1School of
Computing, Dublin City University, Ireland
2Dogs for the Disabled, Cork, Ireland
3Insight Centre for Data Analytics, Dublin City University, Ireland
<EMAIL_ADDRESS>
(2023)

###### Abstract.

Epilepsy or the occurrence of epileptic seizures, is one of the world’s most
well-known neurological disorders affecting millions of people. Seizures
mostly occur due to non-coordinated electrical discharges in the human brain
and may cause damage, including collapse and loss of consciousness. If the
onset of a seizure can be forecast then the subject can be placed into a safe
environment or position so that self-injury as a result of a collapse can be
minimised. However there are no definitive methods to predict seizures in an
everyday, uncontrolled environment. Previous studies have shown that pet dogs
have the ability to detect the onset of an epileptic seizure by scenting the
characteristic volatile organic compounds exuded through the skin by a subject
prior a seizure occurring and there are cases where assistance dogs, trained
to scent the onset of a seizure, can signal this to their owner/trainer. In
this work we identify how we can automatically detect the signalling
behaviours of trained assistance dogs and use this to alert their owner. Using
data from an accelerometer worn on the collar of a dog we describe how we
gathered movement data from 11 trained dogs for a total of 107 days as they
exhibited signalling behaviour on command. We present the machine learning
techniques used to accurately detect signalling from routine dog behaviour.
This work is a step towards automatic alerting of the likely onset of an
epileptic seizure from the signalling behaviour of a trained assistance dog.

Epileptic seizure, seizure-alert dogs, animal signalling behaviour, wearable
accelerometer, feature selection, machine learning

††journalyear: 2023††copyright: rightsretained††conference: The 38th
ACM/SIGAPP Symposium on Applied Computing; March 27-March 31, 2023; Tallinn,
Estonia††booktitle: The 38th ACM/SIGAPP Symposium on Applied Computing (SAC
’23), March 27-March 31, 2023, Tallinn, Estonia††doi:
10.1145/3555776.3577656††isbn: 978-1-4503-9517-5/23/03††article: 4††price:
15.00††ccs: Applied computing Life and medical sciences††ccs: Hardware Sensor
applications and deployments††ccs: Computing methodologies Supervised learning

##  1\. Introduction

Epilepsy, the occurrence of epileptic seizures, is one of the world’s most
common neurological disorders affecting millions of people daily. Seizures
occur due to non-coordinated electrical discharge in the brain and may cause
damage including collapse and loss of consciousness or worse. If a the onset
of a seizure can be detected then the subject can be placed into a safe
environment or positioned so the risk of self-injury as a result of a collapse
can be minimised. However there are no reliable and practical methods to
predict the onset of seizures in an everyday, uncontrolled and routine
environment.

Previous studies have shown that dogs, either household pets or trained
assistance dogs, have the ability to forecast the onset of an epileptic
seizure by detecting the scent of a particular kind of volatile organic
compound associated with seizure onset (Maa et al., 2021). Researchers have
also demonstrated the presence of seizure-associated odours secreted from the
skin and then recording the reactions and signalling behaviour of a pet dog
when it senses such odours.

In the present study we investigate how to automatically detect the signalling
behaviour of trained assistance dogs who detect a seizure onset and who
showcase this through their signalling. For this we use a wearable sensor
including an accelerometer on the collar of the dog. We develop a data
processing pipeline to pre-process the sensor data and use this as input to a
machine learning algorithm that recognises the assistance dog’s signalling
behaviour. The accuracy of the detection of the signalling behaviour is such
that this one step in the overall goal of using assistance dogs in practice,
is demonstrated as being feasible.

The rest of this paper is organised as follows. In the next section we present
some background and literature including the challenges in this work in
Section 2.6. We then give an overview of our data processing pipeline followed
by a description of the sensor data we gathered from trained dogs, including
the annotation of this as data for training a machine learning algorithm. In
Section 4 we describe our data pre-processing including feature extraction,
data analysis and the algorithms we used for classification and prediction. In
Section 5 we present our results including a discussion before our final
conclusions.

##  2\. Background and Literature Review

###  2.1. Epilepsy

Epilepsy is a neurological brain disorder affecting millions of people daily.
It is characterised by the occurrence of spontaneous and seizures where the
subject loses control of motion and may jerk and thrash around uncontrollably.
In simple terms a seizure can be defined as an internal electrical storm in
the brain. It is the consequence of abnormal, excessive discharges by nerve
cells. The seizures in epilepsy are mostly related to a brain injury or are
hereditary, but more often the cause is completely unknown (Powell et al.,
2021). Seizures result from a non-coordinated neuronal electrical discharge
and recurring seizures occur in the cerebral cortex (Rohan et al., 2020) and
may cause serious damage, including loss of consciousness and collapses and
self-inflicted injury where the subject can bang their head on the ground
(Litt and Echauz, 2002), for example.

Some of the common causes of epilepsy are inborn or developmental including
genetic abnormalities or structural problems in the brain like malformed veins
or areas which have not developed normally. One in three people will be
diagnosed with epilepsy by the age of 60. The unexpected and uncontrolled
seizures not only make subjects feel vulnerable but also impact their quality
of life. While there are medical interventions which can reduce the
seriousness or the number of seizures, many of these pharmacological
interventions have undesired side-effects and getting the dosage as well as
the combination of medicines correct, is difficult.

###  2.2. Seizure Onset Prediction

The quality of life of those who experience epileptic seizures would be
significantly improved if the onset of seizures could be predicted in advance.
In such a case, a warning mechanism for the onset of a seizure would allow
subjects time to remove themselves from a potentially dangerous situation,
e.g. from a staircase, to somewhere safer like lying on a bed, a sofa or even
lying on the floor. Furthermore, seizure prediction offers the possibility of
intervening in the dynamics of the brain before or at the start of a seizure
and ultimately aborting seizures before they occur (Yu et al., 2016).

Although research has been carried out on predicting epileptic seizures since
the 1970s, there are still no reliable or practical methods available to
predict seizure onset in subjects with epilepsy. Therefore there is a pressing
need to develop approaches to seizure onset detection. The introduction of any
technique to detect and accurately predict the onset of a seizure, giving
sufficient time for a person to take control of the situation would also
encourage a healthier and more normal lifestyle as the risk of self-industry
during a seizure would be reduced. Since seizures are unpredictable, those who
are susceptible to seizures must live a lifestyle which avoids activities
which contain any serious kind of risk, certainly when they are operating
without a carer or assistant (Bhuvaneswari et al., 2020).

###  2.3. Assistance Dogs for Predicting Seizure Onset

Humans have studied animal behaviour for hundreds of years including during
domestication (Naik et al., 2020). While most dogs in domestic settings are
pets and have evolved to naturally live with human families offering
companionship, assistance dogs are trained to undertake a variety of tasks to
help individuals. Sometimes these individuals have physical or even cognitive
disabilities that impair their quality of life.

One form of assistance that assistance dogs can offer is anticipating the
onset of seizures. Early work in this area was reported more than 2 decades
ago, work which analyse the reaction of dogs to odours and that work found
that the dogs could predict when a seizure was imminent based on smell alone
(Litt and Echauz, 2002). This early work was followed by work which proved
that epileptic seizures are associated with a specific kind of odour which
comes from a naturally occurring volatile organic compound that we secrete
through our skin, and which is detectable by pet dogs (Brown and Goldstein,
2011). More recent work by Powell et al. found that this scent which is an
indicator of seizure onset can be detected even by untrained pet dogs (Powell
et al., 2021) while work in (Maa et al., 2021) isolated a unique seizure scent
which acts as a biomarker for seizure onset indication.

In addition to the potential for scent-based detection of seizure onset,
visual cues from a subject’s behaviour or postural changes (e.g. changes in
motor activity and or mood, or experiencing of auras) prior to an impending
seizure can also be a form of signalling from a subject which could be
detected by an assistance dog (Catala et al., 2018). Whether scent-based or
behaviourally based, this past work indicates that assistance dogs could
effectively be trained to recognise a seizure onset and then give an overt
signal by their behaviour, to warn his/her owner or caregiver.

###  2.4. Wearable Accelerometers on Animals

The use of wearable accelerometers by people to measure physical activity,
sedentary time, energy expenditure and sleep-related behaviours, is a well-
established market and the billion-dollar pet industry is now getting involved
in this wearables market. The global pet wearable market is expected to grow
at a compound annual growth rate of 13.5% before 2025, with enormous growth in
demand from Asia after initial popularity in the West. The typical pet
wearables on the market are dog activity and fitness trackers, like FitBark,
PetPace, and Whistle. These devices employ accelerometers, GPS, and vitality
sensors to measure the activity and sleep patterns of pets (Zamansky et al.,
2019).

Accelerometers are used in animal science in various contexts such as
estimation of energy expenditure and assessment of behaviours in wildlife
(Kumpulainen et al., 2021). Existing activity trackers are useful for
evaluating simple canine behaviours according to past research. Depending on
the sensitivity of the measurement unit, the spontaneous activities of a dog,
such as locomotion, postural change and movement of body in each posture, can
be differentiated by analysis of accelerometer data (Yam et al., 2011).
Accelerometer data can be used to extract more detailed behaviours to assess
an animal’s health and wellbeing. Previously these devices were used to detect
musculoskeletal problems, ear infections and mental health conditions such as
separation distress (den Uijl et al., 2017). A number of systematic reviews of
how to process data from wearables showcases that various attributes must be
considered for accurate data collection and processing. These include
placement of the sensor, sampling frequency, filters used on the gathered
data, epoch length and non-wear-time (Migueles et al., 2017; Riaboff et al.,
2022).

###  2.5. Machine Learning on Accelerometer Data

The accelerometer data from collar-mounted canine activity monitors can be
used to measure a dog’s activity levels including step count and distance
travelled. Recent advances in machine learning and embedded computing have
made it possible to classify animal behaviour in a more nuanced and accurate
manner. Through a collar-mounted accelerometer, one study confirmed that
activity monitors using validated algorithms can accurately detect important
canine behaviours (Chambers et al., 2021). Validated algorithms have wide
practical applications when used in commercially available canine activity
monitors. Based on the accelerometer data, algorithms accurately detected
eating and drinking behaviour, other behaviours such as licking, petting,
rubbing, scratching, and sniffing were also accurately detected (Chambers et
al., 2021).

In this work we set out to collect information from a number of volunteers
participants and their dogs on signalling behaviours or “tricks” that their
dogs can perform while wearing an accelerometer to collect data on their
movement. We also gather information on the various characteristics of the
dogs including age, size, breed and the range of signalling behaviours it can
perform. Video footage of the dogs performing those behaviours will also be
recorded for short instances several times during the data logging period.
When we manually annotate this video footage for the occurrence of such
behaviours it gives us a ground truth of positive instances. We can then
combine the behaviour detection results across a variety of dog breeds and
determine the level of accuracy for automatic detection of seizure onset based
on the analysis of accelerometer data from accelerometers worn on the collars
of trained assistance dogs.

###  2.6. Challenges

In order to investigate how to detect canine signalling behaviour we collect
data from an always-on accelerometer combined with videographic recordings of
dogs while performing their signalling behaviour. We use this as training data
for a machine learning algorithm.

Using a triaxial accelerometer comes with hidden challenges. One problem is
that accelerometer data may not be directly comparable across dogs. The shape
and size of the dog that the device is attached to raises issues such as its
relative height with respect to the ground which can alter its relative
position and result in data being dissimilar across different dogs. Thus the
dogs’ physical characteristics which include size, age, weight and breed are
important. Dogs which are physically strong and well-built may not be as agile
as lighter breeds of the same size who might perform certain signalling
activities more rapidly.

Previous studies have indicated that orientation inconsistency of a wearable
sensor is important (Majumder and Deen, 2020). For placement on a dog, the
best placement for an activity sensor has been determined to be ventral
attachment to the neck collar because this placement also makes it possible to
detect behaviours that do not involve movement of the whole body, such as
scratching or eating (Hussain et al., 2022). However, when worn by participant
dogs on their collars, the collar and thus the sensor can rotate around the
neck, or be attached initially in different orientations. In some scenarios
accelerometers may be mounted on a rigid structure such as a fixed harness
worn by the dog in order to maintain a standard sensor orientation. While this
addresses the issue of sensor orientation, wearing a harness 24/7 is
uncomfortable for an assistance dog and this is not an option for our eventual
use case so in this work we fix the sensor on the dog collar even if the
collar may rotate around the dog’s neck.

To isolate any bias from inconsistency of the sensor orientation we compute
the average sum of the three axial values. Thus, for the purpose of our work
we treat rotation of the sensor and the collar as systematic biases that are
associated with the devices. A more in-depth and principled analysis of this
assumption is beyond the scope of the current study and the reader is referred
to (Riaboff et al., 2022) for further discussion on this topic.

##  3\. Data Gathering and Annotation

Data for this work was gathered from a group of 11 volunteer dogs with an age
ranging from 1 to 9 years. All dogs were deemed to be agile and healthy by
their owners and this study was approved by the Research Ethics committee of
Dublin City University School of Computing. All volunteer dog owners read a
plain language statement and signed and returned an informed consent form.
Volunteers were recruited by information leaflets circulated among dog
training social media groups and almost all volunteer owners were themselves
dog trainers. Information on each volunteer dog was documented including name,
the tricks or signalling behaviours that it could perform on command and we
purposely included a range of breeds, and ages among the volunteer dogs as
shown in Table 1. The unprocessed sensor data is publicly available at
(Smeaton et al., 2022). In screening volunteer dogs for participation, the
ability to perform certain agile activities such as spins and jumps were
mandatory and dog owners were informed of this prior to recruitment.

Table 1. Characteristics of Volunteer Dogs Number | Breed | Size | Age | Data logging | Behaviours | No of videos
---|---|---|---|---|---|---
1 | Collie Mix | Medium | 3 Yrs | 9 days | Spins, Jumps and Rollovers | 20
2 | Springer/Terrier MIX | Medium-Small | 8 Yrs | 4 days | Spins, Jumps | 16
3 | Spaniel Mix | Medium | 4 Yrs | 9 days | Spins, Jumps | 27
4 | Poodle | Medium-Small | 2 Yrs | 9 days | Spins, Jumps | 9
5 | Lurcher | Medium | 5 Yrs | 3 days | Spins, Jumps | 14
6 | Lurcher | Medium | 2 Yrs | 9 days | Spins, Jumps | 9
7 | Collie Mix | Medium | 3 Yrs | 8 days | Spins, Jumps | 8
8 | Hungarian Vizla | Medium | 3 Yrs | 12 days | Spins, Jumps | 6
9 | Border Collie | Medium | 2 Yrs 4 Months | 19 days | Spins, Jumps | 11
10 | Australian Retriever | Medium | 1 Yr | 14 days | Spins, Jumps | 9
11 | Golden Retriever | Medium | 7 Months | 11 days | Spins, Jumps | 3
|  |  | Average age 3.5 Yrs | Total 107 days |  | Total 132

The AX3 accelerometer (Axivity, York, UK https://axivity.com/) is a popular
device for wearable sensing research. A pre-configured and full-charged device
with parameters preset for data gathering as described below, was sent to each
volunteer and placed securely on the dog’s collar and left in place to gather
data in the dog’s natural environment and routine behaviour. The AX3 data
logger, shown in Figure 1, is a wearable sensor which includes a microelectro-
mechanical systems (MEMS) accelerometer and flash memory. Multidimensional
data collection and configurable resolution/frequency are two of the reasons
why the AX3 sensor was chosen. This allowed the sensor to be set at a sampling
rate up to 100Hz if needed.

Based on literature describing the use of the AX3 in other animal studies, a
sampling rate of 12.5Hz and 12 bits for data resolution was considered optimal
for our experiment (Riaboff et al., 2022), trading battery live vs. sensing
frequency. It is possible to represent sensor data captured as either single-
dimensional data or multi-dimensional data which constitutes the x, y, and z
directional dimensions. The energy consumption for the wearable device is an
important parameter due to the fact that dog activity should be tracked over a
prolonged period of time. Reducing sample rates, logging fewer accelerometer
axes, or reducing features to record can increase battery life and for our
configuration, battery life for full charges is between 6 and 8 weeks. A
recent study examined dairy cows’ laying, standing, and feeding behaviour and
found that fewer logging axes did not affect the accuracy of the
categorisation (Vázquez Diosdado et al., 2015). By choosing optimal sensor
measurement settings such as sampling rate, logging axis, and feature
selection, the energy consumption in a wearable device can be reduced
significantly. A number of studies have examined these features of
accelerometer data for human activity classification, but none have examined
how they affect behaviour recognition accuracy in dog tracking.

Figure 1. AX3 Accelerometer

Each dog performs a range of different activities as part of their natural
behaviour during the data logging period including standing, sitting,
sleeping, lying, jumping and running, and all the this while wearing an
accelerometer on their collar. Occasionally during the data logging period and
at a time of her/his choice, the owner would record a video on a smartphone of
their dog performing a range of signalling behaviours on command. A total of
132 such videos were shared with us through an online platform as soon as the
video recordings were made. This allowed us to retrospectively identify and
synchronise the video recordings with the accelerometer data. When the logging
period was completed and the sensor returned to us, data was downloaded and
the device was cleaned and shared with the next participant to repeat the
process. Figure 2 outlines the data flow of the data gathering process.

Figure 2. Data flow during data gathering

Table 2 shows the most common of the video-recorded dog behaviours in our
investigation across all dogs indicating the activity performed and the time
to complete it as averaged across all participant dogs. Not all dogs were able
to complete all behaviours which meant we had to decide on a single signalling
behaviour which was

  1. (1)

common across many/most/all assistance dogs or that an assistance dog could be
trained to do;

  2. (2)

rarely occurring in a dog’s natural behaviour, so not sitting for example as
that is a fairly commonplace activity;

  3. (3)

distinct enough in terms of motion characteristics to make it more easily
distinguishable from other dog movements.

Table 2. Activity Table Number | Activity type | Average completion time (all dogs)
---|---|---
1 | Spin clockwise | 1.5 sec
2 | Spin anti-clockwise | 1.2 sec
3 | Stand | 2.2 sec
4 | Jump | 1 sec
5 | Sit | 1 sec
6 | Rollover | 2 sec
7 | Idle | N/A
8 | Other Movements | N/A

Behaviours observed during the recorded sessions across all dogs included
clockwise spins, anti-clockwise spins, sitting, standing, jumping and
rollovers as indicated in Figure 3. From this we can see that spins, where a
dog spins in a circle with all four legs on the ground, is a signalling
behaviour that meets our requirements of being easy to train, occurs rarely in
day-to-day activities and is distinct enough from other dog movements to be
automatically detected. There is also a sufficient number of example
recordings in the videos to provide enough data on which to train a machine
learning algorithm to recognise this behaviour automatically.

Figure 3. Numbers for each signalling behaviour manually annotated as training
data in volunteer’s videos

The data recorded on the AX3 accelerometers was downloaded at the end of each
particiant’s logging period when the sensor was returned to us and the raw
CWA-file was exported for further processing. Data was validated to confirm
that the actual sample rate of the accelerometers matched the sample rate
predefined in the OMGUI software
(https://github.com/digitalinteraction/openmovement/wiki/AX3-GUI) used to
configure and manage the AX3. The recorded data was then trimmed to remove
data logging while the AX3 was in transit from us to the dog (by post) and
from the participant back to us (by post) and cases where it was idle for
prolonged periods. The volunteers had been briefed prior to sensor placement
about the range of behaviours that should be performed and recorded, and the
raw sensor data was manually inspected in order to annotate it with a
corresponding activity label by viewing the synchronously recorded videos.
This was also the method used to annotate sensor data.

As part of the process of annotating video data, activity definitions were set
before analysing the videos. The accelerometer signal was not visible to the
annotators while annotating. Video recordings were viewed locally and manually
annotated frame-by-frame by two observers blinded to the sensor data and using
a controlled vocabulary. To calculate the annotation format, the initial and
final times of each activity were recorded. All the annotated information was
saved which included the number and labels of movements and their time
intervals, each label corresponding to a single dog behaviour. Sensor data was
imported into CSV format and synchronised with the videos for each dog. Once
annotated, data from each dog’s recording session was exported into a separate
data file containing features for time, sensor data, and the annotated
signalling behaviour.

Several important points were identified during the annotation stage where
sensor orientation inconsistency was observed and noted. Several examples of
orientation inconsistencies were identified. In Figure 4 for example we see
spikes in the intensity of motion in the y-axis and to a much lesser extent in
the x- and z-axes of the accelerometer (blue and red lines respectively in
Figure 4) while the jump indications were observed to a much higher extent in
the accelerometer’s y-axis (green line in Figure 4. Yet for the same dog but
recorded a day earlier, the recordings for the same jumping behaviour is shown
in Figure 5 indicating the detection of the jumping motion by the x-axis (blue
line in Figure 5) and inversely by the z-axis (red line in Figure 5) as well
as by the y-axis (green line in Figure 5). This indicates that for the same
jumping behaviour by the same dog and of approximately the same intensity
(through the graphs in Figures 4 and 5 are not to the same scale), the dog’s
collar must have rotated around her neck and hence the sensor’s orientation
had changed. This demonstrates the challenge of managing the inconsistency of
sensor orientation.

Figure 4. Jump indication on the y-Axis

Figure 5. Jump indication on the x-axis and z-axis

Finally the common annotations from the two observers were merged and exported
to a single file. This method was followed by annotating for behaviour for
each dog individually, by each annotator independently. Where annotators
disagreed, i.e. where one annotator indicated a signalling behaviour and the
other did not, then that annotation was discarded and it was only the
commonly-agreed annotations that were retained. Once the annotation process
was completed, the frame-by-frame output and time labels were exported to
another CSV file. Finally this data could be used in various combinations to
carry out the investigations. The number of annotations for each of the major
dog activities is shown in Figure 3.

##  4\. Data Preprocessing

Before developing algorithms for the classification of dog signalling
behaviour, it is important to understand the nature of the signals produced by
the triaxial accelerometer. These consist of three separate data streams that
represent time series for acceleration on each axis x, y and z. Complementary
to the three axes, an additional time series, the named sum can been obtained
by computing the magnitude of the acceleration described into signal vector
magnitude (SVM), a time-series independent of sensor orientation and invariant
to any movement of the collar around the dog’s neck. SVM is also known as Amag
and shown below, is a useful metric for calculating movements which are not
axis-specific and is used extensively in pre-processing raw accelerometer data
(Riaboff et al., 2022).

$Amag=\sqrt{x^{2}+y^{2}+z^{2}}$

Some examples of raw accelerometer recordings are in Figure 6 indicating
indicates a variety of movements in all directions being captured on the
sensor in the top graph and an idle stage without any movements on the sensor
indicating the dog is still in the bottom graph. Raw accelerometer data was
recorded for all the participant dogs at 12Hz, giving 12 values of raw
accelerometer data for each second. In order to obtain the most important and
discriminative statistical features of the data we applied feature engineering
where the sensor data was used to derive several features. An important
advantage of feature engineering is that it allows the most relevant and
important features to be extracted from a collection of data. Our feature
extraction process is outlined below.

Figure 6. Sample recordings of dog activity. The top graph shows movement in
all 3 directions while the bottom graph shows a dog at rest with no movement.

###  4.1. Feature Extraction

Feature sets are reduced sets of original data that represent main
characteristics and behaviours and are thus abstractions of raw data. A
feature vector is a subset of larger input data that contains important clues
for activity recognition and is used as an input into clustering or
classification algorithms. In our work we used the Time Series Feature
Extraction Library (TSFEL) written in Python which automatically extracts over
500 features from a time series (Barandas et al., 2020). These include Fast
Fourier Transform (FFT), Wavelet Transform (WT), Autocorrelation, Mean,
Median, Sum of Absolute Differences, Entropy, Peak to Peak Distance, Slope,
Zero crossing rate, Interquartile range, ECDF percentile count and ECDF slope,
Spectral maximum peaks, Maximum power spectrum, etc. TSFEL can handle
multidimensional time series, as is the case here.

In feature selection it is advantageous to use smaller feature sets that are
fast to compute when there is a shortage of computational resources in the
applications. While working with larger sets of complex time-series features
can be slower to compute, it can result in higher performance in downstream
applications where accuracy is more important. Feature extraction is usually
applied to windows taken from a time series and testing is required to
determine how the signal window size affects subsequent classification
performance.

###  4.2. Relevant Features

The characteristics of a given classification task may require different
features to be extracted from data, and some variables in the original data
may be irrelevant or redundant. A feature selection technique is used to
identify the most informative features and to limit the computational demands
when applying the recognition system to new observations. The feature
selection stage selects a smaller subset of the original features to identify
those most useful features. The choice of features acquired from a data set
and the window length over which they are computed are key factors.

In the work addressed here, recognising dog signalling behaviour depends
directly on features extracted for motion analysis. The same signalling
behaviour movement can be performed by different dogs in different ways
depending on their size, speed, and age, resulting in considerable variability
in the features derived from the data collected by collar-worn sensors. It is
therefore important to identify features with high discriminative ability in
order to achieve effective classification. An effective feature set should
show little variation between repetitions of the same movements performed by
the same individual dog, but should differ between different dogs performing
the same task. To achieve this we extract only features related to motion and
energy. This will enhance the accuracy of predicting a movement across a
variety of sensor readings.

A dictionary of features is developed which categorises them into three types:
statistical, temporal, and spectral as indicated in Table 3. Some of these
have been extensively investigated in previous studies and have proven useful
in recognising activity and motion (Vázquez Diosdado et al., 2015). A window
size of 12 data values was pre-selected to correspond to approximately 1
second duration since the frequency in the AX3 is 12.5Hz, and a window overlap
of 50% was set for our feature extraction.

Table 3. Motion features for classification from time series Statistical features | Temporal features |  Spectral features
---|---|---
Kurtosis; | Zero crossing; |  Maximun frequency;
Skewness; | Pairwise correlation; |  Median frequency;
Mean; | Auto Correlation; |  Cepstral Coefficients;
Standard Deviation; |  |  Power Spectrum;
Interquartile Range; |  |  Power Bandwidth;
Root Mean Square; |  |  Fundamental
Median Absolute Deviation; |  |  Frequency Spectrum;

The final step in data preprocessing and feature extraction is to address the
difference in magnitude among the features obtained from TSFEL as several
machine learning algorithms emphasise more heavily those features with a
higher magnitude (Barandas et al., 2020). A normalising step is thus applied
to features before classification to prevent this unwanted effect. To account
for different scale factors and units, all described features are normalised
to zero mean and unit variance before proceeding to algorithm development and
analysis.

###  4.3. Algorithm Development

In order to verify if the selected features selected are informative, we use
different classification methods for various dog signalling behaviours.
Because the optimisation criteria are convex, Support Vector Machine
classifiers are desirable in this scenario. The support vector machine is a
linear classifier based on supervised learning, but it can also be used as a
nonlinear classifier by using different kernel functions (Cervantes et al.,
2020). In Support Vector Machine learning, the basic idea is to find the
separation hyperplane that can divide the training data set correctly and have
the largest set intervals. We also trained and tested the data set on other
classification algorithms including logistic regression, KNN, Random forest
and Naive Bayes. All our results were validated by 10-fold cross validation.

Our experiments were conducted using the Python programming language and the
interface was Jupyter Labs. A standard M1 processor Macintosh with 16GB of CPU
memory was used and processing time varied across different classification
models.

##  5\. Experimental Results

In prediction and classification of assistance dog signalling behaviour,
conventional machine learning algorithms based on accelerometer data can be
used. To distinguish between different dog behaviours and features extracted,
supervised learning is used in this work. Due to the availability of annotated
training data, supervised approaches are thus straightforward to implement.

In our experiments we ran different machine learning classification algorithms
on the features extracted from accelerometer data from the 11 assistance dogs.
The algorithms we used included Support Vector Machines, Logistic Regression,
kNN, Random Forest and Naive Bayes. A classifier for each of the 11
participant dogs was trained and tested on data from that dog only, using a
support vector machine classification model and the results for the detection
of dog spins (in either direction) are indicated in Table 4.

Table 4. Support vector machine classification results for individual classifiers for 11 participant dogs  Number | Dog Breed | Precision | Recall | F1 Score | Support | Accuracy
---|---|---|---|---|---|---
1 | Collie Mix | 0.67 | 0.71 | 0.69 | 17 | 0.78
2 | Springer/Terrier Mix | 0.67 | 1.00 | 0.8 | 4 | 0.98
3 | Spaniel | 0.46 | 0.5 | 0.55 | 11 | 0.52
4 | Poodle | 0.38 | 0.56 | 0.45 | 9 | 0.65
5 | Lurcher | 0.74 | 0.8 | 0.77 | 14 | 0.79
6 | Lurcher | 0.17 | 0.5 | 0.25 | 2 | 0.7
7 | Collie Mix | 0.64 | 0.5 | 0.56 | 14 | 0.65
8 | Hungarian Vizla | 1.00 | 0.77 | 0.87 | 13 | 0.77
9 | Border Collie | 0.67 | 1.00 | 0.8 | 2 | 0.94
10 | Australian Retriever | 1 | 0.93 | 0.97 | 15 | 0.99
11 | Golden Retriever | 0.75 | 0.43 | 0.55 | 7 | 0.76
| Average | 0.64 | 0.73 | 0.67 |  | 0.78

The results in Table 4 show a range of performances depending on the dog. Of
the 4 evaluation measures included, recall is the most important for our use
case. When an assistance dog detects a seizure onset from its owner and
exhibits signalling behaviour, we always want that behaviour to be detected.
Our use case can tolerate false positives where the algorithm incorrectly
interprets a behaviour from the sensor data as being a signal because we can
discard a false alert, but we do not want to miss a signalling behaviour,
either true or false, in case it is true. Recall scores range from 0.5 to 1.0
across different dogs with an average of0.73 and while some of the per dog
performance figures may be disappointing we are not interested in training a
machine learning classifier to recognise signalling behaviour from each dog
independently, we want to pool the annotations to train a single classifier
that operates across all dogs. Nevertheless, the performance of the classifier
for some dogs, notably dogs numbered 2 (Springer/Terrier), 9 (Border Collie)
and 10 (Australian Retriever) was good. All of these are at the larger end of
the medium sized breeds with faster movements.

The combined training data of all 11 participant dogs was used to train and
test a classifier for dog spinning behaviour in either direction using the set
of 5 classification model techniques mentioned earlier. Validation tests were
performed using 10-fold cross validation algorithms individually for each
participant and with the combined training data. The results based on the
combined training data from all 11 participants are indicated in Table-5 using
5 different machine learning algorithms. Our best overall results were
achieved using the random forest and the naive Bayes algorithms with naive
Bayes having the best recall score at 0.98. Performance results for the
support vector machine, logistic regression and kNN were not as good. It
should be noted that these results were also calculated using an 10-fold cross
validation to further strengthen their validity and demonstrate an acceptable
level of performance for this task when used in real world scenarios.

Table 5. Classification Results for classifiers using training data taken from all dogs  Number | Classifier | Precision | Recall | F1 Score | Support | Accuracy
---|---|---|---|---|---|---
1 | Support vector machine | 0.84 | 0.62 | 0.71 | 52 | 0.71
2 | Logistic Regression | 0.74 | 0.67 | 0.70 | 50 | 0.70
3 | K-NN | 0.33 | 0.08 | 0.12 | 52 | 0.84
4 | Random Forest | 0.95 | 0.77 | 0.85 | 52 | 0.96
5 | Naive Bayes | 0.95 | 0.98 | 0.97 | 52 | 0.98

##  6\. Conclusions

The combination of data from wearable sensors and machine learning is an
exciting new frontier and opportunity. The contribution of this work is the
design and evaluation of a behaviour recognition system based on data from a
wearable accelerometer that identifies the presence of a pre-trained
signalling behaviour. The set of features automatically extracted from the
sensor data ensures acceptable recall in terms of classification results. The
results presented here suggest that using a realtime equivalent of the AX3
accelerometer can identify dogs’ signalling behaviour and for an assistance
dog trained to alert when sensing the onset of an epileptic seizure, this can
then trigger an alert to the owner or carer.

One possibility suggested to further improve the quality of the classification
is to use additional wearable sensors besides accelerometry. (Arablouei et
al., 2022) has worked with accelerometry fused with GNSS or GPS data but this
or any other form of location information would not be useful in our use case
since an assistance dog’s signalling behaviour could occur anywhere, including
indoors.

Future research is needed to continue to develop and refine the model to
identify additional signalling behaviours that could be performed easily by
dogs and for which they could easily be trained. These could include
rollovers, beg-ups, and other distinct behaviours.

## References

  * (1)
  * Arablouei et al. (2022) Reza Arablouei, Ziwei Wang, Greg J Bishop-Hurley, and Jiajun Liu. 2022.  Multi-modal Sensor Data Fusion for In-situ Classification of Animal Behavior Using Accelerometry and GNSS Data.
  * Barandas et al. (2020) Marília Barandas, Duarte Folgado, Letícia Fernandes, Sara Santos, Mariana Abreu, Patrícia Bota, Hui Liu, Tanja Schultz, and Hugo Gamboa. 2020\.  TSFEL: Time series feature extraction library.  _SoftwareX_ 11 (2020), 100456.
  * Bhuvaneswari et al. (2020) M Bhuvaneswari, J Shafiq Mansoor, T.S Aravinth, M Arun Kumar, and R Sangeetha. 2020\.  An Automatic System for Detection and Suppression of Epileptic Seizure. In _2020 4th International Conference on Electronics, Communication and Aerospace Technology (ICECA)_. 1524–1529.  https://doi.org/10.1109/ICECA49313.2020.9297423
  * Brown and Goldstein (2011) Stephen W Brown and Laura H Goldstein. 2011.  Can seizure-alert dogs predict seizures?  _Epilepsy Research_ 97, 3 (2011), 236–242.
  * Catala et al. (2018) Amelie Catala, Hugo Cousillas, Martine Hausberger, and Marine Grandgeorge. 2018.  Dog alerting and/or responding to epileptic seizures: A scoping review.  _PloS ONE_ 13, 12 (2018), e0208280.
  * Cervantes et al. (2020) Jair Cervantes, Farid Garcia-Lamont, Lisbeth Rodríguez-Mazahua, and Asdrubal Lopez. 2020.  A comprehensive survey on support vector machine classification: Applications, challenges and trends.  _Neurocomputing_ 408 (2020), 189–215.
  * Chambers et al. (2021) Robert D Chambers, Nathanael C Yoder, Aletha B Carson, Christian Junge, David E Allen, Laura M Prescott, Sophie Bradley, Garrett Wymore, Kevin Lloyd, and Scott Lyle. 2021.  Deep learning classification of canine behavior using a single collar-mounted accelerometer: Real-world validation.  _Animals_ 11, 6 (2021), 1549.
  * den Uijl et al. (2017) Ingrid den Uijl, Constanza B Gómez Álvarez, David Bartram, Yoni Dror, Robert Holland, and Alasdair Cook. 2017.  External validation of a collar-mounted triaxial accelerometer for second-by-second monitoring of eight behavioural states in dogs.  _PloS ONE_ 12, 11 (2017), e0188481.
  * Hussain et al. (2022) Ali Hussain, Sikandar Ali, Abdullah, and Hee-Cheol Kim. 2022\.  Activity Detection for the Wellbeing of Dogs Using Wearable Sensors Based on Deep Learning.  _IEEE Access_ 10 (2022), 53153–53163.  https://doi.org/10.1109/ACCESS.2022.3174813
  * Kumpulainen et al. (2021) Pekka Kumpulainen, Anna Valldeoriola Cardó, Sanni Somppi, Heini Törnqvist, Heli Väätäjä, Päivi Majaranta, Yulia Gizatdinova, Christoph Hoog Antink, Veikko Surakka, Miiamaaria V Kujala, et al. 2021\.  Dog behaviour classification with movement sensors placed on the harness and the collar.  _Applied Animal Behaviour Science_ 241 (2021), 105393.
  * Litt and Echauz (2002) Brian Litt and Javier Echauz. 2002.  Prediction of epileptic seizures.  _The Lancet Neurology_ 1, 1 (2002), 22–30.
  * Maa et al. (2021) Edward Maa, Jennifer Arnold, Katherine Ninedorf, and Heather Olsen. 2021.  Canine detection of volatile organic compounds unique to human epileptic seizure.  _Epilepsy & Behavior_ 115 (2021), 107690.  https://doi.org/10.1016/j.yebeh.2020.107690
  * Majumder and Deen (2020) Sumit Majumder and M. Jamal Deen. 2020.  A Robust Orientation Filter for Wearable Sensing Applications.  _IEEE Sensors Journal_ 20, 23 (2020), 14228–14236.  https://doi.org/10.1109/JSEN.2020.3009388
  * Migueles et al. (2017) Jairo H Migueles, Cristina Cadenas-Sanchez, Ulf Ekelund, Christine Delisle Nyström, Jose Mora-Gonzalez, Marie Löf, Idoia Labayen, Jonatan R Ruiz, and Francisco B Ortega. 2017.  Accelerometer data collection and processing criteria to assess physical activity and other outcomes: a systematic review and practical considerations.  _Sports Medicine_ 47, 9 (2017), 1821–1845.
  * Naik et al. (2020) Hemal Naik, Renaud Bastien, Nassir Navab, and Iain D Couzin. 2020.  Animals in Virtual Environments.  _IEEE Transactions on Visualization and Computer Graphics_ 26, 5 (2020), 2073–2083.  https://doi.org/10.1109/TVCG.2020.2973063
  * Powell et al. (2021) Neil A Powell, Alastair Ruffell, and Gareth Arnott. 2021\.  The untrained response of pet dogs to human epileptic seizures.  _Animals_ 11, 8 (2021), 2267.
  * Riaboff et al. (2022) L Riaboff, L Shalloo, AF Smeaton, S Couvreur, A Madouasse, and MT Keane. 2022\.  Predicting livestock behaviour using accelerometers: A systematic review of processing techniques for ruminant behaviour prediction from raw accelerometer data.  _Computers and Electronics in Agriculture_ 192 (2022), 106610.
  * Rohan et al. (2020) Tanbin Islam Rohan, Md. Salah Uddin Yusuf, Monira Islam, and Shidhartho Roy. 2020.  Efficient Approach to Detect Epileptic Seizure using Machine Learning Models for Modern Healthcare System. In _2020 IEEE Region 10 Symposium (TENSYMP)_. IEEE, 1783–1786.  https://doi.org/10.1109/TENSYMP50017.2020.9230731
  * Smeaton et al. (2022) Alan F. Smeaton, Hitesh Raju, and Ankit Sharma. 2022\.  Accelerometer and annotation data for trained assistance dogs performing signalling behaviour on instruction.  https://doi.org/10.6084/m9.figshare.20459367.v1
  * Vázquez Diosdado et al. (2015) Jorge A Vázquez Diosdado, Zoe E Barker, Holly R Hodges, Jonathan R Amory, Darren P Croft, Nick J Bell, and Edward A Codling. 2015\.  Classification of behaviour in housed dairy cows using an accelerometer-based activity monitoring system.  _Animal Biotelemetry_ 3, 1 (2015), 1–14.
  * Yam et al. (2011) PS Yam, V Penpraze, David Young, MS Todd, AD Cloney, KA Houston-Callaghan, and John J Reilly. 2011.  Validity, practical utility and reliability of Actigraph accelerometry for the measurement of habitual physical activity in dogs.  _Journal of Small Animal Practice_ 52, 2 (2011), 86–91.
  * Yu et al. (2016) Pen-Ning Yu, Shokofeh A. Naiini, Christi N. Heck, Charles Y. Liu, Dong Song, and Theodore W. Berger. 2016\.  A sparse Laguerre-Volterra autoregressive model for seizure prediction in temporal lobe epilepsy. In _2016 38th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC)_. 1664–1667.  https://doi.org/10.1109/EMBC.2016.7591034
  * Zamansky et al. (2019) Anna Zamansky, Dirk van der Linden, Irit Hadar, and Stephane Bleuer-Elsner. 2019.  Log My Dog: Perceived Impact of Dog Activity Tracking.  _Computer_ 52, 9 (2019), 35–43.  https://doi.org/10.1109/MC.2018.2889637


# The WQN algorithm for EEG artifact removal in the absence of scale
invariance 111Funding: D.H.’s research is supported by grants ANR NEUC-0001,
PSL and CNRS pre-maturation, and by the European Research Council (ERC) under
the European Union’s Horizon 2020 research and innovation program (grant
agreement Nº 882673).

Matteo Dora<EMAIL_ADDRESS>Stéphane Jaffard<EMAIL_ADDRESS>David
Holcman<EMAIL_ADDRESS>
###### Abstract

Electroencephalogram (EEG) signals reflect brain activity across different
brain states, characterized by distinct frequency distributions. Through
multifractal analysis tools, we investigate the scaling behaviour of different
classes of EEG signals and artifacts. We show that brain states associated to
sleep and general anaesthesia are not in general characterized by scale
invariance. The lack of scale invariance motivates the development of artifact
removal algorithms capable of operating independently at each scale. We
examine here the properties of the wavelet quantile normalization algorithm, a
recently introduced adaptive method for real-time correction of transient
artifacts in EEG signals. We establish general results regarding the
regularization properties of the WQN algorithm, showing how it can eliminate
singularities introduced by artefacts, and we compare it to traditional
thresholding algorithms. Furthermore, we show that the algorithm performance
is independent of the wavelet basis. We finally examine its continuity and
boundedness properties and illustrate its distinctive non-local action on the
wavelet coefficients through pathological examples.

[1]organization=IBENS UMR8197, École Normale Supérieure – PSL, postcode=75005,
city=Paris, country=France [2]organization=Univ Paris Est Creteil, Univ
Gustave Eiffel, CNRS, LAMA UMR8050, F-94010 Creteil, France, postcode=F-94010,
city=Créteil, country=France

##  1 Introduction

For nearly a century, brain activity has been measured using
electroencephalography (EEG), a technique that uses electrodes placed on the
scalp of a patient to record the electrical activity of the brain [1]. This
physiological signal reflects the collective activity of neuronal populations
[2]. By analyzing the statistical and spectral properties [3] of the EEG, it
is possible to identify transient brain oscillations within a frequency range
between zero and hundreds of Hz that reflects key cognitive events, such as
specific responses to sensory stimulations, learning and memory, sleep stages,
meditation, coma, and more.

Because of its non-invasive nature, EEG recordings have been widely adopted in
the clinical setting for screening or monitoring tasks. An example is general
anaesthesia (GA), a procedure consisting in placing the brain into an
artificial but reversible coma state, which can now be routinely monitored in
real-time during surgery by recording the EEG signal from few electrodes. This
setting provides a continuous feedback about the depth of anaesthesia,
allowing for accurate control of the anaesthetic dose required to keep the
patient in a safe unconscious state.

However, the amplitude of the EEG signal typically varies in the microvolt
range, making the EEG highly susceptible to contamination by artifacts
originating from several types of sources. The issue is particularly
pronounced in the clinical setting, where there is a limited control over the
environment and the potential for artifact contamination. Artifacts include
noise from electrical equipments, muscle contractions, eye movements, or small
displacements of the electrodes, that alter the EEG signal.

Eliminating artifacts from the EEG signal is thus a crucial concern. Wavelet-
based methods [4, 5, 6] have been used extensively in this regard by
thresholding coefficients to remove artifacts from the EEG signal, taking
advantage of the different properties of the artifact and the physiological
EEG [7, 8, 9] distributions. In this direction, we recently introduced an
empirical method, the WQN algorithm [10, 11], designed to remove artifacts
from single-channel EEGs for real-time applications in clinical monitoring.
This adaptive approach allows attenuating transient artifacts in the EEG by
normalizing the wavelet coefficient distribution during the artifact so that
it matches the one obtained from a preceding uncontaminated interval.

While our previous studies have demonstrated the high effectiveness of the WQN
algorithm in removing transient artifacts from EEG [10, 11], the underlying
reasons for its efficiency remain unclear. In the current article, we provide
a comprehensive analysis of the WQN algorithm in two directions: First, we
study the statistical properties of the EEG signal, and second, we derive
several properties of the WQN algorithm.

The manuscript is organized as follows. In Section 2, we use multifractal
analysis techniques [12, 13, 14] in order to characterize the scaling
behaviour of different classes of EEG signals and artifacts, showing that
scale-invariance is not always observed. This result justifies the need to
modify wavelet coefficients separately for each scale, as proposed by the WQN
algorithm. In Section 3, we describe the WQN algorithm and study its
properties. In particular, we show that the WQN algorithm cannot introduce
unwanted singularities in the signal. We then show that the WQN algorithm is
robust with respect to changes of the wavelet basis. We also put in evidence
its non-local action on the wavelet coefficients and explore some of its
consequences; finally we demonstrate how it performs with respect to
pathological cases where signals are perturbed by different types of
singularities and random noise.

##  2 Scaling properties of the EEG signal

In the spectral analysis of the EEG signal, the scale invariance property is
characterized by a power-law decay present in the power spectrum

$S(f)\sim 1/f^{a},$ (1)

with exponent $a$. A scale invariant behavior can provide valuable insights
into the underlying mechanisms of neuronal networks [15]. In this section we
investigate in which contexts the EEG signal shows such scale-invariant
behaviour. We recall that power spectrum of the activity of local neuronal
ensembles can be fitted by a power law [16], however the EEG signal often
contains additional oscillatory activity revealed by the presence of specific
brain waves, such as the $\alpha$ wave in the range 8–12 Hz, which cannot be
filtered without significantly altering the properties of the signal.

###  2.1 Wavelet multifractal analysis framework

To estimate the scaling properties of the EEG, we resort to the setting
supplied by multifractal analysis [13, 17, 18], which extends the notion of
scale invariance to signals that cannot be characterized by a single scaling
exponent; this approach allows to overcome the limitation to second-order
statistics of the power spectrum by replacing the Fourier transform with a
multiresolution tool such as the wavelet transform. Practically, in the
wavelet-based multifractal framework, the scale invariance of a signal $x(t)$
is measured through a wavelet scaling function $\eta_{x}(q)$, which
encapsulates the scaling exponents associated with the moments of order $q$
($q>0$) of the wavelet coefficients. This approach can be seen as a
generalization of traditional spectral scale-invariance analysis (as in eq.
1), which is recovered for $q=2$ [18].

To define the scaling function $\eta_{x}$, we start with the wavelet
decomposition for functions defined on $\mathbb{R}$ [19, 20, 21, 13]. We use a
mother wavelet $\psi$, such that an orthonormal basis of $L^{2}(\mathbb{R})$
is obtained by dilations and translations of $\psi$:

$\left\\{\psi_{j,k}\left(t\right)=2^{-j/2}\,\psi\left(2^{-j}\,t-k\right),\,j,k\in\mathbb{Z}\right\\}.$
(2)

For a given signal $x(t)$, its decomposition on the wavelet basis is given by

$x(t)=\sum_{j,k\in\mathbb{Z}}c_{j,k}\,\psi_{j,k}(t),$ (3)

where $c_{j,k}=\left<x,\psi_{j,k}\right>$ are the wavelet coefficients222The
scalar product is defined by $\left<f,g\right>=\int_{\mathbb{R}}f(t)\,g(t)dt$.
of $x$ at scale $j$. We define the wavelet structure functions of $x$ based on
its normalized wavelet coefficients, as follows:

$\forall q>0,\qquad S_{x}(j,q)=2^{j}\sum_{k}\left|2^{-j/2}c_{j,k}\right|^{q},$
(4)

and finally the wavelet scaling function $\eta_{x}(q)$ is implicitly defined
by

$\forall q>0,\qquad S_{x}(j,q)\sim 2^{\eta_{x}(q)j},$ (5)

in the limit of small scales, i.e. when $j\rightarrow-\infty$. The
mathematical interpretation is

$\eta_{x}(q)=\liminf_{j\rightarrow-\infty}\frac{\log_{2}(S_{x}(j,q))}{j}.$

This quantity provides information about the global regularity of $x$. The
particular value $q=2$ allows to recover the information supplied by the
traditional spectral analysis [18], i.e. the Hurst exponent
$\alpha=\eta_{x}(2)+1$ for eq. 1. Recall that the Sobolev spaces
$H^{s,q}(\mathbb{R})$ are composed of functions whose fractional derivatives
belong to $L^{q}(\mathbb{R})$; $\eta_{x}$ has the following function space
interpretation:

$\eta_{x}(q)=q\cdot\sup\\{s:f\in H^{s,q}(\mathbb{R})\\}.$

In practice, we shall compute the wavelet scaling function $\eta_{x}$ using a
log-log plot regression, i.e. fitting $\log_{2}S_{x}(j,q)$ as a linear
function of $j$, in the range of scales where the EEG signal is most
informative (typically in the range 0.1–100 Hz). Note that the multifractal
framework is pertinent only if the log-log plot is approximately linear in a
sufficiently large range of scales, i.e. if the data exhibit an average
scaling invariance. We will examine below whether this holds or not for
different types of EEG signals.

Another relevant parameter which has proved useful in the setting of scale
invariance is the uniform Hölder exponent $H^{\text{min}}_{x}$ which can be
defined as follows. The largest normalized wavelet exponent at scale $j$ is
computed as

$D_{x}(j)=\sup_{k}2^{-j/2}|c_{j,k}|,$

then the scaling exponent $H^{\text{min}}_{x}$ is implicitly defined by

$D_{x}(j)\sim 2^{H^{\text{min}}_{x}j}\quad\mbox{ when }j\rightarrow-\infty,$

i.e.

$H^{\text{min}}_{x}=\liminf_{j\rightarrow-\infty}\frac{\log_{2}(D_{x}(j))}{j}.$

Its function space interpretation is obtained through the Hölder spaces
$C^{s}(\mathbb{R})$:

$H^{\text{min}}_{x}=\sup\\{s:x\in C^{s}(\mathbb{R})\\}.$

Similarly to the scaling function, the coefficient $H^{\text{min}}_{x}$ can be
estimated by log-log regression of $D_{x}(j)$ versus $2^{j}$. As above, this
approach is only applicable when the log-log plot is well described by a
straight line on a relevant range of scales.

###  2.2 Multifractal analysis of the EEG signal

We now investigate the scaling properties of EEG signals under this
multifractal framework. To test whether EEGs present a general scale-invariant
structure, we analyze signals acquired in different conditions:

  * 1.

EEG recorded during active task execution [22];

  * 2.

EEG for subjects at rest, with eyes closed [23, 24];

  * 3.

EEG recorded during sleep [25, 24];

In general, the high frequency content of EEGs is higher during mental
activity and decreases with rest and sleep, as slow and more regular patterns
progressively emerge. We visualize this tendency in figure 1, where we present
a comparison of three EEG signals: during task (fig. 1-A1), rest (fig. 1-A2),
and sleep (fig. 1-A3). The emergence of the slow patterns associated with
sleep is appreciable already in the time domain (compare A1 to A3). The power
spectrum of the EEG during task execution can be well fitted by a power law
(fig. 1.B1), while rest and sleep EEG show significant deviations (fig.
1.B2–B3). In particular, EEG during rest is characterized by a strong alpha
wave ($8\text{\,}\mathrm{Hz}12\text{\,}\mathrm{Hz}$), manifested as a peak in
the power spectrum in the alpha-band frequencies (fig. 1.B2). The EEG during
sleep presents a predominance of lower frequencies in the delta
($1\text{\,}\mathrm{Hz}4\text{\,}\mathrm{Hz}$) and theta
($4\text{\,}\mathrm{Hz}8\text{\,}\mathrm{Hz}$) bands (fig. 1.B3), which
characterize non-rapid eye movement sleep. Thus, both rest and sleep EEG
examples present a clear deviation from a typical $1/f^{a}$ scale-invariant
behaviour.

We thus adopted the multifractal framework described above to test whether it
is possible to characterize these types of EEG under a more general definition
of scale invariance. For each EEG sample, we computed the wavelet structure
functions $S_{x}(j,q)$ and estimated the uniform Hölder exponent
$H^{\text{min}}_{x}$ using Daubechies wavelets with 3 vanishing moments [19]
(the result confirms that this choice yields sufficiently smooth wavelets to
analyze such data). In fig. 1.C we show the structure functions $S_{x}(j,q)$
for $q=1,2$. For the EEG during task, $\log_{2}S_{x}(j,q)$ can be well fitted
by a linear function of $j$ (fig. 1.C1, dashed lines), making it possible to
define a multifractal scaling via the scaling function $\eta_{x}(q)$. Note
that we restricted the fit to the frequency range
$0.1\text{\,}\mathrm{H}\mathrm{z}50\text{\,}\mathrm{H}\mathrm{z}$ which
contains the most relevant physiological information and is not affected by
filtering and acquisition limitations (indicated by the shaded area in fig.
1.C). Both EEG during rest and sleep show significant deviations from linear
behaviour of the structure functions in the log-log plots (fig. 1.C2–C3),
which prevents possible estimation of the respective scaling functions.
Similarly, while we can reasonably define a $H^{\text{min}}_{x}$ exponent for
the EEG during active task execution (fig. 1.D1), the result is unreliable for
the other two cases (fig. 1.D2–D3).

To conclude, we report here that scale invariance of the EEG signal depends on
the context of acquisition and is not a general feature. In particular, while
scale invariance can be defined for EEG during mental activity, the same
cannot be said of those brain states which are characterized by more regular
rhythms and patterns such as sleep or rest. This result is not surprising as
the EEG reflects processes occurring at various timescale (from milliseconds
[26] to minutes), generated by processes of different natures: ionic channels,
neuronal spiking and bursting, transient changes of the membrane potential, or
oscillations that reflect communication between different brain regions [27].
Depending on the context, the EEG signal can thus become dominated by patterns
at specific timescales which break a possible underlying scale invariance.
This observation justifies, in the context of EEG artifact correction to use
wavelet algorithms such as WQN, which act on each scale independently.

Figure 1:  Multifractal analysis of the EEG signal. A. Samples of EEG signals
recorded during task execution (A1), rest (A2), and sleep (A3). B. Power
spectral density plots for the three cases (log-log). C. Wavelet structure
functions for $q=1,2$, as described in eq. 4. The shaded area indicates the
range of scales corresponding to a frequency range
$0.1\text{\,}\mathrm{Hz}50\text{\,}\mathrm{Hz}$. The dashed lines indicate the
linear fit that can be used to estimate the scaling function (see eq. 5). D.
Estimations of uniform Hölder exponent $H^{\text{min}}$. EEG during task
execution shows a good fit (D1), while rest and sleep give unreliable results
(D2–3).

###  2.3 Wavelet coefficient distributions for EEG artifacts

We now characterize the EEG and artifactual signals by examining the
distribution of their wavelet coefficients. We considered EEGs containing two
common types of artifacts, ocular (EOG) and muscular (EMG), and compare them
to uncontaminated EEGs. Ocular artifacts are caused by eye movements, which
can alter the recorded electric potential since the eye acts as a dipole (with
a difference of potential between the cornea and the fundus) while muscular
artifacts are derived from the interfering electrical activity produced by
muscle contraction. We decomposed two-second signals from the Denoise-Net
dataset [28] using a Daubechies wavelet with two vanishing moments to obtain
the coefficient distributions for each scale, shown in fig. 2.A.

The distributions of wavelet coefficients for the uncontaminated EEG can be
approximated by a normal distribution (fig. 2.A, first row, dashed line) with
minor deviations. In the case of ocular artifacts (fig. 2.A, second row), the
non-zero coefficients are concentrated at large scales (low frequency) with
distinctive asymmetric distribution. Wavelet coefficients for muscular
artifacts (fig. 2.A, third row) show a deviation from the normal distribution,
but are spread across multiple scales. We investigate a possible scale
invariant behaviour by plotting the variance of the wavelet coefficients
versus the scale $j$ (fig. 2.B), corresponding to the wavelet scaling function
$S(j,q=2)$.

We notice that the variances of the uncontaminated EEGs can be approximated by
a linear relation in log-log scale, thus showing some form of scale
invariance, as shown in fig. 1. The variance of ocular artifacts decays
rapidly with smaller scales, while muscular artifacts are characterized by an
almost constant variance across all scales.

To conclude, we have shown that EEG artifacts present significant deviations
from the uncontaminated EEG signal both in distribution of the wavelet
coefficients and their persistence across scales. Moreover, significant
differences in the scaling behaviour can be observed between distinct families
of artifacts, such as ocular and muscular.

Figure 2:  Characterization of ocular and muscular artifacts compared to EEG
in the wavelet domain. A. Distribution of the wavelet coefficients at
different scales for EEG, EOG (ocular artifacts), and EMG (muscular
artifacts). The EEG distributions can be fitted by Gaussian distribution with
minor deviations (dashed line, standard deviation $\sigma$). B. Plot of the
variance versus scale. Variances for the pure EEG can be well approximated by
a linear function (in log-log scale), compatibly with scale invariance shown
in fig. 1, while artifacts show deviations.

##  3 Properties of the WQN algorithm

We briefly recall the wavelet quantile normalization (WQN) algorithm [10]. WQN
is an adaptive method that allows to attenuate transient artifacts based on
statistics estimated from clean regions of the signal. The algorithm was
especially designed to be applied on EEG signals in the context of brain
monitoring, where signal corruption is a consequence of artifacts generated by
motion of electrodes, muscular activity or eye motion that alter the
physiological signal generated by neuronal activity [29, 1].

First, the signal is decomposed on a wavelet basis according to eq. 3. We
assume here that the artifacted intervals are well identified and isolated.
They usually consist of a small portion of the total EEG signal. Once an
artifact segment has been localized, it is associated to a clean reference
segment where we expect the underlying signal to be similar but uncontaminated
by artifacts. Note that these uncontaminated statistics depend on the patient
and timing so that they cannot be acquired in advance.

In the original implementation, such reference intervals are defined by
considering a temporally adjacent uncontaminated portion of the signal of
roughly the same length as the contaminated portion. For each decomposition
scale $j$, we define the wavelet coefficients associated with the artifacted
and reference intervals by $c^{\text{(art)}}_{j}$ and $c^{\text{(ref)}}_{j}$
respectively. In practice, since we work with finite intervals, the wavelet
decomposition is carried on up to scale $M$ (i.e. $j=1,\dots,M$) guaranteeing
the presence of a sufficient number of coefficients for the largest scale
(e.g. $\text{Card}\left\\{c^{\text{(art)}}_{M,k}\right\\}>30$).
In the second step, the wavelet coefficients $c^{\text{(art)}}_{j}$ of the
artifacted sample are modified in order to fit the statistics of the reference
coefficients $c^{\text{(ref)}}_{j}$ at the corresponding scale $j$. This
normalization is performed, scale by scale, by computing the empirical
cumulative density functions (CDF) $F^{\text{(ref)}}_{j}$,
$F^{\text{(art)}}_{j}$ of the coefficients amplitude for artifacted and
reference signal respectively, defined as

$F_{j}(x)=\frac{1}{N_{j}}\sum_{k=1}^{N_{j}}{1}_{\left|c_{j,k}\right|<\,x},$
(6)

where $N_{j}$ indicates the number of coefficients at scale $j$ and
${1}_{\left|c\right|<\,x}$ is the indicator function which takes value 1 if
$\left|c\right|<\,x$ and 0 otherwise. The wavelet coefficients
$c^{\text{(art)}}_{j}$ are then modified so that the distribution of their
amplitude matches that of the reference segment, via the mapping $T_{j}$
defined as

$T_{j}(x)=F^{\text{(ref)}\,-1}_{j}\left(F^{\text{(art)}}_{j}(x)\right),$ (7)

where $F^{\text{(ref)}\,-1}_{j}$ indicates the generalized inverse of
$F^{\text{(ref)}}_{j}$ (in the sense of completed graphs for discontinuous
increasing functions) see [11] and Fig. 3. Finally, the normalization function

$\displaystyle\lambda_{j}(c)=sign(c)\cdot\min\left\\{\left|c\right|,\,\left|T_{j}\left(c\right)\right|\right\\},$
(8)

maps a coefficient $c$ from $c^{\text{(art)}}_{j}$ to its possibly attenuated
value. The corrected coefficients are thus defined by

$c^{\text{(corr)}}_{j,k}=\lambda_{j}\left(c^{\text{(art)}}_{j,k}\right).$ (9)

Equation 8 ensures that the norm of wavelet coefficients is never increased, a
key requirement to guarantee the regularity of the algorithm, as we will show
in section 3.2. The corrected version of the signal is obtained by replacing
the artifacted coefficients $c^{\text{(art)}}_{j,k}$ by the corrected
coefficients $c^{\text{(corr)}}_{j,k}$ and then inverting the discrete wavelet
transform.

Figure 3:  Correction of EEG signals perturbed by ocular and muscular
artifacts. A1. EEG signal corrupted by ocular artifact (EOG) and its
reconstruction by WQN. A2. Left: cumulative density function and its mapping
from artifact (orange) to reference (black). Right: Projections on the wavelet
basis showing attenuation of the artifactual components on different scales.
B1–B2. Similar presentation of WQN correction for an EEG signal contaminated
by muscular artifact (EMG).

###  3.1 Illustration of the WQN algorithm

We present two examples to illustrate how the WQN algorithm works and how it
adapts to artifacts affecting different scales. To this aim, we have
contaminated an EEG signal by adding an ocular artifact (EOG) and a muscular
artifact (EMG), as illustrated in fig. 3 (A1–B1, first row). In fig. 3.A2-B2
(left) we represent how the wavelet coefficients are transported from the
artifacted cumulative density function to the reference one. In the case of
ocular artifact (fig. 3.A2) most of the transport occurs in the first scales
(level 5), corresponding to low frequency components, while shorter scales
(levels 1–5) which are not affected by the EOG artifact are left almost
unmodified. The reconstruction of the signal, which can be seen as the sum of
the corrected wavelet projections, is shown in fig. 3.A1 (blue). The energy of
the EOG artifact is mostly concentrated in the low frequencies, corresponding
to the approximation coefficients $c_{5}^{app}$ (the coarsest scale of the
wavelet decomposition). Appropriately, most of the correction takes place at
this scale (fig. 3.A2, first row).

In a second example (fig. 3.B) we present an EEG signal corrupted by a
muscular artifact (EMG), which is characterized by a wide frequency signature
(as shown in fig. 2), with high frequency components significantly more
powerful than the EEG signal. In fig. 3.B2, we show how the WQN algorithm
adapts to this different artifactual signature by attenuating the artifactual
component in both shorter scales (levels 1–3) and larger scales (level 5). In
conclusion, although the EMG and EOG artifacts are characterized by different
statistics, the adaptive approach of the WQN algorithm makes it effective at
reconstructing the original signal by performing appropriate corrections on a
scale-by-scale basis.

###  3.2 Regularity properties of the WQN algorithm

In this section, we use the wavelet decomposition to derive a regularization
property of the WQN algorithm. We recall that wavelets are unconditional bases
of most classical function spaces [12], such as Sobolev spaces $H^{s,p}$ for
$1<p<\infty$ or closely related Besov spaces $B^{s,q}_{p}$ for $0<p,q<\infty$
see e.g. [20, 13]. This implies that these spaces have a wavelet
characterization which bears on the moduli of the wavelet coefficients, and
which is an increasing function of each of these moduli. For example, a
function $f$ belongs to $B^{s,q}_{p}$ if its wavelet coefficients satisfy the
condition

$\left(\sum_{n}(2^{(1/p-1/2-s)m}|c_{m,n}|)^{p}\right)^{1/p}\in l^{q},$ (10)

and (10) yields a norm (or a quasi-norm, when $p$ or $q$ is less than 1) which
is equivalent to the Besov norm.
By construction, the WQN algorithm is wavelet decreasing (i.e. it does not
increase the size of the wavelet coefficients) and thus has the following
regularity property: For $0<p,q\leq\infty$, and for any $s\in\mathbb{R}$, it
maps functions of a Besov or Sobolev space into the same space. Since the
mapping is not linear, this property does not imply that it is continuous on
the corresponding functional space; this question is relevant as the
continuity of a denoising algorithm is a prerequisite to guarantee its
numerical robustness. In order to investigate this problem, it is useful to
compare the present algorithm with the classical wavelet thresholding and
wavelet shrinkage algorithms, which are conceptually simpler, and where the
same problem arises.

####  3.2.1 Comparison with wavelet thresholding and wavelet shrinkage

Wavelet thresholding and wavelet shrinkage were introduced to perform
denoising of signals or images without smoothing the signal to be recovered
(in contradistinction with convolution-based techniques). In this context, the
eliminated “noise” is assumed informally to be characterized with statistics
of wavelet coefficients which, at a given scale, are assumed to be stationary,
of small amplitude and with short range correlations. The algorithm is
efficient if the statistical properties of the signal to be recovered strongly
differ from the noise i.e. if its wavelet coefficients form a sparse sequence
(most of them almost vanish), and the other coefficients are large. Note that
this situation is opposite to the one we consider in the present article,
where the artifacts to be eliminated have a sparse signature while the signal
to be recovered presents the statistical properties of such a noise. However,
wavelet thresholding and wavelet shrinkage have also been used in such
contexts: once the splitting has been performed, one keeps the “noisy” part
instead of the sparse one [8, 30, 9]. It is therefore legitimate to compare
their performance with the WQN algorithm.
We briefly recall the wavelet thresholding. Given a threshold level $t>0$,
wavelet thresholding in a given wavelet basis is defined as follows: once an
appropriate normalization of the wavelet coefficients has been chosen, the
wavelet thresholding mapping is the operator $T$ which maps the wavelet
coefficient $c_{m,n}$ to

$d_{m,n}=f_{t}(c_{m,n}).$

where

$f_{t}(x)=x1_{[-t,t]}(x).$ (11)

The mapping $T$ is wavelet decreasing, and therefore maps functions in a Besov
or Sobolev space to the same function space. Nonetheless, the function $f_{t}$
is discontinuous and consequently, the operator $T$ is not continuous on any
of these spaces. Indeed, to show this property, we consider two functions $f$
and $g$ with respectively wavelet coefficients $c^{1}_{m,n}$ and $c^{2}_{m,n}$
which coincide, except for one wavelet coefficient, such that these
coefficients for $f$ and for $g$ respectively are $t-\varepsilon$, and
$t+\varepsilon$. Taking $\varepsilon$ arbitrarily small, the Besov norm of
$f-g$ can be made arbitrarily small, but the Besov norm of $T(f)-T(g)$ is (up
to the normalization factor of the corresponding wavelet coefficient)
$t+\varepsilon$, and therefore does not tend to zero, when
$\varepsilon\rightarrow 0$. This lack of continuity implies numerical
instabilities of the algorithm which are well documented, see e.g.[31, 32,
33].
We now discuss the wavelet shrinkage algorithm, which is based on the function

$g_{t}(x)=sgn(x)\cdot(|x|-t)^{+}.$ (12)

Once an appropriate normalization of the wavelet coefficients has been chosen,
the wavelet shrinkage operator $U$ is defined as mapping the wavelet
coefficient $c_{m,n}$ to

$e_{m,n}=g_{t}(c_{m,n}).$

The mapping $U$ is wavelet decreasing, and therefore maps functions in a Besov
or Sobolev space to the same function space. But, additionally, in
contradistinction with the previous case, $g_{t}$ is continuous, and therefore
two functions $f$ and $g$ whose wavelet coefficients are close are now mapped
to functions $U(f)$ and $U(g)$, which also have close wavelet coefficients.
More precisely, since $g_{t}$ is Lipschitz, it follows easily that, for a
given Besov or Sobolev space $E$, $\parallel U(f)-U(g)\parallel_{E}\leq
C\parallel f-g\parallel_{E}$, i.e. the mapping $U$ is Lipschitz in the
corresponding Besov space. Refined stability properties of wavelet shrinkage
can be attributed to $g_{t}$ continuity, i.e. (at least implicitly) to the
continuity of the mapping $U$, see [31, 32, 33].
One drawback of both of these algorithms is that they are local in the wavelet
domain, i.e. each wavelet coefficient is modified independently of the other
ones, and therefore, they do not preserve the statistics of wavelet
coefficients [11]. This phenomenon had no negative impact when wavelet
thresholding and wavelet shrinkage were used for their initial purpose, i.e.
to restore the “sparse” part of the signal, but it becomes a major drawback
when it is used in the opposite direction of recovering the “noisy” component,
in which case restoring the right statistics of the signal can be a major
issue. For instance, where the eliminated artefact was localized, wavelet
thresholding or shrinkage set the wavelet coefficients to zero, thus leading
to inhomogeneities in the restored signal. One of the purposes of the WQN
algorithm is to circumvent this drawback by restoring everywhere the correct
anticipated statistics of wavelet coefficients [10]. As a consequence, it is
not local in the wavelet domain. The value attributed to a coefficient depends
on the entire statistic of coefficients at a given scale, and therefore the
analysis of the regularity properties of the algorithm is more involved than
for wavelet thresholding and wavelet shrinkage; nonetheless, a preliminary
investigation of its main features will be performed in subsection 3.5.

###  3.3 Robustness with respect to wavelet basis

To examine the performance of the WQN algorithm, we use five classical wavelet
bases, currently used in signal processing (sym5, db5, coif3, bio3.5 and dmey)
[21]. After we added an electrode moving artifact on an EEG signal, we use the
WQN to remove the artifact (Fig. 4), resulting in the green curves in the
different sub-figures. To quantify this performance, we computed the Average
Root Mean Squared Error on the ensemble of data consisting in two types of
artifacts (EOG ad EMG). We can conclude that the RMSE is quite independent of
the wavelet bases with a mean around 0.03 for both types of artifacts (fig.
4.A).
Finally, we tested the effect of increasing the number of vanishing moments of
the wavelets. Again, we found that there was no noticeable consequence on the
corrected artifact with the three following wavelet bases: Daubechies, Symlets
and Coiflets. To conclude, the different wavelets bases have little influence
on the corrected artifacts.

###  3.4 Boundedness properties of the WQN algorithm in functional spaces

As already mentioned above, since the WQN algorithm does not increase the size
of the wavelet coefficients, it maps functions of a Besov or Sobolev space
into the same space. One can also consider other function spaces which are
based on histograms of wavelet coefficients at each scale and thus encapsulate
more functional information than Besov spaces do [12, 34]. The maximal
information which is invariant under the change of (smooth) wavelet basis, is
encapsulated through the wavelet profile of $f$, which is defined as follows.
For a function $f$, we define

$F_{m}(\alpha)=\text{Card}\left\\{n:\;\;|c_{n,m}|\geq
2^{(\alpha+1/2)m}\right\\};$

and the wavelet profile $\nu_{f}(\alpha)$ is

$\nu_{f}(\alpha)=\lim_{\varepsilon\rightarrow
0}\;\left[\limsup_{j\rightarrow\infty}\left(\frac{\log(F_{m}(\alpha+\varepsilon))}{\log(2^{j})}\right)\right].$

This definition formalizes the following heuristic: There are about
$2^{-\nu_{f}(\alpha)m}$ wavelet coefficients larger than $2^{(\alpha+1/2)m}$.
The corresponding function spaces are defined similarly: Let $\nu(\alpha)$ be
a nondecreasing function which takes values in $\\{-\infty\\}\cup[0,1]$. A
function $f$ belongs to the space $S^{\nu}$ if its wavelet coefficients
satisfy. $\forall\alpha\in\mathbb{R},$ $\forall\varepsilon>0,$

$\forall C>0,\;\exists M\;\forall m\leq M\;\;\;\;F_{m}(\alpha)\leq
2^{-(\nu(\alpha)+\varepsilon)m}.$

Since the several wavelet algorithms that we considered are wavelet
decreasing, it follows that these algorithms map functions which belong to a
$S^{\nu}$ space to the same $S^{\nu}$ space.

Figure 4:  Correction of electrode movement artifacts using the WQN algorithm.
We show the time-plot of the corrections (green) of an added artifact (orange)
using different wavelet bases: sym5, db5, coif3, bio3.5 and dmey [35]. Figure
5:  WQN performance under different wavelet bases. A. Average Root Mean
Squared Error (RMSE) of WQN for bootstrapped EEG signals contaminated by EOG
and EMG, for the following wavelet bases: Symlet with 5 vanishing moments
(sym5), Daubechies with 5 vanishing moments (db5), Coiflet with 6 vanishing
moments (coif3), biorthogonal spline with 3 and 5 vanishing moments in
synthesis and analysis wavelet respectively (bior3.5), discrete Meyer wavelet
(dmey). B. WQN performance for different number of vanishing moments,
calculated on the EOG contaminated dataset.

###  3.5 Continuity properties

After having considered the issue of the boundedness of the wavelet transport
algorithm on several classes of function spaces, we now turn to the problem of
its continuity. We start by a simple remark which will allow to position the
problem correctly. Assume that the unaltered data on which the histograms of
wavelet coefficients are recorded has a given finite length $L$. Then one
computes the number $N_{m}\sim L2^{-m}$ of wavelet coefficients at scale
$2^{m}$, which constitute the reference signal on which the altered data will
be mapped. The transport algorithm maps the wavelet coefficients of the
altered signal on this finite set of cardinality $N_{m}$. Now assume that, for
the altered signal, two wavelet coefficients of consecutive size are extremely
close; if their sizes are exchanged, the coefficients $c_{m,n}$ and $c_{m,l}$
on which they are mapped will also be exchanged; it follows that, no matter
how close the two starting functions are picked, one wavelet coefficient of
their image will differ by the value $c_{m,n}-c_{m,l}$; it follows that,
strictly speaking, the WQN algorithm is not continuous on any Besov or Sobolev
space. However, this theoretical argument does not necessarily constitute a
drawback in applications if we make the assumption that the coefficient
distribution for the reference (unaltered) data is continuous and is sampled
with enough precision; indeed, in practice consecutive wavelet coefficients
will be mapped to very close values, and the discontinuity of the mapping
would be of no practical consequence, since its “jumps” would be of very small
size. Note that this assumption is satisfied by the data we consider since the
histograms of wavelet coefficients follows a generalized Gaussian distribution
(Fig. 2). At this point, another phenomenon has to be taken into account: even
if the repartition function of coefficients for the reference data is
continuous (so that the “target” coefficients are very close), it is still
possible that the altered signal has a large number of coefficients which are
close to each other, so that a very small perturbation in the size of
coefficients would exchange two coefficients of very different ranks. This
argument shows that a continuity result for the functional operator underlying
the algorithm cannot follow from making only the assumption of a regularity of
the target probability density function (PDF). However, though such situations
yield mathematical counterexamples, they are not met in practice, since the
data on which the algorithm is applied also display smooth PDFs (Fig. 2), and
therefore do not exhibit large clusters of coefficients taking almost the same
value.

##  4 WQN algorithm on pathological examples

Figure 6:  Comparison of WQN and thresholding methods on pathological
examples. A. Reference coefficients perturbed by additive square artifact
$x_{\text{sq}}$. The WQN remapping causes swapping of coefficients, while
thresholding methods act locally. B. Signal with additive triangle
perturbation $x_{\text{tr}}$. C. Signal with a smooth perturbation
$x_{\text{cos}}$. Non-local effects of WQN are similar to the triangle case.

To further highlight the properties of the WQN algorithm, we present its
application to pathological examples and compare its behaviour with soft (eq.
12) and hard (eq. 11) thresholding. We initially study how the algorithm
operates on a single level of wavelet coefficients to highlight the non-local
action of the WQN, which we previously discussed in section 3.5. We consider
as a reference signal (in the wavelet space) a ramp $c_{\text{ref}}(t)=t$, as
the simplest signal having a simple invertible CDF. We then perturb the
wavelet coefficients by adding three types of artifacts to the reference
coefficients: a square artifact $c_{\text{sq}}(t)=(H(t-t_{1})-H(t-t_{2}))$
where $H$ is the Heaviside step function, a triangle

$c_{\text{tr}}(t)=\begin{dcases}\frac{2}{t_{2}-t_{1}}t,&\text{ if }\,t_{1}\leq
t\leq\frac{t_{1}+t_{2}}{2},\\\ 1-\frac{2}{t_{2}-t_{1}},&\text{ if
}\,\frac{t_{1}+t_{2}}{2}\leq t\leq t_{2},\\\ 0,&\text{
otherwise,}\end{dcases}$

and a cosine artifact

$c_{\text{cos}}(t)=\begin{dcases}\frac{1+\cos\left(\pi(2(t-t_{1})/T-1)\right)}{2},&\text{
if }\,t_{1}\leq t\leq t_{2},\\\ 0,&\text{ otherwise,}\end{dcases}$

with $T=t_{2}-t_{1}$. Considering
$c_{\text{art}}\in\\{c_{\text{sq}},c_{\text{tr}}\,c_{\text{cos}}\\}$, the
final artifacted coefficients are given by

$c(t)=c_{\text{ref}}(t)+c_{\text{art}}(t).$ (13)

In fig. 6 we show the correction of coefficients $c(t)$ when applying the WQN
quantile remapping and compare it with soft and hard thresholding. In the case
of the square artifact (fig. 6.A), due to the non-locality of the WQN
algorithm (see section 3.5), the mapping of the coefficients performed by WQN
throught the inverse CDF (quantile function) can swap the values of group of
coefficients (compare fig. 6.A first and second column). For wavelet
thresholding, we consider the threshold to be equal to the maximum value of
the unperturbed coefficients (which equals 1 in these examples). In the case
of the square artifact, this results in perfect isolation of the artifacted
coefficients, which are set to zero and 1 by hard and soft thresholding
respectively. Interestingly, while the thresholding methods act locally, they
do not result in a lower mean squared error with respect to WQN. The cases of
the triangle and cosine perturbation are presented in fig. 6.B and C
respectively. In both cases, the WQN remapping still shows non-local swapping
of coefficients, although the quantile function is continuous. Coefficients
with the same amplitude exist in the unperturbed ($t<t_{1},t>t_{2}$) and
perturbed ($t_{1}<t<t_{2}$) intervals, thus thresholding methods cannot
perfectly isolate the perturbed coefficients. We note that the difference in
smoothness between the triangle and cosine artifacts has no effect on the non-
local action of the WQN algorithm, which produces similar reconstructions in
the two cases. However, in both cases thresholding methods do not result in a
lower mean squared error with respect to WQN.

As a final example, we evaluated how the WQN algorithm performs on a sine wave
perturbed by white noise. This toy example was used to mimic a rhythmic brain
activity perturbed by random artifact such as an epileptic seizure. The
artifacted signal is given here as the sum

$x_{\text{art}}(t)=\sin(\omega t)+\sigma w(t),$ (14)

where $w$ is a centered Brownian noise of unit variance. In fig. 7.A we show
the signal reconstructed by WQN as we vary the noise amplitude
$\sigma=1,2,5,10$. This example shows that WQN allows a robust reconstruction
of rhythmic signals even in the presence of strong perturbations
($\sigma=10$). We quantify the noise reduction achieved by WQN in fig. 7.B.

Figure 7:  WQN reconstruction of a sine wave perturbed by additive white
noise. A. Reconstruction examples for $\sigma=1,2,5,10$. B. Mean squared error
(MSE) of the artifacted and the WQN-reconstructed signal with respect to the
original sine wave, i.e. $\operatorname*{\mathbb{E}}[(\sin(\omega
t)-x(t))^{2}]$, plotted versus the artifactual noise amplitude $\sigma$.

##  5 Conclusion and final remarks

We explored here the properties of EEG signals for different classes of brain
states and artifacts. We presented evidence that, for particular states such
as sleep and general anaesthesia, EEGs are not in general fully scale
invariant; the deviations are characterized by alternating dominance of
specific frequency bands. Moreover, while the wavelet coefficients statistics
at different scales are well approximated by Gaussian distributions in the
case of clean EEG, those of common artifactual signals such as EOG and EMG are
described by generalized Gaussians.

We then studied the properties of the WQN algorithm for EEG artifact removal.
We showed how the WQN algorithm transports the wavelet coefficients from the
distribution of an artifacted EEG signal into a distribution compatible with
the uncontaminated EEG, allowing to restore the signal statistics. We found
here that the WQN algorithm smoothens the discontinuities introduced by
artifactual signals and provided insight into its regularity, continuity, and
boundedness properties. We also highlighted how the remapping of wavelet
coefficients through the quantile function can produce non-local effects, as
opposed to traditional thresholding methods which always operate locally on
the wavelet coefficients. Indeed, WQN can transport a wavelet coefficient to
another position in time possibly located far away from the original position.
This effect can thus generate local distortion and some information contained
in the correlation structure of the wavelet coefficients can be lost. The
classical functional space such as Besov spaces are unable to “detect”
correlations between locations of wavelet coefficients, but other functional
spaces, such as “oscillations spaces” can be more appropriate to detect them.
In particular, they are not invariant under the “shuffling” of wavelet
coefficients [36]. Further investigations are needed in that direction.

Lastly, pathological cases can arise where the unperturbed physiological
signal exhibits irregularities such as spikes and waves. For instance, this is
particularly evident in certain forms of epilepsy or during deep anesthesia
[37]. When an artifact occurs in correspondence of these dynamics, it would be
interesting to test how WQN algorithm can map a singular signal into a
reference that also contains singularities. It is of interest to evaluate the
effectiveness of the WQN algorithm when an artifact occurs in correspondence
of these dynamics, analyzing the ability of the WQN algorithm in effectively
mapping a singular signal, affected by artifacts, to a reference signal that
likewise exhibits singularities. Given its clinical relevance, such analysis
warrants further exploration.

## Code and data availability

EEG, EMG, and EOG data used in the present article can be obtained from
publicly available datasets [22, 23, 24, 25, 28]. The Python code reproducing
all results and figures is available on Zenodo
(https://doi.org/10.5281/zenodo.8127712).

## Competing interests

The Authors declare no Competing Financial Interests.

## References

  * [1] D. L. Schomer, F. L. Da Silva, Niedermeyer’s electroencephalography: basic principles, clinical applications, and related fields, Lippincott Williams & Wilkins, 2012.
  * [2] E. R. Kandel, J. H. Schwartz, T. M. Jessell, S. Siegelbaum, A. J. Hudspeth, S. Mack, et al., Principles of neural science, Vol. 4, McGraw-hill New York, 2000\.
  * [3] M. Le Van Quyen, A. Bragin, Analysis of dynamic brain oscillations: methodological advances, Trends in neurosciences 30 (7) (2007) 365–373.
  * [4] P. Flandrin, M. Amin, S. McLaughlin, B. Torrésani, Time-frequency analysis and applications, IEEE signal processing magazine 30 (6) (2013) 19.
  * [5] D. L. Donoho, J. M. Johnstone, Ideal spatial adaptation by wavelet shrinkage, biometrika 81 (3) (1994) 425–455.
  * [6] I. M. Johnstone, B. W. Silverman, Wavelet threshold estimators for data with correlated noise, Journal of the royal statistical society: series B (statistical methodology) 59 (2) (1997) 319–351.
  * [7] V. Krishnaveni, S. Jayaraman, L. Anitha, K. Ramadoss, Removal of ocular artifacts from EEG using adaptive thresholding of wavelet coefficients, Journal of Neural Engineering 3 (4) (2006) 338–346.  doi:10.1088/1741-2560/3/4/011.
URL https://iopscience.iop.org/article/10.1088/1741-2560/3/4/011

  * [8] G. Inuso, F. La Foresta, N. Mammone, F. C. Morabito, Wavelet-ICA methodology for efficient artifact removal from Electroencephalographic recordings, in: 2007 International Joint Conference on Neural Networks, IEEE, 2007, pp. 1524–1529.  doi:10.1109/IJCNN.2007.4371184.
URL http://ieeexplore.ieee.org/document/4371184/

  * [9] M. Chavez, F. Grosselin, A. Bussalb, F. D. V. Fallani, X. Navarro-Sune, Surrogate-based artifact removal from single-channel eeg, IEEE transactions on neural systems and rehabilitation engineering 26 (3) (2018) 540–550.
  * [10] M. Dora, D. Holcman, Adaptive single-channel eeg artifact removal for real-time clinical monitoring, IEEE Transactions on Neural Systems and Rehabilitation Engineering (2022).
  * [11] M. Dora, S. Jaffard, D. Holcman, The wqn algorithm to adaptively correct artifacts in the eeg signal, Applied and Computational Harmonic Analysis 61 (2022) 347–356.
  * [12] S. Jaffard, Beyond Besov spaces - part 1: Distributions of wavelet coefficients, Journal of Fourier Analysis and Applications 10 (3) (2004) 221–246.  doi:10.1007/s00041-004-0946-z.
  * [13] S. Jaffard, Y. Meyer, R. D. Ryan, Wavelets: tools for science and technology, SIAM, 2001.
  * [14] P. Ciuciu, P. Abry, C. Rabrait, H. Wendt, Log wavelet leaders cumulant based multifractal analysis of evi fmri time series: Evidence of scaling in ongoing and evoked brain activity, IEEE Journal of Selected Topics in Signal Processing 2 (6) (2008) 929–943.  doi:10.1007/s00041-004-0946-z.
  * [15] M. V. Tsodyks, H. Markram, The neural code between neocortical pyramidal neurons depends on neurotransmitter release probability, Proc. Natl. Acad. Sci. USA 94 (1997) 719–723.
  * [16] G. Buzsáki, K. Mizuseki, The log-dynamic brain: how skewed distributions affect network operations, Nature Reviews Neuroscience 15 (4) (2014) 264–278.
  * [17] S. Jaffard, P. Abry, S. Roux, Function spaces vs. scaling functions: tools for image classification, Mathematical Image processing (Springer Proceedings in Mathematics) M. Bergounioux ed. 5 (2011) 1–39.
  * [18] P. Abry, H. Wendt, S. Jaffard, G. Didier, Multivariate scale-free temporal dynamics: From spectral (Fourier) to fractal (wavelet) analysis, Comptes Rendus de l’Académie des Sciences 20 (5) (2019) 489–501.
  * [19] I. Daubechies, Ten lectures on wavelets, SIAM, 1992.
  * [20] Y. Meyer, Ondelettes et Opérateurs, Hermann, Paris, 1990, english translation, Wavelets and operators, Cambridge University Press, 1992.
  * [21] S. G. Mallat, A theory for multiresolution signal decomposition: the wavelet representation, IEEE transactions on pattern analysis and machine intelligence 11 (7) (1989) 674–693.
  * [22] B. Blankertz, G. Dornhege, M. Krauledat, K.-R. Müller, G. Curio, The non-invasive berlin brain–computer interface: fast acquisition of effective performance in untrained subjects, NeuroImage 37 (2) (2007) 539–550.
  * [23] K. T. Sweeney, H. Ayaz, T. E. Ward, M. Izzetoglu, S. F. McLoone, B. Onaral, A methodology for validating artifact removal techniques for physiological signals, IEEE transactions on information technology in biomedicine 16 (5) (2012) 918–926.
  * [24] A. L. Goldberger, L. A. Amaral, L. Glass, J. M. Hausdorff, P. C. Ivanov, R. G. Mark, J. E. Mietus, G. B. Moody, C.-K. Peng, H. E. Stanley, Physiobank, physiotoolkit, and physionet: components of a new research resource for complex physiologic signals, circulation 101 (23) (2000) e215–e220.
  * [25] S. F. Quan, B. V. Howard, C. Iber, J. P. Kiley, F. J. Nieto, G. T. O’Connor, D. M. Rapoport, S. Redline, J. Robbins, J. M. Samet, et al., The sleep heart health study: design, rationale, and methods, Sleep 20 (12) (1997) 1077–1085.
  * [26] B. Hille, Ionic channels in excitable membranes. current problems and biophysical approaches, Biophysical Journal 22 (2) (1978) 283–294.
  * [27] G. Buzsaki, Rhythms of the Brain, Oxford University Press, 2006.
  * [28] H. Zhang, M. Zhao, C. Wei, D. Mantini, Z. Li, Q. Liu, Eegdenoisenet: A benchmark dataset for deep learning solutions of eeg denoising, arXiv preprint arXiv:2009.11662 (2020).
  * [29] W. O. Tatum, B. A. Dworetzky, W. D. Freeman, D. L. Schomer, Artifact: recording eeg in special care units, Journal of Clinical Neurophysiology 28 (3) (2011) 264–277.
  * [30] K. T. Sweeney, T. E. Ward, S. F. McLoone, Artifact removal in physiological signals—practices and possibilities, IEEE transactions on information technology in biomedicine 16 (3) (2012) 488–500.
  * [31] A. Antoniadis, Wavelet methods in statistics: some recent developments and their applications, Statistics Surveys 1 (none) (2007) 16 – 55.  doi:10.1214/07-SS014.
URL https://doi.org/10.1214/07-SS014

  * [32] A. Chambolle, R. A. De Vore, N.-Y. Lee, B. J. Lucier, Nonlinear wavelet image processing: variational problems, compression, and noise removal through wavelet shrinkage, IEEE Transactions on Image Processing 7 (3) (1998) 319–335.
  * [33] T. Pang, Z. Shen, Sparse estimation: An mmse approach, Constructive Approximation 57 (2) (2023) 727–758.
  * [34] J. Aubry, F. Bastin, S. Dispa, S. Jaffard, Topological properties of the sequence spaces s-nu, Journal of Mathematical Analysis and Applications 321 (1) (2006) 364–387.  doi:10.1016/j.jmaa.2005.08.036.
  * [35] S. G. Mallat, A theory for multiresolution signal decomposition: the wavelet representation, in: Fundamental Papers in Wavelet Theory, Princeton University Press, 2009, pp. 494–513.
  * [36] S. Jaffard, Beyond besov spaces, part 2: oscillation spaces, Constructive approximation 21 (2004) 29–61.
  * [37] A. Rigouzzo, L. Khoy-Ear, D. Laude, N. Louvet, M.-L. Moutard, N. Sabourdin, I. Constant, Eeg profiles during general anesthesia in children: A comparative study between sevoflurane and propofol, Pediatric Anesthesia 29 (3) (2019) 250–257.