Antoine Picot

Sylvie Charbonnier

Alice Caplier

Monitoring drowsiness on-line using a single encephalographic channel

Introduction

Drowsiness is the transition state between awakening and sleep during which a decrease of vigilance is generally observed. This can be a serious problem for tasks that need a sustained attention, such as driving. According to a report of the American National Highway Safety Traffic Administration [START_REF] Royal | National survey of distracted and drowsy driving attitudes and behavior[END_REF] driver drowsiness is annually responsible for about 56,000 crashes which is the reason why more and more researches have been developed to build automatic detectors of this dangerous state. Both behavioural and physiological modifications occur during drowsiness. Reaction time is slower, vigilance is reduced and information processing is less efficient, which can generate abnormal driving. Moreover, as drowsiness is the transition between awakening and sleep, it induces an increase of the number and the duration of blinks and yawns. Changes in cerebral activity also happen and can be observed thanks to electroencephalography. Researches on driver state monitoring has begun about thirty years ago and are still very active. The driver state monitoring systems can be classified into three kinds of system: those focusing on the vehicle behaviour, those focusing on the driver physical behaviour and those focusing on the driver physiological behaviour. The first systems developed were the ones using sensors monitoring the vehicle behaviour [START_REF] O'hanlon | A psychophysiological evaluation of devices for preventing lane drift and run-off-road accidents[END_REF][START_REF] Klein | Relationship between truck ride quality and safety of operations: methodology development[END_REF]. The main features studied are the steering wheel movements, the lateral position of the car on the road, the standard deviation of lateral position (SDLP) and the time to line crossing (TLC). The purpose is to detect an abnormal behaviour of the car, due to the driver drowsiness. The problems encountered by this kind of methods are that the features used depend on the shape of the road and how one drives, which may change a lot from one driver to another [START_REF] Renner | EEG-Based Drivers' Drowsiness Monitoring Using a Hierarchical Gaussian Mixture Model[END_REF].

To overcome these problems, researches have focused on systems using sensors monitoring drivers' awareness. One widespread technique to monitor the driver state is the use of a video camera. Indeed, a lot of information can be extracted from the driver face to monitor fatigue such as gaze, frequency and duration of eye blinking and yawning or percentage of eyelid closure. A lot of examples using camera to monitor the driver state can be found in the literature [START_REF] Grace | A drowsy driver detection system for heavy vehicles[END_REF][START_REF] Ji | Real time visual cues extraction for monitoring driver vigilance[END_REF][START_REF] Vural | Drowsy driver detection through facial movement analysis[END_REF]. These kinds of systems focus on the drivers' visual attention. Face, mouth and eye tracking algorithms are used to detect the face. Once the face, the eyes and the mouth are located, it is easy to detect eye blinking or yawning and calculate their frequency and duration. Frequency and duration of yawning or eye blinking too high indicate a decrease of attention. The gaze can be calculated with the eyes and the face position [START_REF] Smith | Determining driver visual attention with one camera[END_REF] or using a stereoscopic camera [START_REF] Ji | Real time non-intrusive monitoring and prediction of driver fatigue[END_REF]. Then, it allows the driver to be warned when he is not looking at the road. However, many differences can be observed between drivers, which makes it hard to monitor fatigue with only one feature [START_REF] Karrer | The role of individual differences in driver fatigue prediction[END_REF]). An interesting way of merging the different features (eye blinking, yawning, gaze...) is used by [START_REF] Ji | A probabilistic framework for modeling and real-time monitoring human fatigue[END_REF]. They use probabilistic networks which allow all features to contribute to the decision of the level of attention. Moreover, external factors (weather, hour of the day, etc...) can contribute in these networks to determine the level of attention. However, video features are not the best indicators of drowsiness. According to Dinges [START_REF] Dinges | An overview of sleepiness and accidents[END_REF], the best indicators of fatigue are the physiological indicators. The electroencephalogram (EEG) and the electro-oculogram (EOG) are mainly used to study drowsiness. Yet, several researches have focused on other physiological indicators such as the electrocardiogram (ECG) to monitor drivers' heart rate [START_REF] Törnros | Heart rate measures as drowsiness indicators[END_REF] or the drivers' temperature [START_REF] Quanten | Thermoregulatory changes at driver sleepiness[END_REF]. The EOG is the measurement of the resting potential of the retina. It gives an accurate measurement of eyes movements. Many features can be extracted from this information such as eyelid opening and closing parameters, blinks frequency, blinks amplitude, blinks duration... According to [START_REF] Galley | Blink parameter as indicators of driver's sleepiness -possibilities and limitations[END_REF], EOG is a relevant measure to monitor fatigue since some extracted features are really sensitive to drowsiness. One of the most efficient features extracted is the PERCLOS (PERcentage of eyelid CLOSure). This feature has been defined by Wierwille [START_REF] Wierwille | Research on vehicle-based driver status/performance monitoring[END_REF]. It is the percentage of eyelid closure over the time. Knippling (Knippling, 1998) showed that PERCLOS is a good indicator of drowsiness that increases with fatigue. Electroencephalography measures the electrical activity of the brain from electrodes placed over the scalp. Drowsiness appears into the EEG spectrum by an increase of activity in the frequency bands [START_REF] Hanley | The meaning and use of the area under a receiver operating characteristic (roc) curve[END_REF][9][10][11][12]Hz (alpha band) and [4][5][START_REF] References Akerstedt | Subjective and objective sleepiness in the active individual[END_REF][7][START_REF] Hanley | The meaning and use of the area under a receiver operating characteristic (roc) curve[END_REF]Hz (theta band) predominantly in the parietal and central regions of the brain. In the same time, a decrease of activity in the band [12][13][14][15][16][17][18][19][20][21][22][23][24][25][26]Hz (beta activity) can also be observed, as beta activity increases with cognitive tasks and active concentration. This has been shown in several studies (Santamaria & Chiapa, 1987;[START_REF] References Akerstedt | Subjective and objective sleepiness in the active individual[END_REF][START_REF] Kay | Changes in airway resistance during sleep onset[END_REF]. EEG is so efficient in detecting drowsiness that it is often used as a reference indicator. In this case, the reference is built by expert doctors who visually observe the proportion of alpha and theta activity on a short-time window as in [START_REF] Muzet | Awake pilot study no.2: Testing steering grip sensor measures[END_REF]. The analysis is done off-line and is time consuming. In order to make the analysis of drowsiness in an automatic way, the EEG power spectrum can be computed using Fast Fourier Transform or using wavelets Transform but none of these techniques seems better than the other. The number of EEG channels used to monitor drowsiness fluctuates from a few to about thirty. The advantage of using a large number of EEG channels is to obtain spatial information on how the EEG energy is shifting from one frequency band to another [START_REF] Makeig | Independent component analysis of electroencephalographic data[END_REF][START_REF] Lin | Estimating driving performance based on eeg spectrum analysis[END_REF]. However, using only a few EEG channels is faster and easier to compute. The features obtained from EEG can be used in many ways. Several studies proposed to monitor some ratios between different EEG power bands. De Waard and Brookhuis (De Waard & Brookhuis, 1991) suggest to monitor the ratio between the alpha activity and the theta activity on the beta activity [(alpha+theta)/beta] (beta band: [12-26]Hz) which increases with driver drowsiness. Drowsiness decision can be made using EEG features used as inputs to a classifier. Ben [START_REF] Ben Khalifa | Alertness states classification by SOM and LVQ neural networks[END_REF]) used a connexionist method to detect drowsiness. They proposed to use a single parieto-occipital channel (P4-O2) to avoid ocular frontal artefacts. Relative power in the [1-23]Hz band is computed using a short term fast Fourier transformation and divided in twenty three sub-bands of 1Hz. These features are then used as entries in a linear vector quantization (LVQ) neural network to detect artefacted states and in a self organizing map (SOM) neural network to classify states into two states: awake and drowsy. Lin et al. (Lin et al., 2005a) applied Principal Component Analysis (PCA) on the whole spectrum of a 2-channels EEG. A fast Fourier transform is then computed on the principal component to extract the power in 10 EEG bands from 1 to 40Hz. These features are used for training a linear regression model which estimates the driving performance, i.e. the deviation between the car and the centre of the road, assuming that driving performances are correlated with driver's alertness. They also used the two principal components obtained using an independent component analysis (ICA) on a 33-channels EEG instead using a PCA on a 2-channels EEG with their linear regression model to improve their estimation of driving performance. Most recently, Rosipal et al. (Rosipal et al., 2007) used hidden Gaussian mixture model to monitor drowsiness. A hierarchical Gaussian mixture model (hGMM) with two mixture components at the lower hierarchical level is used. Each mixture models the data density distribution of one of the two drowsiness cornerstones/classes represented by 4-second long EEG segments with low and high drowsiness levels. The spectral content of each EEG segment is transferred into a compact form of autoregressive model coefficients. Their study is performed on a large number of drivers. However, even if some on-line techniques exist, they seem hard to implement: they either need to be trained on numerous data previously recorded on a driver or they require a large number of EEG channels. Yet, the fact that EEG is used as a drowsiness reference by expert doctors in physiological studies proves the utility of an on-line driver drowsiness detection system using EEG. In this paper, we propose to develop a system to detect driver drowsiness that uses as few EEG channels as possible that works on-line and that does not need to be trained. Indeed, it seems that a relatively simple drowsiness detection system for drivers would be easier to implement in a car. The method proposed is described in section 2. The results obtained on a significant database of 40 EEG recordings from 20 drowsy drivers are shown and discussed in section 3. Finally, conclusions and perspectives are presented.

Drowsiness detection method

After a short presentation on EEG and how drowsiness appears in EEG, the drowsiness detection method is presented step by step.

EEG and drowsiness

An electroencephalogram is a measurement of the electrical activity of the brain from electrodes placed on the scalp [START_REF] Blinowska | Electroencephalography(EEG)[END_REF], according to the international 10-20 system shown in fig. 1. Drowsiness is characterized by an increase of alpha and theta activities, predominantly in the parietal (P) and central (C) regions of the brain, and a slowdown of blinks and eye movements [START_REF] References Akerstedt | Subjective and objective sleepiness in the active individual[END_REF][START_REF] Kay | Changes in airway resistance during sleep onset[END_REF][START_REF] Muzet | Awake pilot study no.2: Testing steering grip sensor measures[END_REF]. Even if different scales of drowsiness classification exist, none of them are standardized and there are no standardized rules to differentiate the levels of drowsiness (as the Rechtschaffen and Kales rules [START_REF] Rechtschaffen | A Manual of Standardized Terminology, Techniques and Scoring System for Sleep Stages of Human Subject[END_REF]) for the study of sleep). This may be due to the quite recent interest on drowsiness compared to the sleep analysis and the difficulty to collect drowsiness data. More than 10s Slow Table1. OSS Criteria There are two kinds of scales: subjective sleepiness scales like the Karolinska Sleepiness Scale (KSS) [START_REF] References Akerstedt | Subjective and objective sleepiness in the active individual[END_REF] which allows drivers to directly evaluate their own drowsiness and Objective Sleepiness Scales (OSS) which is used by expert doctors to evaluate drivers' drowsiness after driving. The OSS used in this study is the five levels scale from 0 (awake) to 4 (very drowsy) developed by [START_REF] Muzet | Awake pilot study no.2: Testing steering grip sensor measures[END_REF]. Decisions are made every 20s and depend on the length of alpha and theta bursts as well as on the speed of eyes movements and blinks. The different criteria are presented in table 1. This scale is used by doctors to make decisions.

Objective

Method Principle

The purpose is to design a drowsiness detection algorithm which can work on-line, inspired by the OSS. The overview of the detection method is shown in fig. 2. First the EEG power spectrum is computed using a Short Time Fourier Transform (STFT) to calculate the relative power into the different EEG bands every second. Then, the relative power of the alpha band is median filtered using a sliding window to reject abnormal values. A Means Comparison Test (MCT) is computed at last to compare the energy to a reference level, learnt at the beginning of the recording while the patient is not supposed to be drowsy. MCT is normalized. A common threshold of detection can be proposed taking into account the acceptable level of false alarms and validated using experiments which has been presented in [START_REF] Picot | On-line automatic detection of driver drowsiness using a single EEG channel[END_REF]. Concomitantly, a Variances Comparison Test (VCT) is computed on the raw EEG data to detect high amplitude artefacts. Information on the occurrence of artefacts can be used as an index of reliability on the "drowsy decision".

EEG Power Spectrum

The EEG power spectrum is computed using a Short Time Fourier Transform (STFT). The power spectrum is computed every second on a window of two seconds using Welch's periodogram method [START_REF] Welch | The use of fast fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms[END_REF]. The overlapping window between the previous and next value is 1 second. Then, the relative powers in each band are calculated as the ratio of the power in one band and the power of the whole EEG spectrum. Only the range [1-30]Hz is used because activity below or above this range is likely to be artefactual. For example the α relative power is calculated as follows:

Median filtering

Median filtering is used to smooth the αrel signal and to reject abnormal values. The median is the value separating the higher half of a population from the lower half. Compared to the mean, the median is known for its robustness towards outliers, as far as the number of outliers is lower than half the length of the population. Here, the median of the relative powers is calculated every second, before performing MCT, using a sliding window of 10s.

Means Comparison Test

The method of MCT is inspired by Ragot [START_REF] Ragot | Validation de données et diagnostic[END_REF] and is applied on the relative powers in the alpha band. A moving window is compared to a fixed reference window as shown in fig. 3. x and whose variances are s1 2 and s2 2 . Then, the variable:

        + + + = 2 1 2 1 2 2 2 2 1 1 2 1 n 1 n 1 2 - n n s n s n x - x t (2) follows a 2 - n n 2 1 +
liberty degrees Student law. The equality of the two means can be tested by a bilateral test with a confidence threshold λ: -t1λ/2 < t < t1λ/2 . If the two populations have the same length n (i.e. n1=n2=n) and their theoretical variances are equal, (2) can be formulated as:

n s s x - x t 2 2 2 1 2 1 + = (3)
So, the variable t follows a n-1 liberty degrees Student law. Moreover, if the populations are large enough, i.e. if n1 and n2 are equal or greater than 20, the test is performed with the variable:

2 2 2 1 2 1 2 1 n s n s x - x u + = (4)
which then follows a centred reduced normal law. The means equality is then tested by a bilateral test with a confidence threshold λ: -u1λ/2 < u < u1λ/2 . Furthermore, the theoretical variances equality is no longer necessary.

Here, the test is computed on the relative power signals calculated every second. The reference, calculated from the fixed window, is calculated at the beginning of a recording with the assumption that, before driving, the driver is not drowsy. The mean value calculated during this period provides a reference for the not drowsy stage. The mean values calculated on line are compared to this "not drowsy" reference value. The length of the fixed reference window is n1=60s and the one of the moving window is n2=30s, so n1 and n2 are greater than 20 as αrel is calculated every second. Hence, the formula (4) can be used.

The threshold λ fixes the percentage of false alarms expected. The higher the threshold, the lower the percentage of false alarms. In this study, λ is empirically chosen and discussed in section 3.

Artefact detection

An artefact is an electric perturbation of the EEG signal due to patient movements or measurement problems. Artefacts pollute the whole EEG band and it is quite impossible to extract reliable EEG information when an artefact occurs. There are several types of artefacts. They may be due to ocular movements, face muscles movements or measurement devices problems such as electrode unstuck The method of VCT is inspired by Ragot [START_REF] Ragot | Validation de données et diagnostic[END_REF] and is directly applied on the raw EEG. The principle is the same than the MCT: the variance of a moving window is compared to the variance of a reference window. Let us consider two independent populations with normal distributions. Their lengths are n1 and n2 and their experimental variances are s1 2 and s2 2 . Then, the variable:

2 2 2 1 s s F = (5)
follows a Fisher law with k1=n1-1 and k2=n2-1 liberty degrees. The variances equality is then tested by a bilateral test with a confidence threshold λart: Here, the test is computed on the raw EEG data every 10s. A moving window is compared to a fixed reference window as shown in fig. 3. The reference variance is calculated on an artefact-free time window lasting one minute, chosen at the beginning of the recording,.

F λart /2 < F < F1-λart /2.
Then, a moving window of 10s is compared to this reference every 10s. The threshold λart is empirically chosen and discussed in section 3.

The goal is to detect the artefact, not to reject it. As the system is working on a minimal number of EEG channels, it is not possible to recover lost EEG information of the artefact since missing information cannot be found somewhere else. Nevertheless, detecting the occurrence of an artefact provides information on the signal quality: whenever an artefact is detected, the concomitant relative power extracted from EEG should not be used to evaluate drowsiness.

Method relevance

The whole algorithm can be applied on-line. However, the sliding window of 10s used for median filtering induces a delay of 5s and the sliding window of 30s used for the MCT induces a delay of 15s. The artefact detection is computed in parallel with a sliding window of 10s which induces a delay of 5s. So, the decision provided by the algorithm is delayed by 20s from the signals recording. This latency in the decision will be taken into account when comparing the results to the expert's decisions.

The general purpose of this algorithm is the detection of drowsiness. The MCT detects α bursts, which are indicators of drowsiness as seen in section 1. The reference is calculated on a fixed window chosen at the beginning of the signal, supposing that the driver is completely awake when he starts driving. So, the mean calculated on the moving window is compared to a wakefulness reference. If the bilateral test is higher than the threshold, the driver is then considered as drowsy, otherwise he is considered as awake. Fig. 6 shows how the signal is processed in the detection system. First the relative power in the α band (b) of the EEG (a) is calculated. Then, it is smoothed by median filtering (c). A MCT is performed (d) on the filtered signal and is thresholded to make the decision awake or drowsy (e). Moreover, a VCT is calculated on the raw EEG signal to detect the occurrence of highamplitude artefacts polluting the whole EEG band. The detection of these high-amplitude artefacts does not allow rejecting them but provides information on the quality of signal around this point. It means that if artefacts are found on a part of the signal, decisions on drowsiness in this part tend to be less reliable than if not.

The point with detecting α bursts in EEG signal is the difficulty to define a common threshold for all drivers because of the large inter-individual differences [START_REF] Karrer | The role of individual differences in driver fatigue prediction[END_REF].

Here, the level of αrel power in the "awake" state is learned on each driver from the reference window. Moreover, the output of MCT is a variable following a centred reduced normal law. So, the threshold used in the bilateral test has statistical meaning and is the same for all drivers.

In the same way, as the output of the VCT is a variable following a Fisher distribution, the threshold used to detect high-amplitude artefacts has a statistical meaning and is the same for all drivers.

Results and discussion

Database

The database used for the evaluation of the method has been provided by the CEPA (Centre The method proposed in this chapter provides a binary decision [awake; drowsy] while the database has been manually labelled using five levels. Moreover, the expert classified non overlapping intervals of 20s (epochs) while the automatic system makes a decision every second. To compare our results with the expert's decision, the following validation technique was used. The five expert decision levels were converted into a binary decision by considering as drowsy any decision superior or equal to 1 in the expert's scale as shown in fig. 9. This figure shows the expert decision on a five levels scale (a) and on a binary scale (b) and the drowsiness detection obtained using our system (c). Furthermore, each 20s epoch classified by the expert was directly compared to the system decision: if during the 20s interval, the system classifies at least 1s as "drowsy", then the decision for the epoch was "drowsy". Else it was "awake".

Technical validation

Epochs were then compared one by one and classified according to the table of confusion 2.

Expert decision

Automatic decision

awake drowsy awake True Negative (TN)

False Negative (FN) drowsy False Positives (FP)

True Positive (TP) Table 2

. Table of confusion

The true positive rate (TPrate) or detection rate is the ratio between the number of true "drowsy" automatic decisions and the number of "drowsy" expert decisions. The false positive rate (FPrate) is the ratio between the number of false "drowsy" automatic decisions and the number of "awake" expert decisions. They are calculated according to ( 6) and ( 7).

FN TP

TP TP rate + = (6) TN FP FP FP rate + = (7)
The results are displayed as Receiver Operating Characteristic (ROC) curves [START_REF] Hanley | The meaning and use of the area under a receiver operating characteristic (roc) curve[END_REF], plotting TPrate in function of FPrate. The purpose is to have the highest TPrate with the lowest FPrate.

Results using alpha relative power

Results without artefact detection

The drowsiness detection algorithm was applied on the whole database, with a decision threshold λ (defined in section 2.4) varying from 1.5 to 5, on each of the four EEG channels. The results presented in fig. 10 are those obtained when the MCT is applied on the alpha relative power without considering artefact detection. "Star" markers correspond to the P3 channel, "circle" markers to the F3 channel, "square" markers to the C3 channel and "triangle" markers to the O1 channel. The head at the bottom right of Fig. 10 reminds the reader of the position of each channel. For each channel, the results represented with the markers the further on the right corresponds to the smallest λ and those with the markers the further on the left to the biggest λ. It is coherent: increasing λ diminishes the FPrate while decreasing the TPrate.

Fig. 10. Results obtained using different EEG channels

It is obvious from fig. 10 that the results are better when the P3 parietal channel is used, which is in concordance with results from the literature: drowsiness is characterized by an increase of α activity predominately in the parietal region of the brain. Indeed, results obtained with EEG recorded by C3, F3 or O1 are only slightly better than those that would be obtained with a random classifier. In the following sections, only the results obtained on P3 channel are shown.

The optimal result on P3 are TPrate=82,1% and TPrate=19,2% with λ=3. However changing the threshold does degrade the performances (TPrate=85,1% and TPrate=23,5% with λ=1,5 and TPrate=76,9% and TPrate=14,8% with λ=5), which proves that the method is not sensitive to the threshold value. The last signal (c) shows the result of the artefact detection (dotted line): when "zero" no high amplitude artefact is detected and when "one", an artefact is detected. The dotted line boxes underline high-amplitude artefacts.

Results with artefact detection

First, this example shows that the detected artefacts correspond to high-amplitude electric perturbations of the EEG signal. As high-amplitude artefacts have not been evaluated by an expert on the dataset, it is not possible to quantify the performances of the artefact detection method. Nevertheless, a visual check of all the recordings shows that all the apparent highamplitude artefacts have been detected. Fig. 12 shows the number of artefacts detected in the database (total number and corresponding percentage on the database) in function of the value of the threshold λart (used for artefact detection). 12 that if the threshold is too small (λart<5), detection is really sensitive and a lot of points are rejected. Visually, this means a lot of false alarms. When increasing λart, the number of artefacts detected decreases quickly till λart=6, which visually seems an appropriate threshold. Indeed, for this threshold value, all the visible high-amplitude artefacts are detected without false alarms. At this point, one can see that high-amplitude artefacts represent only a small part of the dataset: about 2%.

Results using other features than alpha relative power

Results presented in section 3.3 are now compared to results obtained with other features proposed in the literature. The results are displayed in fig. 14. Results from section 3.3, obtained using MCT on the median filtered αrel signal, are represented by "star" markers. "Square" and "circle" markers represent results obtained using MCT and median filtering respectively on θrel and βrel signals. Note that β activity increases with cognitive tasks and active concentration, so drowsiness is characterized by a decrease of the β activity. So, the detection algorithm using β as the main feature consider the driver as "drowsy" when the output of the MCT is lower than the threshold -λ (varying from -5 to -1). "Triangle" markers correspond to results obtained with the combined signals αrel|θrel. Decisions are made independently on αrel and on θrel and then merged with a logical OR. The optimum threshold λ=3 was used for αrel. Displayed results are obtained with a threshold varying from 1,5 to 5 for θrel detection. The idea to use both αrel and θrel is inspired by table 1, where it is assumed that drowsiness is characterized by an increase of the activity in one of the two frequency bands α and θ. "Hexagram" markers represent the (α+θ)/β feature. This feature has been suggested by De Waard and Brookhuis (De Waard & Brookhuis, 1991). As α and θ activity are supposed to increase with drowsiness whereas β activity is supposed to decrease, this feature should be increasing with drowsiness. In this case, MCT is computed on the sum of the signals αrel and θrel divided by the βrel signal and only one threshold is used to make the decision.

Fig. 14. ROC curves using different features for drowsiness detection

The best results are obtained with the drowsiness detection algorithm applied on the αrel signal (TPrate=82,4%, FPrate=18,3%). Since the algorithm was tested with the same threshold on data recorded from 20 different persons, this tends to show that the method can be applied on different persons without adapting the tuning parameter.

The results obtained with αrel|θrel show that θrel is not relevant to detect drowsiness since the number of false positive increases tremendously when this information is added. This is confirmed by the results obtained with θrel only. The (α+θ)/β ratio and βrel give correct results (TPrate=76,2% and TPrate=32,1% for (α+θ)/β and TPrate=75,9% and TPrate=24,1% for βrel) but worse than the results obtained with the αrel information only.

Now, if we compare the results obtained with the literature, the results obtained are as good as those found when using a trained algorithm. Lin et al. (Lin et al., 2005a) proposed to monitor driving performance, i.e. the capacity to maintain the car in the middle of the road computing a linear regression model on a 2-channels EEG. They obtained a correlation of r=0,88 between their model and the driving performances when the model is trained and tested on the same session. The correlation decreases to r=0,7 when trained and tested on different sessions. So, this method needs to be tuned for each driver as the model estimated for one driver does not work so well on another. Lin et al. increase these results using ICA on a 33-channel EEG (Lin et al., 2005a) to compute their linear regression model and obtain a correlation of r=0,88 between their estimation and the driving performances on the testing session. Nevertheless, this model needs to be trained on a large amount of data and has been tested on 5 drivers only. Ben Khalifa et al. [START_REF] Ben Khalifa | Alertness states classification by SOM and LVQ neural networks[END_REF]) obtained 92% of true drowsiness detections by training a neural network on a the EEG spectrum of the P4-O2 channel but this result is obtained on the training set and decreases to 76% of true detections on the validation set. Moreover, these results are obtained on a set of only 4 drivers. At least, Rosipal et al. (Rosipal et al., 2007) obtained about 77% of true detections of drowsiness states by using hGMM on the spectral content of EEG transferred into a compact form of autoregressive model coefficients. This study has been performed on a large number of drivers and needs a period of training.

The advantage of the method proposed in this paper is that it does not need to be trained or adapted. The same threshold can be used for all drivers. Moreover, as the method has been tested on huge dataset, the results can be considered significant. Results are improved using this fuzzy logic approach since FPrate is increased and FPrate is decreased. The results obtained with this method are TPrate=84,6% and FPrate=17,9% (TPrate=82,4% and FPrate=18,3% with λ=3 when using only αrel information). This means that the βrel information is relevant to detect drowsiness when combined with αrel. Moreover, compared to the method proposed in section 3.5, there is no need to select an appropriate detection threshold for αrel and βrel. The fuzzy approach increased the detection reliability.

Conclusion

An algorithm to detect on line drivers' drowsiness from a P3 EEG channel has been presented in this paper. This algorithm is based on a means comparison test (MCT) applied on the EEG relative power calculated in the alpha band and in the beta band. The results of the MCT test are then merged using fuzzy logic. The algorithm can operate on-line with a short delay and does not need to be tuned. Performances obtained on a large data set recorded from 20 different drivers are 84,6% of true detection and 17,9% of false detection using the parietal EEG channel only. A high-amplitude artefact detection system has been developed and combined to the drowsiness detector. It enables periods of time when the EEG signal is unreliable to be detected on line. No decision is made by the drowsiness detector while the artefact detector classifies the EEG signal as unreliable. The artefact detector is tuned by a single threshold whose value is independent of the driver.

The next step of this work is to add an ``eye blinks and yawn'' detection system thanks to a high frame rate camera and to merge the decisions to obtain a highly reliable automatic drowsiness detector. Fuzzy logic could be a first step to merge this information.

Fig. 1 .

 1 Fig. 1. 10-20 system EEG is described in term of rhythmic activity and transients. The rhythmic activity is divided into frequency bands: delta (δ) activity ([0.5-4]Hz), theta (θ) activity ([4-8]Hz), alpha (α) activity ([8-12]Hz), beta (β) activity ([12-26]Hz) and gamma (γ) activity (over 26Hz). Most of the time, only the range [1-30]Hz is used because activity below or above this range is likely to be artefactual (under standard clinical recording techniques).Drowsiness is characterized by an increase of alpha and theta activities, predominantly in the parietal (P) and central (C) regions of the brain, and a slowdown of blinks and eye movements[START_REF] References Akerstedt | Subjective and objective sleepiness in the active individual[END_REF][START_REF] Kay | Changes in airway resistance during sleep onset[END_REF][START_REF] Muzet | Awake pilot study no.2: Testing steering grip sensor measures[END_REF]. Even if different scales of drowsiness classification exist, none of them are standardized and there are no standardized rules to differentiate the levels of drowsiness (as the Rechtschaffen and Kales rules[START_REF] Rechtschaffen | A Manual of Standardized Terminology, Techniques and Scoring System for Sleep Stages of Human Subject[END_REF]) for the study of sleep). This may be due to the quite recent interest on drowsiness compared to the sleep analysis and the difficulty to collect drowsiness data.

Fig. 2 .

 2 Fig. 2. Drowsiness detection method

  in the bands α ([8-12]Hz) , θ ([4-8]Hz) and β([12-26]Hz) are named αrel, θrel and βrel respectively.

Fig. 3 .

 3 Fig. 3. Illustration of the windows for the MCT The classical MCT has quite restrictive conditions due to the fact that the theoretical variances are unknown. Let us consider two independent populations of length n1 and n2, whose means are 1x and 2 x and whose variances are s1 2 and s2 2 . Then, the variable:

  . The proposed solution to detect high-amplitude artefacts corresponding to electrode unstuck is to apply a Variances Comparison Test (VCT) on the EEG. Artefact values are around mV whereas EEG signal values are around µV. High-amplitude artefacts pollute the whole EEG band as shown in the red box on the spectrogram in fig. 4.

Fig. 4 .

 4 Fig. 4. Example showing a high-amplitude artefact

Fig. 6 .

 6 Fig. 6. Signal processing from EEG to drowsiness detection High-amplitude artefacts pollute the EEG signals and generate isolated high abnormal values on the whole EEG band of the power spectrum. A median filter is used to smooth the α relative power signal and to reject abnormal isolated values to avoid false detection.

  d'Etudes de Physiologie Appliquée) from Strasbourg (France) using the driving simulator PAVCAS ("Poste d'Analyse de la Vigilance en Conduite Automobile Simulée"). PAVCAS is a moving base driving simulator composed of a mobile base with four liberty degrees (vertical and longitudinal movements, swaying and pitching) and a real-time interactive visualization unit. The visualisation unit reproduces very well the driving conditions on a freeway by day or night. Images are shown on five screens in front of the vehicle and are arranged in semicircle.The database is composed of forty recordings from twenty subjects. Each subject was recorded while driving for 90 minutes, a first time perfectly rested and a second time suffering from sleep deprivation (the subject had slept 4 hours only) in diurnal conditions. The database is thus composed of 60 hours of driving data. Each recording includes four EEG channels (left frontal (F3), central (C3), parietal (P3) and occipital (O1)), one EOG channel and a video of the driver's face. Objective sleepiness was evaluated on each recording by an expert doctor using the scale described in section I. Data acquisition of physiological signals was performed at 250Hz.

Fig. 9 .

 9 Fig. 9. Comparison between expert decision (a and b) and system decision (c)

  An example of artefact detection is shown on fig. 11. The first signal (a) is the EEG raw data. The signal (b) is the result of the VCT. The dotted line corresponds to the threshold λart=6.

Fig. 11 .

 11 Fig. 11. Example of high amplitude artefact detection It can be seen in fig.12 that if the threshold is too small (λart<5), detection is really sensitive and a lot of points are rejected. Visually, this means a lot of false alarms. When increasing λart, the number of artefacts detected decreases quickly till λart=6, which visually seems an

Fig. 12 .

 12 Fig. 12. Number of artefact detected in function of λart

Fig. 13 .

 13 Fig. 13. ROC curve for drowsiness detection when artefact detection is usedIt is obvious in fig.13that results are slightly improved: TPrate is a bit increased while FPrate is a bit decreased. Using the threshold λ=3, results increase from TPrate=82,1% and FPrate=19,2% to TPrate=82,4% and FPrate=18,3%. So, artefact detection decreases the number of false decisions. The fact that results are only slightly increased can be explained by the fact that high-amplitude artefacts represent about 2% only of the dataset. Artefact detection will be taken into account in the following sections.

Fig. 15 .

 15 Fig. 15. Experimental probabilities of being drowsy in function of threshold λ

Fig. 17 .

 17 Fig. 17. Performances obtained using αrel, βrel and (αrel, βrel) merged with fuzzy logic

Acknowledgements

The authors are grateful to the Centre d'Etudes de Physiologie Appliquée (CEPA) in Strasbourg (FR) for providing the data and their help, as well as the Laboratoire

Results merging alpha and beta relative powers

From the previous section, the best results are obtained when αrel or βrel are used as features, which naturally gives the idea to merge these two features to increase the decision reliability The technique used to merge αrel and brel is fuzzy logic, which is based on the theory of fuzzy sets developed by Zadeh [START_REF] Zadeh | Fuzzy sets[END_REF]. Let us consider µDr(αrel) and µDr(βrel), the membership functions, which represent the membership degree of αrel and βrel, independently considered, to the "drowsy" state. The purpose is to make a decision Dr(αrel, βrel) using both µDr(αrel) and µDr(βrel),. The driver is considered to be drowsy if both the decision made using αrel and the decision made using βrel is "drowsy". This is expressed thanks to the t-norm product as follows:

Note that µAw(.) is the membership function of the "awake" state and is the complementary of µDr(.). The denominator is used here to normalize Dr(αrel,βrel) between 0 and 1. One has to define the membership function µDr(αrel) and µDr(βrel). A study of the probability of being drowsy in function of the MCT's threshold λ on αrel and βrel, P(dr|αrel) and P(dr|βrel), is displayed in fig. 15. The "square" markers line displays the experimental P(dr|αrel) and the "circle" markers line displays the experimental P(dr|-βrel). Probabilities are calculated as the percentage of true drowsiness detections obtained on periods when the relative power is over λ. The membership function µDr (αrel) and µDr (βrel) are then designed from the results presented in fig. 15. As αrel and βrel have a very similar behaviour, the same membership is used for αrel and βrel. This membership function is presented in fig. 16. The driver is considered as "drowsy" when Dr(αrel,βrel) is larger than 0.5 The results obtained with this method are shown in fig. 17 with the "circle" markers. They are compared to the results obtained using the MCT on αrel only ("square" markers).


Antoine Picot

Sylvie Charbonnier

Alice Caplier

Drowsiness detection based on visual signs: blinking analysis based on high frame rate video

In this paper, an algorithm for drivers' drowsiness detection based on visual signs that can be extracted from the analysis of a high frame rate video is presented. A study of different visual features on a consistent database is proposed to evaluate their relevancy to detect drowsiness by data-mining. Then, an algorithm that merges the most relevant blink features (duration, percentage of eye closure, frequency of the blinks and amplitude-velocity ratio) using fuzzy logic is proposed. This algorithm has been tested on a huge dataset representing 60 hours of driving from 20 different drivers. The main advantage of this algorithm is that it is independent from the driver and it does not need to be tuned. Moreover, it provides good results with more than 80% of good detections of drowsy states.

I. INTRODUCTION

Drowsiness is the transition between awake state and sleep during which one's abilities to observe and analyse are strongly reduced. According to the National Highway Traffic Safety Administration (NHTSA), about 100,000 crashes are the direct result of driver drowsiness each year. This is the reason why more and more researches are made to build automatic detectors of this dangerous state.

According to Renner and Mehring [START_REF] Renner | Lane departure and drowsiness -two major accident causes -one safety system[END_REF], drowsiness can be detected both in brain activity which refers to the ability to process the information and in eyes activity which refers to the perception ability. This study focuses on the visual signs of drowsiness: blinks.

Electrooculogram (EOG), which is the measurement of the eye electrical muscles activity, has been widely used in the literature to estimate drowsiness ( [START_REF] Caffier | Experimental evaluation of eyeblink parameters as a drowsiness measure[END_REF], [START_REF] Galley | Blink parameter as indicators of driver's sleepiness -possibilities and limitations[END_REF]). EOG is the most reliable technique to detect and characterize blinks due to its high sample rate (from 250Hz to 500Hz) and it is used as a reference by expert doctors to evaluate drowsiness ( [START_REF] Gillberg | Sleepiness and performance of professional drivers in a truck simulator -comparisons between day and night driving[END_REF], [START_REF] Muzet | AWAKE pilot study no.2: Testing steering grip sensor measures[END_REF]). Unfortunately, EOG requires at least three electrodes placed on the driver's skin, which is not pleasant. So, for obvious ergonomic reasons, research community has focused on the use of video to track the driver's eyes and face and thus detect whether he is drowsy or not.

Several features can be extracted from the video to estimate drowsiness such as blink duration, blink frequency, the PERcentage of eye CLOSure (PERCLOS) as well as driver's gaze or facial expressions ( [START_REF] Knipling | PERCLOS: a valid psychophysiological measure of alertness as assessed by psychomotor vigilance[END_REF], [START_REF] Ji | Real-time eye, gaze, and face pose tracking for monitoring driver vigilance[END_REF], [START_REF] Boverie | Diagnostic fusion for in vehicle driver vigilance assessment[END_REF]). Most of these systems focus on too long blink durations to detect drowsiness ( [START_REF] Horng | Driver fatigue detection based on eye tracking and dynamk, template matching[END_REF]). Johns [START_REF] Johns | The amplitude-velocity ratio of blinks: a new method for monitoring drowsiness[END_REF] suggested to use the amplitude-velocity ratio as an universal drowsiness indicator. Most recently, Damousis et al. [START_REF] Damousis | A fuzzy expert system for the early warning of accidents due to driver hypo-vigilance[END_REF] proposed a fuzzy expert system merging blink duration, A recent study by Picot et al. [START_REF] Picot | Comparison between EOG and high frame rate camera for drowsiness detection[END_REF] has shown that the use of a high frame rate video (200f ps) can replace the EOG signal for the extraction of several blink features with the same accuracy. A drivers' drowsiness detection system based on these particular blinking features is proposed here. This system is presented on EOG because expertized high frame rate drowsiness video are still not available. Nevertheless, a huge dataset of 60 hours of driving EOG data is used here to develop this system. The database is first used to select the most relevant blinking features. A drowsiness detection method based on fuzzy logic is then proposed and validated on the database.

The outline of this paper is the following. In section II, a study of the different features for drowsiness estimation is processed by data-mining. This study leads to the proposition of a drowsiness detection system based on fuzzy logic. Results obtained on a database of 20 different drivers are presented and discussed in section III.

II. DROWSINESS DETECTION SYSTEM

The most relevant blinking features for drowsiness estimation are first selected from a large expertized dataset by data-mining. Fuzzy logic is then used to merge these different features and detect drowsiness.

A. Video based blink detection and characterization

The method to extract blink features from video has been presented in a recent study [START_REF] Picot | Comparison between EOG and high frame rate camera for drowsiness detection[END_REF]. The general principle of this method is shown in Fig. 1.

Twelve features extracted with this method were compared to those extracted from the EOG. It has been shown that the following features: duration at 50% (D50) which is the duration calculated from the half rise amplitude to the half fall amplitude, percentage of eye closure at 80% (P 80) which is the percentage where the eyes were closed at least 80%, the frequency (F ) which is the number of blinks during a defined period, the amplitude-velocity ratio( A/P CV ) which is the ratio between the amplitude of the blink and the maximum speed during the closing period and the energy of eye closure (E C ) which the maximum energy of moving edges during the closing period, can be estimated from a 200Hz video with the same accuracy as from EOG. In the following, a drowsiness detection using these features is presented. The analysis is done using features extracted from EOG but, considering the results obtained in [START_REF] Picot | Comparison between EOG and high frame rate camera for drowsiness detection[END_REF], the detection could be implemented using features extracted from video and would reach the same

The relevancy of these different features to estimate drowsiness is studied here by data-mining to find which features can be used to design the drowsiness detection system.

B. Material

The database used in this study is consistent. It includes 60 hours of driving data from 20 different drivers. Each subject has been recorded while driving on a simulator for 90 minutes, a first time perfectly rested and a second time suffering from sleep deprivation (the subject had slept for 4 hours only). Each recording includes four EEG channels and one EOG channel. Only the EOG channel is used in this paper. The EOG is recorded at 250Hz. Data acquisition was performed by the CEPA (Centre d' Études de Physiologie Appliquée), Strasbourg, FR.

Objective sleepiness was evaluated every 20s on each recording by an expert doctor using the scale described in table I [START_REF] Muzet | AWAKE pilot study no.2: Testing steering grip sensor measures[END_REF]. As only drowsiness visual signs are studied here, the five expert decision levels are converted into a binary decision by considering as "drowsy" any decision superior or equal to 2. Otherwise, the driver is considered as "awake". This represents 6461 samples classified as "awake" and 1096 as "drowsy". Two thirds of this database is randomly selected. It is used to study the relevancy of the selected features. This represents 4384 "awake" samples and 728 "drowsy" samples which is large enough to extract information about the different features by data-mining. The last third (2077 "awake" samples and 368 "drowsy" samples) is used to validate the drowsiness detection system presented in II-D.

C. Feature selection

The purpose of this section is to study the relevancy of each feature mentioned in section II-A with respect to drowsiness. The mean value of each feature is computed every second on a sliding window. Three sliding windows' lengths (20s, 40s and 60s) have been studied. The drowsiness detection is processed on each feature by using a threshold. Two ways of choosing the threshold are discussed here. The first one is to use a fixed common threshold s F for all drivers. The second one is to compute the threshold according to equation (1) for each driver. In this equation, m is the mean value of the studied feature, σ its variance and k a multiplying constant. The values m and σ are learnt at the beginning of the recording when the driver is considered as awake. The threshold is noted s 1min when it is learnt on the first minute and s 5min when it is learnt on the first five minutes. This method has been inspired by the work of Svensson [START_REF] Svensson | Blink behaviour based drowsiness detection method development and validation[END_REF] on EOG.

s = m + k • σ (1)
This method provides a binary decision [awake; drowsy] every second while the database is expertized by time slices of 20s. To compare the results obtained with expert's decision, the following technique is used. Each 20s slice classified by the expert is directly compared to the system decision: if during the 20s interval, the system classifies at least 10s as "drowsy", then the decision for the slice is "drowsy". Else it is "awake".

The results are compared using a Receiver Operating Characteristic (ROC) curves, plotting true positive rate (T P rate ) in function of the false positive rate (F P rate ). The T P rate is calculated as the ratio between the number of "drowsy" samples correctly classified by our system and the number of "drowsy" expertized samples. The F P rate is computed as the ratio between the number of false "drowsy" detection on the number of "awake" expertized samples. The purpose is to have the highest T P rate with the lowest F P rate .

Fig. 2 shows the results for each feature and different thresholds. In this figure, the "square" markers correspond to the results obtained with a fixed threshold s F , the "circle" markers to those obtained with a threshold s 1min learnt during the first minute and the "triangle" markers to those obtained with a threshold s 5min learnt during the first five minutes. These results are obtained by varying the different threshold values in a coherent range for each feature. The results obtained for D50, P80, F and A/PCV only are shown. Indeed, the E C feature provides the worst results (T P rate =69,3% and F P rate =32,3% with an optimal threshold). Moreover, the results have shown that the best results are obtained with the 20s sliding window. So, only the results obtained with a sliding window of 20s are presented here.

It is obvious from Fig. 2 that results are better when a fixed threshold is used for all the features. This means that the first minutes of the recording provides a poor estimation

TABLE II BEST RESULTS AND THRESHOLD OBTAINED FOR EACH FEATURE

of the "awake" state for the blinking features. Results obtained with s 5min are also much better for all the features than those obtained with s 1min which means that the first minute of video recording is not long enough of the driver ocular behaviour to be correctly assessed. Nevertheless, results obtained using s 5min are lower than those obtained using s F . In the following sections, only fixed thresholds will be considered.

The best results as well as the best threshold s F for each feature are presented in table II. This table shows that the best feature to detect drowsiness is P 80 (T P rate = 83, 2% and F P rate = 21, 3%). This is in concordance with the literature [START_REF] Knipling | PERCLOS: a valid psychophysiological measure of alertness as assessed by psychomotor vigilance[END_REF], [START_REF] Galley | Blink parameter as indicators of driver's sleepiness -possibilities and limitations[END_REF]. Moreover, D50, P80, F and A/PCV only will be used to design the drowsiness detection system as they all provides good results.

D. Drowsiness detection using fuzzy logic

The method proposed to detect drowsiness use fuzzy logic to merge the four features considered. The first step of the method is to transform each feature into a fuzzy variable whose value is between 0 and 1. Each fuzzy variable is defined as the membership degree to the "drowsy" state: the closer to 1 the value of the variable, the higher the probability of being drowsy. The membership function is defined by equation ( 2). Coefficients a and b are respectively defined by equations ( 3) and ( 4) depending on the threshold s F , selected in II-C for each feature, in order to standardize the membership function. The values -0.25 and +0.25 have been chosen to centre the membership function on the considered threshold s F .

f a,b (x) =    0 if x ≤ a x-a b-a if a ≤ x ≤ b 1 if x ≥ b (2) a = s F -0, 25 • s F (3) b = s F + 0, 25 • s F (4)
The membership degree of each feature D50, P80, F and A/PCV is respectively noted µ D50 , µ P 80 , µ F and µ A/P CV . Then, the fuzzy variables are merged at each time i using the equation ( 5). The fusion is computed at each instant as the average of the different membership degrees.

µ f usion (i) = 1 4 j µ j (i), j ∈ {D50, P 80, F, A/P CV } (5)
The result is noted µ f usion and its value is between 0 and 1. The relevancy of the detection depends on the value of µ f usion : the closer to 1 µ f usion , the more relevant the detection. The value 1 corresponds to the "drowsy" state and 0 to the "awake" state. To make the decision, the driver is considered as "drowsy" if µ f usion is greater than 0,5. Otherwise, he is considered as "awake". Actually, this function provides a compromise between the different features.

III. RESULTS AND DISCUSSION

The decisions of the drowsiness detection system are made every second. The system decision is directly compared to the expert's one by 20s slices, in the same way as in the section II-C. The method has been tested on the third of the database which has not been used to choose the different thresholds. The different results are shown in Fig. 3. The "square" markers correspond to results obtained with the data fusion method, the "triangle" markers to those obtained with D50 only, the "star" markers to those obtained with P80 only, the "hexagram" markers to those obtained with F only and the "diamond" markers to those obtained with A/PCV only. Table III presents the T P rate and the F P rate related values.

First, results obtained using each feature alone proves that the thresholds used are robust. The results obtained on the last third of the database are indeed very similar to those found during the study of the different features on the two other thirds of the database. So, the method is independent of the driver since the same thresholds give good results on 20 different drivers. Moreover, results obtained with the fusion method are better than those obtained with only one feature. The fuzzy fusion gives T P rate = 81, 7% which is as good as the average T P rate found with the features processed independantly (around 80%). Indeed, good detections on the different features are corresponding which explains that the results of the fusion provides a T P rate equivalent as those found for each feature independently. Furthermore, the F P rate = 13, 4% obtained with the fusion is obviously better than those found for each feature. This means that the false alarms obtained with the different features are not occuring at the same time. Thus, when the different features are merged, the false alarms produced by one of the feature are eased by the others which explains the decrease of the F P rate . So, the fusion of the different features make the detection more robust: it decreases the number of false alarms without changing the number of correct detections. These results are confirmed by a McNemar test comparing the performances obtained by fusion to those obtained with the features alone.

The results presented are also good compared to those found in the literature. Johns [START_REF] Johns | The amplitude-velocity ratio of blinks: a new method for monitoring drowsiness[END_REF] obtained 75% of correct detections and 30% of false alarms while detecting drowsiness with the amplitude-velocity ratio only. This is coherent with the results found for the A/PCV feature (T P rate = 78.6% and F P rate = 23.9%). A fuzzy expert system merging blink duration, blink frequency and blink interval to estimate drowsiness has been propopsed by Damousis et al. [START_REF] Damousis | A fuzzy expert system for the early warning of accidents due to driver hypo-vigilance[END_REF]. They obtained 90% of correct detections with a fuzzy expert system which is slightly better than our results. Nevertheless they also had 30% of false alarms which is worse than the F P rate obtained by the proposed method. This means that the features used in these paper for the drowsiness detection are more relevant. However, it is necessary to keep some reserves since these results have not been obtained on the same database. Nevertheless, the database used in this paper is large enough (60 hours of driving from 20 different drivers) to make a good assessment of the performances of our system.

IV. CONCLUSIONS AND FUTURE WORKS

An universal algorithm for the automatic detection of drivers drowsiness based on visual signs has been presented here. This algorithm is based on the fuzzy fusion of several blinking features that can be extracted from a high frame rate video. The different features used in the system has been selected by data-mining on a consistent database. The algorithm has been tested on a large number of drivers and does not need to be tuned. Its validation on EOG signals gives good results with 81,7% of correct detections and only 13,4% of false alarms. These results should be equivalent when using a high frame rate camera since the features used here have been cautionly chosen and can be extracted from EOG and high frame rate video with the same accuracy.

The next step of this work is the achievement of an acquisition campaign of drowsy data with a high frame rate camera so as to validate the whole method proposed on video data. Once validated, this video based drowsiness detection system should be added to a EEG based drowsiness detection system ( [START_REF] Picot | Recent Advances in Biomedical Engineering, chapter Monitoring drowsiness on-line using a single encephalographic channel[END_REF]) to merge the decisions to obtain a highly reliable automatic drowsiness detector.

Fig. 1 .

 1 Fig. 1. Blink detection method

Fig. 2 .

 2 Fig. 2. Results obtained with the different features Feature Results Threshold D50 T Prate = 81.4% s F = 130ms F Prate = 24.9% P80 T Prate = 84, 2% s F = 0.75% F Prate = 22, 3% F T Prate = 78.3% s F = 0, 1Hz F Prate = 18, 6% A/PCV T Prate = 78.5% s F = 17.25 F Prate = 25.3%

Fig. 3 .

 3 Fig. 3. Results of the drowsiness detection Fuzzy D50P80 F A/PCV logic T Prate 81.7% 80.7% 82.8% 77.0% 78.6% F Prate 13.4% 25.1% 21.4% 19.3% 23.9%

TABLE III RESULTS

 III OF THE DROWSINESS DETECTION

ACKNOWLEDGEMENT

The authors are grateful to the Centre d' Études de Physiologie Appliquée (CEPA) in Strasbourg (FR) for providing the data and their help.


Alain Y Kibangou

Cyrille Siclet

Laurent Ros

ZF OFDM Receiver for Underwater Communications

Keywords: Underwater acoustic communications, OFDM, Doppler distorsion, Hankel matrix, Vandermonde decomposition, Total least squares

In this paper, we propose a new scheme of transmission and reception of OFDM ( Orthogonal Frequency Division Multiplexing) signals for underwater wireless communications. The transmitter makes use of a kind of superimposed pilot signals, which are used by the receiver for estimating both the channel parameters (paths gain and delay) and the Doppler. Then, the estimated parameters are used for estimating the informative symbols in the zero forcing (ZF) sense. A new estimation scheme based on the resolution of two harmonic retrieval problems is proposed. Unlike standard methods, data resampling and estimation of residual carrier frequency offset are avoided. The efficiency of the proposed scheme is evaluated by means of simulations.

I. INTRODUCTION

Nowadays, a great interest is dedicated to study fleets of autonomous underwater vehicles (AUV) working together to reach a common objective such as a gradient search and following for source detection. For this purpose derivation of distributed, multi-vehicle co-ordination schemes and development of efficient underwater acoustic communication protocols are needed. For achieving high data rate and large system capacity, Orthogonal Frequency Division Multiplexing (OFDM) has been claimed to be an efficient communication technology. It allows designing low complexity receivers to deal with highly dispersive channels. This facts motivates the use of OFDM in underwater environments.

Underwater acoustic channels are wideband in nature due to the small ratio of carrier frequency to the signal bandwidth, which introduces frequency-dependent Doppler shifts [START_REF] Li | Multicarrier communication over underwater acoustic channels with nonuniform doppler shifts[END_REF]. They also exhibit several propagation paths.

In order to adequately recover the transmitted information, algorithms at the receiver must include estimation and compensation of the Doppler scaling factor, channel estimation, and information symbols estimation. Several approaches have been suggested in the literature for estimating the Doppler scaling factor. They are based on the use of preamble and postamble of a packet consisting of multiple OFDM blocks Alain Y. Kibangou is with GIPSA-Lab, System control Department, University Joseph Fourier, CNRS, 961 rue de la Houille Blanche, B.P. 46-38402 Saint Martin d'Hères, France (e-mail: Alain.Kibangou@gipsa-lab.grenobleinp.fr) Cyrille Siclet and Laurent Ros are with GIPSA-Lab, Image and Signal Department, University Joseph Fourier, CNRS, 961 rue de la Houille Blanche, B.P. 46-38402 Saint Martin d'Hères, France (e-mail: Cyrille.Sicletu@gipsalab.grenoble-inp.fr, Laurent.Ros@gipsa-lab.grenoble-inp.fr) This work was supported by the ANR and the EU funded projects CONNECT and FeedNetBack [START_REF] Li | Multicarrier communication over underwater acoustic channels with nonuniform doppler shifts[END_REF] or by exploiting correlation induced by the cyclic prefix [START_REF] Kim | Parameter study of OFDM underwater communications system[END_REF]. Then, the received signal is resampled by using a sampling period related to the estimated Doppler scaling factor. It is also necessary to estimate and compensate for the residual carrier frequency offset (CFO) since the Doppler can vary between consecutive OFDM blocks inside a given packet.

In this paper, the received data are processed block-byblock. We make use of high resolution harmonic estimation methods to estimate both Doppler scaling factor and channel parameters (path gains and delays). The advantage of the proposed scheme is to avoid data resampling and residual CFO estimation and compensation. The estimated channel parameters and Doppler are then used to build a channel matrix that is used for estimating the informative symbols following a zero forcing scheme. Notation: We denote by V (x, N) the L × N Vandermonde matrix with x ∈ C L as generator. Its first column is constituted by 1s.

II. SIGNAL MODEL

The proposed receiver is based on block-by-block processing. We assume that consecutive OFDM symbols are separated by a sufficient guard interval in order to avoid intersymbol interference. Moreover, the transmitter and the receiver clocks are synchronized.

In the transmission scheme considered in this paper, the signal s(.), results on the superposition of two OFDM signals s 0 (.) and s 1 (.), i.e. s(t) = s 0 (t) + s 1 (t). The information contained in s 0 (.) is known to the receiver and will serve to estimate the channel parameter, whereas s 1 (.) contains the informative symbols unknown to the receiver. The transmitted bandpass signal is then given by:

s(t) = 1 ∑ i=0 s i (t), with s i (t) = Re ∑ k∈K i c k,i e j2π f k,i t Π T i (t -τi )
where:

• K i = {-K i /2, • • • , K i /2 -1}, K i
being the number of subcarriers used for transmitting the learning sequence

{c k,0 }, if i = 0, or the informative sequence {c k,1 }, if i = 1. • f k,i = f c + k∆ f i ,
are the subcarrier frequencies, ∆ f i being the frequency spacing between two consecutive subcarriers and f c is the (central) carrier frequency.

• Π T i (t) = 1 if t ∈ [0, T i ]
and 0 otherwise. • τ0 = 0 whereas τ1 , which will be defined later, stands for the delaying operation of the informative signal s 1 (.) with respect to the learning signal s 0 (.).

The signal bandwidth equals

B = K 1 ∆ f 1 = K 1 /T 1 .
We assume that less subcarriers are used for transmitting the pilot sequence than the informative one, i.e. K 0 < K 1 . Moreover, we have ∆ f 0 = B/K 0 , meaning that the frequency spacing in s 0 is more important than in s 1 (.).

The signal s(.) is transmitted to the receiver through an underwater acoustic channel modeled as:

h(t, τ) = P ∑ p=1 A p (t)δ (t -τ p (t)),
where A p (t) and τ p (t) are respectively the real-valued gain and the delay associated with the pth propagation path. In addition, the following assumptions are adopted:

• All paths are affected by a similar Doppler scaling factor a such that τ p (t) = τ p -at.

• The path delays τ p , the gains A p , and the Doppler scaling factor a are constant over the block duration T = T 0 +T 1 + T g -τ1 .

• The maximal path delay τ p is known and the guard interval T g is sufficiently large to avoid interference between consecutive blocks. We have

T g > max{τ p } = τ max .
The analytic representation of the received signal y(t) is then given by:

y(t) = 1 ∑ i=0 ∑ k i ∈K i P ∑ p=1 c k i ,i A p Π T i ((1 + -τ p -τi ) ×e -j2π f k i ,i τ p e j2π f k i ,i (1+a)t . ( 1
)
Assuming that the bounds of the path delays, τ min 1 and τ max , and that of the Doppler scaling factor are known 2 . We get:

τ min ≤ τ p ≤ τ max , λ min ≤ 1 + a ≤ λ max .
Therefore, for t ∈ ( T 0 +τ max λ min , T ) and τ1 ≥ λ min T -τ max -T 1 , we get the signal part conveying the informative data:

y 1 (t) = ∑ k 1 ∈K 1 h k 1 (t)c k 1 ,1 (2)
with

h k 1 (t) = P ∑ p=1 A p Π T 1 ((1 + a)t -τ p -τ1 )e -j2π f k 1 ,1 τ p e j2π f k 1 ,1 (1+a)t
(3) By sampling the received signal at the period T e ≥ 1 B , equation (2) can be written in matrix form as follows:

y = Hc (4)
where y, H, and c contain respectively the samples of received signal corresponding to t ∈ ( T 0 +τ max λ min , T ), the channel matrix, whose entries are defined in (3), and the informative symbols. 1 The knowledge of the minimal path delay is difficult. Therefore, we get the absolute minimum τ min = 0.

2 λ max and λ min are related to the maximal and minimal velocities of the underwater vehicles, which can be a priori known.

Therefore, if the channel parameters and the Doppler factor are known, the ZF receiver is obtained by solving

ĉ = H H H -1 H H y. (5)
In the following section, we derive a joint channel and Doppler estimation method.

III. DOPPLER AND CHANNEL ESTIMATION

Let us select t 0 and t 1 such that

τ max λ min ≤ t 0 < t 1 ≤ T 0 + τ min λ max . ( 6
)
For t ∈ [t 0 ,t 1 ], the received signal can be written as:

y(t) = ∑ k 0 ∈K 0 c k 0 ,0 B k 0 e j2π f k 0 ,0 (1+a)t (7)
with

B k 0 = P ∑ p=1 A p e -j2π f k 0 ,0 τ p . ( 8
)
In discrete-time, we get:

y n = ∑ k 0 ∈K 0 c k 0 ,0 B k 0 e j2π f k 0 ,0 (1+a)nT e . ( 9
)
Inside the observation window defined by the bounds ( 6), the received signal can be viewed as a mixture of K 0 complex valued harmonics with constant magnitudes. As a consequence, the magnitudes c k 0 ,0 B k 0 and the phases 2π f k 0 ,0 (1 + a)T e can be obtained by solving a harmonic retrieval problem. From [START_REF] Papy | Exponential data fitting using multilinear algebra: The single-channel and multi-channel case[END_REF], one can note that the magnitude B k 0 can also be viewed as a mixture of P harmonics. So, for estimating both the Doppler scaling factor, the channel gains and delays, we have to solve two harmonic retrieval problems. In the literature, the harmonic retrieval problem have been extensively studied [START_REF] Kay | Modern spectral estimation, Theory and Application[END_REF]- [START_REF] Papy | Exponential data fitting using multilinear algebra: The single-channel and multi-channel case[END_REF]. In this paper, we derive new estimators by making use of HTLS (Hankel Total Least Squares) [START_REF] Papy | Exponential data fitting using multilinear algebra: The single-channel and multi-channel case[END_REF], [START_REF] Van Huffel | Enhanced resolution based on minimum variance estimation and exponential data modeling[END_REF], a high resolution method. In the sequel, we first give the main guidelines for solving the harmonic retrieval problem before deriving our estimators.

A. Harmonic retrieval using HTLS [START_REF] Papy | Exponential data fitting using multilinear algebra: The single-channel and multi-channel case[END_REF] Let us consider the samples z q such that z

q = K ∑ k=1 α k e jϕ k q ,q = q 0 , • • • , q 0 + N -1. By setting N = L + M -1, L > Q, M ≥ Q,
we first build the Hankel matrix Z ∈ C L×M having respectively as first column and last row the vectors

Z .1 = z q 0 z q 0 +1 • • • z q 0 +L-1 T and Z L. = z q 0 +L-1 z q 0 +L • • • z q 0 +N-1 .
It admits the two following decompositions (10) and (11):

Z = S diag(α α α)T T . ( 10
) with α α α = (α 1 , • • • , α K ) T , S T = V (φ φ φ , L), T T = V (φ φ φ , M), φ φ φ = e jϕ 1 , • • • , e jϕ K T
, and

Z = UΣ Σ ΣV H , ( 11
)
where U ∈ C L×K , V ∈ C M×K , and Σ Σ Σ ∈ C K×K contain respectively the left and right singular vectors and the singular values of Z. From (10) and (11), we deduce that U and S generate the same subspace. Hence, it exists a nonsingular matrix Q ∈ C K×K such that U = SQ. The Vandermonde matrix S possesses the shift-invariance property expressed as: S diag(φ φ φ ) = S, where the line on the bottom (resp. on the top) of a matrix stands for deleting the first (resp. the last) row. We also have:

U = SQ, Ū = SQ
By combining these two above equations, we get:

Ū = UQ -1 diag(φ φ φ )Q. ( 12
)
Hence, the poles e jϕ k are the eigenvalues of Φ Φ Φ = Q -1 diag(φ φ φ )Q. Therefore, we have to estimate Φ Φ Φ and then deducing its eigenvalues. Φ Φ Φ is estimated by solving Ū = UΦ Φ Φ in the total least squares sense. We get:

Φ Φ Φ = -W 12 W -1 22 ( 13
)
where

W = W 11 W 12 W 21 W 22 ∈ C 2K×2K
is the matrix containing the right singular vectors of U Ū . Once the poles estimated, the magnitude α k can be obtained as the least squares solution of the vectorized Vandermonde decomposition (10):

vec(Z) = (T S)α α α, ( 14
)
where denotes the Khatri-Rao product, the poles in S and T being replaced by their estimated values.

B. Doppler scaling factor estimation

Let us consider the N samples y n , n = n 0 , • • • , N + n 0 -1, of the received signal y n as given in [START_REF] Van Huffel | Enhanced resolution based on minimum variance estimation and exponential data modeling[END_REF]. Herein, the pole and the magnitude of the involved harmonics are φ k 0 = e j2π f k 0 ,0 (1+a)T e and α k 0 = c k 0 ,0 B k 0 e j2π f k 0 ,0 (1+a)n 0 T e . Let us denote by φk 0 the estimated poles using the HTLS algorithm. If the sampling period T e is chosen such that T e ≤ 1 2λ max max{ f k 0 ,0 } then the angle of φk 0 , denoted ∠ φk 0 , belongs to [-π, π]. In fact, the estimated eigenvalues equal the actual poles up to a permutation. For removing such an ambiguity, the estimated eigenvalues can be sorted with their pulsation in an ascending order. Therefore, we deduce the following estimator for the Doppler scaling factor:

â = -1 + 1 K 0 ∑ k 0 ∈K 0 ∠ φk 0 2π f k 0 ,0 T e (15)

C. Path delays and gains estimation

Let us denote αk 0 the magnitudes obtained by solving (14), where the poles in S and T have been replaced by their estimated values. We get:

Bk 0 = c * k 0 ,0 c k 0 ,0 2 αk 0 e -j2π f k 0 ,0 (1+ â)n 0 T e , k 0 ∈ K 0 . ( 16
)
Note that we can rewrite B k 0 , defined in (8), as

B k 0 = P ∑ p=1 A p e -j2π f -K 0 /2,0 τ p e -j2πk 0 ∆ f 0 τ p , k 0 ∈ K 0 .
Provided, K 0 ≥ 2P, we can solve this harmonic retrieval problem using the HTLS method, the poles and the magnitudes being e -j2π∆ f 0 τ p and A p e -j2π f -K 0 /2,0 τ p respectively.

Let us denote βp the estimated harmonic magnitudes using HTLS. We can then estimate the gains A p as the magnitudes

Âp = βp . ( 17
)
We also deduce

θ p = βp | βp| = e -j2π f -K 0 /2,0 τ p .
The estimation of the path delays is a more complicate task. Indeed, we cannot guarantee that 2π∆ f 0 τ p belongs to [-π, π]. As a consequence, we cannot directly deduce the value of the path delays from the angle of the entries ψp of β β β . However, it exists an integer m 0 such that

∠ ψp = -2π∆ f 0 τ p + 2πm 0 .
For different integers m, we define a set of possible path delays τ p,m such that τ p,m < T g :

τ p,m = - ∠ ψp 2π∆ f 0 + m ∆ f 0 and θ p,m = e -j2π f -K 0 /2 τ p,m .
Let us now consider the square distance

κ p,m = θ p -θ p,m 2 = 1 -e -j2π f -K 0 /2 m -m ∆ f 0 2 (18)
After few manipulations, it can be shown that this square distance is minimal only for τ p,m = τ p , the actual path delay:

τp = arg min κ p,m . (19)

IV. SIMULATION RESULTS

In these simulations, the range of frequency used by the underwater vehicles was [15kHz -28kHz]. We considered Doppler scaling factor lower than 10 -3 , meaning that the maximal relative speed was 1.5 m/s. The carrier frequency was set equal to f c = 21 kHz, whereas the guard interval was T g = 10 ms. We used K 0 = 16 sub-carriers for the learning sequence and K 1 = 512 for the informative one. The duration of the learning sequence was T 0 = 2T g whereas that of the informative sequence was T 1 = 46.54 ms. The sampling frequency was 1/T e = 56.056 kHz. We made use of the QPSK modulation. All the results presented below are averaged values over 100 independent Monte Carlo runs. The additive noise was a complex valued white Gaussian noise. The paths gain and delay of the simulated channel are depicted in figure 1. We evaluate the performance of the derived ZF receiver in terms of BER (Bit Error Rate). Figure 2 and3 depict the mean of the BER for channels with P = 4 and P = 6 for two different values of the Doppler scaling factor. We compare the performance of the equalization method when the Doppler and the channel parameters are estimated by the method proposed herein and with the ideal case, i.e. the actual parameters are used for constructing the channel matrix. In the ideal case, by increasing the SNR, we get better performance for a channel with P = 6 than that with P = 4. Since, the first channel has more propagation paths, the received signal is more redundant and therefore allows improving the informative sequence estimation. When using the estimated parameters, we don't have the same behaviour. Indeed, with the same amount of data, the parameter estimation is degraded by increasing the number of propagation paths. Moreover, the gap in performance between the ideal case and the estimated one increases with the number of propagation paths.

V. CONCLUSION

In this paper, we have presented a new zero-forcing OFDM receiver for underwater acoustic communications. In the proposed scheme, each OFDM symbol is constituted by a the superposition of two OFDM sub-blocks characterized by different number of subcarriers. The first OFDM sub-block is used for estimating the channel parameters and the Doppler scaling factor. The estimation methods are based on high resolution methods for solving harmonic retrieval problems. Then, the estimated parameters are used for retrieving the informative symbols according to the Zero forcing criterion. The proposed scheme is particularly efficient for moderate to high SNR. Robustness to noise, for lower levels of SNR, should be investigated. Extension to multiusers underwater communications is under investigation.

Figure 1 .

 1 Figure 1. Paths gain and delay of the simulated channel.

Figure 2 .

 2 Figure2. BER obtained with the proposed estimation method (a = 6.67 × 10 -4 ).

Figure 3 .

 3 Figure 3. BER obtained with the proposed estimation method (a = 3.33 × 10 -4 ).