Zahia Zidelmala

Ahmed Amirou
email: ahmed.amirou@uha.fr

Djaffar Ould Abdeslam
email: djaffar.ould-abdeslam@uha.fr

Jean Merckle
email: jean.merckle@uha.fr

Z Zidelmal
email: z-zidelmal@mail.ummto.dz

D Ould-Abdeslam

ECG beat classifcation

Keywords: ECG beat classification, Support Vector Machines (SVMs), classification cost

come

Introduction

The ElectroCardioGram (ECG) is a bio-electric signal that records the electrical activities of the heart. The ECG waveform is characterized by the waves: P, QRS and T. The most important wave is the QRS complex characterizing the ventricular contractions. Distances between them (RR intervals) define the rhythm, which is influenced by the emotions and physical activity.

The normal electrical conduction in the heart allows the impulse generated by the Sino-Atrial (SA) node to be rapidly propagated to the myocardium. This stimulated myocardium allows efficient contraction of the heart. In heart conduction disorders, ventricular excitation may not originate in the SA node but from other ectopic centers in the myocardium leading to a Premature Ventricular Contraction (PVC) also called extra-systole or ectopic beat. Typically, PVCs are recognized by irregular RR intervals. The normally generated contractions produce normally shaped QRS complexes, while ectopic heartbeats generate a variety of QRS waveforms, quite differing from the normal ones.

Counting the occurrence of ectopic beats is of particular interest to support the detection of ventricular tachycardia and to evaluate the regularity of the depolarization of the ventricles. For example, the risk of Ventricular Fibrillation (VF) for patients with a structural heart disease is higher with an increased occurrence of PVCs. The automatic detection and classification of the ventricular contractions as normal or premature is a subject of long-term studies. This is the basis of the rhythm analysis, which is usually applied to continuous 24-hour ECG recordings (Holter systems) to identify rhythm disorders.

Several algorithms with signal processing techniques have been proposed for detection and classification of heartbeats. Classical techniques extract heuristic ECG descriptors, such as QRS morphology [START_REF] Ince | A generic and Robust System for Automated Patient-Specific Classification of ECG signals[END_REF][START_REF] De-Chazal | A Patient Adapting Heartbeat Classifier Using ECG Morphology and Heartbeat Interval Features[END_REF][START_REF] Christov | Comparative study of morphological and time-frequency ECG descriptors for heartbeat classification[END_REF] and interbeat R-R intervals [START_REF] Ince | A generic and Robust System for Automated Patient-Specific Classification of ECG signals[END_REF][START_REF] De-Chazal | A Patient Adapting Heartbeat Classifier Using ECG Morphology and Heartbeat Interval Features[END_REF][START_REF] Khazaee | Classification of electrocardiogram signals with support vector machines and genetic algorithms using power spectral features[END_REF][START_REF] Yu | Electrocardiogram beat classification based on wavelet transformation and probabilistic neural network[END_REF]. Other ECG descriptors rely on QRS frequency components [START_REF] Christov | Comparative study of morphological and time-frequency ECG descriptors for heartbeat classification[END_REF][START_REF] Khazaee | Classification of electrocardiogram signals with support vector machines and genetic algorithms using power spectral features[END_REF][START_REF] Talbi | PVC discrimination using the QRS power spectrum and self-organizing maps[END_REF][START_REF] Yu | Electrocardiogram beat classification based on wavelet transformation and probabilistic neural network[END_REF] or matching pursuits for extraction of time-frequency beat descriptor [START_REF] Christov | Comparative study of morphological and time-frequency ECG descriptors for heartbeat classification[END_REF]. Some methods apply QRS template matching procedures [START_REF] Krasteva | QRS Template Matching for Recognition of Ventricular Ectopic Beats[END_REF], wavelet and principal component analysis [START_REF] Ince | A generic and Robust System for Automated Patient-Specific Classification of ECG signals[END_REF][START_REF] Yu | Electrocardiogram beat classification based on wavelet transformation and probabilistic neural network[END_REF][START_REF] Shen | Detection of cardiac arrhythmia in electrocardiograms using adaptive feature extraction and Modified Support Vector Machines[END_REF] and higher order statistics [START_REF] Kutlu | Feature extraction for ECG heartbeats using higher order statistics of WPD coefficients[END_REF][START_REF] Melgani | Classification of electrocardiogram signals with support vector machines and particle swarm optimization[END_REF].

Several discriminative techniques have been developed for ECG beat classification such as Artificial Neural Networks [START_REF] Ince | A generic and Robust System for Automated Patient-Specific Classification of ECG signals[END_REF][START_REF] De-Chazal | A Patient Adapting Heartbeat Classifier Using ECG Morphology and Heartbeat Interval Features[END_REF] or probabilistic neural networks [START_REF] Yu | Electrocardiogram beat classification based on wavelet transformation and probabilistic neural network[END_REF]. Other work used the K th nearest-neighbor rule [START_REF] Christov | Comparative study of morphological and time-frequency ECG descriptors for heartbeat classification[END_REF][START_REF] Kutlu | Feature extraction for ECG heartbeats using higher order statistics of WPD coefficients[END_REF] and genetic algorithms [START_REF] Khazaee | Classification of electrocardiogram signals with support vector machines and genetic algorithms using power spectral features[END_REF]. Among all these methods, Support Vector Machines (SVMs) known as excellent tool for classification and regression problems has shown success in this application field e.g. [START_REF] Khazaee | Classification of electrocardiogram signals with support vector machines and genetic algorithms using power spectral features[END_REF][START_REF] Shen | Detection of cardiac arrhythmia in electrocardiograms using adaptive feature extraction and Modified Support Vector Machines[END_REF][START_REF] Melgani | Classification of electrocardiogram signals with support vector machines and particle swarm optimization[END_REF].

Even though the performance of all these techniques, misclassifications cannot be completely eliminated and, thus, can produce severe penalties. This motivates us to introduce a cost sensitive classifier to heartbeat recognition. This option has a particular interest in medical field where the consequences of misclassification can be substantial. In such situations, rejection may even be more preferable than risking misclassification.

This paper is organized as follows. Section 2 presents data signal and preprocessing techniques. Section 3 provides a detailed description of the feature extraction methodology. The classifier model proposed in this study is developed in section 4. It recalls Bayes rule for binary classification with weighted errors and presents the SVMs in this framework and the learning criterion dedicated for the problem at hand. This proposed method is tested empirically in Section 5. Finally, Section 6 briefly concludes the paper.

Signal processing

ECG data

The annotated ECG records from the MIT-BIH arrhythmia data base [START_REF] Mark | MIT-BIH arrhythmia database[END_REF] has been used in this study. This data base has 48 records. Each record is of length 30 min with 360 Hz sampling frequency. The data were recorded in two channels (modified limb lead II and modified limb lead VI) of surface ECG from long term Holter recorders. They present a variety of waveforms, artifacts, complex ventricular, junctional and supraventricular arrhythmias and conduction abnormalities. Each record is accompanied by an annotation file in which each ECG beat has been identified by expert cardiologists. These labels referred to as 'truth' annotation and are used to develop the classifier and to evaluate its performance in the testing phase.

Since this study is to evaluate the performance of a binary classifier with a reject option, we followed the American Heart Association (AHA) records equivalent annotation [START_REF] Mark | MIT-BIH arrhythmia database[END_REF] and the AAMI recommended practice [START_REF] Mark | AAMI-recommended practice: Testing and reporting performance results of ventricular arrhythmia detection algorithms[END_REF] to form two classes of heartbeats: (i) the positive class representing the ventricular ectopic beats (V) including premature ventricular contractions, fusion of ventricular and normal beats and unknown beats; (ii) the negative class representing the normal beats (N), including all normal heartbeats (approximately 70% of the database) and some of the abnormal beats (Left or Right Bundle Branch Block, aberrated atrial premature beat, nodal or atrial premature beat, nodal escape beat, atrial ectopic beat and nodal ectopic beat). We further restricted the N class to contain only heartbeats, which are representative for the predominant rhythm of the patient: Normal beats, Left Bundle Branch Block beats (L) and Right Bundle Branch Block beats (R). In agreement with the AAMI recommended practice, records containing paced beats (102, 104, 107, and 217) were excluded. From this study, records with no PVC beats (11 records) were also excluded leaving 33 records of interest. The forms of predominant heartbeats are displayed on figure 1.

ECG Filtering

The objective of this paper is to classify heartbeats. Before performing this task, several pre-processing steps were performed on the raw data. In fact the electrocardiogram (ECG) from body electrodes are corrupted by noise. Usually, two principal sources of ECG noise can be distinguished: the first one caused by the physical parameters of the recording equipment and the second one representing the bioelectrical activity also called background activity or baseline wander.

Several noise removal techniques were recently developed. ECG is a nonstationary biosignal corrupted by additive noise. An efficient tool for denoising such signal is the wavelet transform. The latter is a time-scale representation which describes a signal by using the correlation with translation and dilatation of a function called mother wavelet. Basic theory can be found in many papers, i.e. [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]. The objective of wavelet based denoising process is to estimate the signal of interest s(i) from the composite one f (i) by discarding the noise e(i) assuming that f (i) = s(i) + e(i) .

The noisy signal is decomposed into six levels by Discrete Wavelet Transform (DWT) using the Daubechies wavelet db4. Passing the obtained coefficients through a threshold δ, certain values are set to zero. The denoised signal is recovered by taking the Inverse Discret Wavelet Transform (IDWT) of the resulting coefficients. For each decomposition level j, the threshold δ j is computed according to the noise level. δ j = σ j 2log(N ) where N is the number of samples and σ j is the noise standard deviation which can be estimated from the median of its detail coefficients σ j = M AD(d j i )/0.6745 where MAD is the Median Absolute Deviation of the corresponding level. More details can be found in [START_REF] Johnstone | Wavelet threshold estimators for data with correlated noise[END_REF].

The baseline drift is removed by zeroing the scaling coefficients of the DWT at level 6. This is equivalent to filtering the ECG signal with a H-Pass filter having a cutoff frequency of 2.8Hz. The ECG signal sampled at f s = 360Hz is used as the initial scaling coefficients of DWT. They lie in the bandwidth [0, f s /2]. The decomposition with the analysis filter bank is repeated over the scaling coefficients. Each new iteration of the filter bank divides the bandwidth of the scaling coefficients by two. After 6 iterations the resulting scaling coefficients correspond to the low frequency baseline drift of the ECG signal. Now, if these coefficients are set to zero, the reconstructed ECG will have all the details of the original ECG, but with low frequency baseline drift removed. Level 6 is also fixed referring to [START_REF] Zidelmal | QRS detection using wavelet coefficients[END_REF][START_REF] Thakor | Estimation of the QRS Complex Power Spectra for Design of a QRS Filter[END_REF][START_REF] Pan | A real-time QRS detection algorithm[END_REF]. In these papers, it was shown that most of the QRS power lies in [START_REF] Talbi | PVC discrimination using the QRS power spectrum and self-organizing maps[END_REF][START_REF] Yu | Electrocardiogram beat classification based on wavelet transformation and probabilistic neural network[END_REF][START_REF] Krasteva | QRS Template Matching for Recognition of Ventricular Ectopic Beats[END_REF][START_REF] Shen | Detection of cardiac arrhythmia in electrocardiograms using adaptive feature extraction and Modified Support Vector Machines[END_REF][START_REF] Kutlu | Feature extraction for ECG heartbeats using higher order statistics of WPD coefficients[END_REF][START_REF] Melgani | Classification of electrocardiogram signals with support vector machines and particle swarm optimization[END_REF][START_REF] Mark | MIT-BIH arrhythmia database[END_REF][START_REF] Mark | AAMI-recommended practice: Testing and reporting performance results of ventricular arrhythmia detection algorithms[END_REF][START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF][START_REF] Johnstone | Wavelet threshold estimators for data with correlated noise[END_REF][START_REF] Zidelmal | QRS detection using wavelet coefficients[END_REF][START_REF] Thakor | Estimation of the QRS Complex Power Spectra for Design of a QRS Filter[END_REF][START_REF] Pan | A real-time QRS detection algorithm[END_REF][START_REF] Vapnik | The Nature of Statistical Learning Theory[END_REF][START_REF] Grandvalet | Support Vector Machines With a Reject Option[END_REF][START_REF] Wegkamp | Classification methods with reject option based on convex risk minimization[END_REF][START_REF] Hearth | Recommended Practice for Testing and Reporting Performance Results of Ventricular ArrhythmiaDetection Algorithms[END_REF][22] therefore, this baseline removal do not affect the region of interest (QRS window).

ECG beat segmentation

In this study, the R peaks were detected using a robust method based on wavelet coefficients that was introduced in [START_REF] Zidelmal | QRS detection using wavelet coefficients[END_REF]. The peaks Q and S are detected using simple peak detection method. The peak Q is assumed to lie within 50ms before R peak, while S peak is assumed to lie within 100ms after R peak. An example of QRS detection and QS segmentation is displayed on figure 2.

Feature extraction

The QRS complex in ECG signal varies with origination and conduction path of the activation pulse in the heart. When the activation pulse originates in the SA node and travels through the normal conduction path, the normal QRS complex has a sharp and narrow deflection and the spectrum contains high frequency components. When the activation pulse originates in the ventricle and does not travel through the normal path, the QRS becomes wide and the high frequency components of the spectrum are attenuated. A set of algorithms from signal conditioning to measurements of average wave amplitudes, durations, morphology, and areas is usually adopted to perform a quantitative description of a heartbeat. In this study, we select features that clearly characterize and discriminate the quality of the beats such as the R-R intervals, the QRS morphology and the frequency components.

RR intervals

RR intervals provide useful information for clinical diagnosis and identification of pro-arrhythmic events associated with heart-rate variation [START_REF] De-Chazal | A Patient Adapting Heartbeat Classifier Using ECG Morphology and Heartbeat Interval Features[END_REF]. We studied the specific deviations of the RR intervals surrounding the heartbeat in order to give indications about the type of the heartbeat. For a beat i, the RR interval is RR i = R i -R i-1 . This value is normalized to make it independent of the heart rate. The normalized RR interval is

RR norm (i) = RR i -RR max(RR)
where RR and max(RR) are respectively the mean and the maximum of the RR intervals in the same recording. We can see on figure 3 the normalized RR intervals obtained with N, L, R and V beats. We observe a unique peak around zero for all QRS complexes with supraventricular origin (N, R, L) considered in the same class according to the AAMI recommended practice [START_REF] Mark | AAMI-recommended practice: Testing and reporting performance results of ventricular arrhythmia detection algorithms[END_REF] in opposite to V beats leading to a high RR variance. As in [START_REF] De-Chazal | A Patient Adapting Heartbeat Classifier Using ECG Morphology and Heartbeat Interval Features[END_REF], we extracted three features for each beat from the RR sequence: the pre-RR interval is the RR interval between a given beat and the previous one, the post-RR interval is the interval between a given beat and the following one and a local average interval obtained by averaging the ten RR intervals surrounding the considered beat.

Morphological descriptors

• QRS duration varies with origination and conduction path of the activation pulse in the heart so, it is a fundamental feature used for classifying beats. In this work, the QRS duration is represented by the time interval between the two peaks Q and S.

• The morphology of the beat is captured by a four Linear Predictive Coding (LPC) coefficients. The basic idea of this technique is that future values of a discrete signal are estimated as a linear function of previous samples. The most common representation is

y n = p k=1 a k y n-k (1)
where a k is the k th linear prediction coefficient, p is the order of the predictor, y n the present predicted sample and y n-k , the k th precedent sample. Note that a 1 is usually equal to 1.

Frequency features

In figure 4 (a), the estimated Power Spectrum Density (PSD) of windowed (N, R, L and V) beats is displayed. Before computing Fast Fourier Transform (FFT), a 180 ms Blackman window is applied to suppress the discontinuities due to possible adjacent P and T waves. We can see on the same figure (b) that the spectra of V beats distinguish from the spectra of N and R beats at all energy levels because of shifting towards low frequencies. Some confusion can be made between L and V spectra but it is dismissed by RR interval feature which separates perfectly the two types of beats. Accordingly, we extract 6 features, the normalized power levels corresponding to 7.5Hz, 10Hz, 12.5Hz, 15Hz, 17.5Hz and 20Hz to discriminate the beats in the frequency domain.

QS power

To complete the frequency information, we consider detail coefficients at level 4 and 5 obtained with discrete wavelet decomposition. At these levels, the power of each QS segment is computed and used as a classification feature. In [START_REF] Zidelmal | QRS detection using wavelet coefficients[END_REF], the authors have clearly showed that d 4 carries the dominant power of normal QRS complexes while d 5 carries the dominant power of ectopic QRS complexes (see figure 5). In each subband, signal variance represents the averaged power in that band as:

σ 2 x = 1 N N n=1 [x(n) -x] 2 (2)
where x is the sample mean of the signal and N is the number of samples in the considered segment. The value of N differs from one beat to another depending on whether it is large (ectopic beat) or narrow (normal beat).

Data normalization

Using all the observations above, each beat ′ i ′ is represented by a 15 element vector x i . The quantities of the features may be quite different, so their normalization is necessary to standardize all the features to the same level before training the SVM classifier. The formula of the normalization is defined as:

x

′ ij = tansig( x ij -x j σ x j )
From figure 6, one can see that rejecting a pattern is a viable option if and only if r < c 2 . Bays rule can than be stated as:

(x) =      +1 if P (Y = 1|X = x) > P o , -1 if P (Y = 1|X = x) < (1 -P o ) , 0 otherwise . (4)
To minimize the empirical counterpart of the risk 3 computationally not feasible, one could replace it by surrogate loss functions. The most popular are the hinge loss motivated by [START_REF] Vapnik | The Nature of Statistical Learning Theory[END_REF] leading to sparse solutions [START_REF] Grandvalet | Support Vector Machines With a Reject Option[END_REF][START_REF] Wegkamp | Classification methods with reject option based on convex risk minimization[END_REF] and the logistic regression model offering ability to estimate the posterior probability P (Y = 1|X = x) and then a good choice of the threshold δ. In this study, P (Y = 1|X = x) have to be accurate only in the vicinity of of P o and 1 -P o (see equation 4). As in [START_REF] Grandvalet | Support Vector Machines With a Reject Option[END_REF], we built a hinge loss function tangent to neglog-likelihood at f + = log[P o /(1 -P o )] and at f -= log[(1 -P o)/P o ]. The decision rule can be expressed in terms of the function f as follows:

d * (x, f ) =      +1 if f (x) > f + , -1 if f (x) < f -, 0 otherwise . ( 5
)
where the output f (x) of the SVM corresponds to the decision d * , minimizer of the risk (3).

Optimization

Introducing the double hinge loss results in an optimization problem that is similar to the standard SVMs problem. Let C be the penalization parameter to optimize in the training phase, D = C(1 -P o ). The optimization problem is:

               min f,b,ξ,η 1 2 f 2 + C n i=1 ξ i + D n i=1 η i , S.c y i (f (x i ) + b) ≥ τ -ξ i i = 1, . . . , n , y i (f (x i ) + b) ≥ -η i i = 1, . . . , n , ξ i ≥ 0 , η i ≥ 0 i = 1, . . . , n . (6)
where τ = (P o log(P o ))/(1 -P o ) -log(1 -P o ). To compute the solution of 6, we use an active set algorithm following a strategy that proved to be efficient for standard SVMs. We refer to [START_REF] Grandvalet | Support Vector Machines With a Reject Option[END_REF] for a detailed discussion.

Results and discussion

The input of the classifier is a set of 15 element vectors x i , representing the ECG beats. To fine tune a global classifier, we have selected 10 records containing the largest number of PVC beats. This is done to provide a relative balance between the positive class and the negative class in the training set. we used the first 5 min of data taken in the selected records. This practice conforms to the AAMI-recommended procedure allowing the usage of at most 5-min section from the beginning of each patients recording for training [START_REF] Hearth | Recommended Practice for Testing and Reporting Performance Results of Ventricular ArrhythmiaDetection Algorithms[END_REF]. From each segment of 5 min, were selected 150 beats leading to a training set of 1500 vectors. During the test phase, only the last 25 min in each type are used. Hence, the testing data are never part of any training data.

To explain the effectiveness of the classification methodology and the impact of the reject option, we simply apply standard statistical data visualization technique by projecting the high-dimensional data onto a 2D space. Figure 7 shows an example of the reject region produced by the SVM classifier when our algorithm is applied to 221 and 214 containing together a negative classe of 3574 beats (1543 of Normal beats and 2003 of L beats) and 652 of V beats considered as a positive classe. On this figure, we can see the reject thresholds given by equation [START_REF] Talbi | PVC discrimination using the QRS power spectrum and self-organizing maps[END_REF]; f + = 0.8473 and f -= -0.8473.

A complete description of a such classifier is given by the error-reject tradeoff (the error rate E against the reject rate R). A typical tradeoff curve is displayed on figure 8(a). For small rejection rate (R < 1.9%), the curve shows a high slope. This means that the majority of the rejected patterns would have been erroneously classified. We point out that the advantage of classifying with rejection can be judged by the error-reject tradeoff. figure 8 positive class.

A lot of studies dealing with heartbeat classification using the same database are present in the literature. All previous works used a two decision rule (No Rejection). The results we obtained with No Rejection (NR) show to be competitive to average results reported by other authors separating beats whose origin is from the SA node from premature ventricular contractions (see Table 1). On the same table, we reported our achieved results corresponding to the Minimal Classification Cost (M C c ). Note however that the main advantage of our approach is the ability to optimize the classification cost.

Table 1: Comparison of the performances of our classifier with others using Artificial Neural Network (ANN) or hybrid Genetic Algorithm-SVM (SVMGA) in terms of sensitivity.

Method

Model classifier S e (%) [START_REF] De-Chazal | A Patient Adapting Heartbeat Classifier Using ECG Morphology and Heartbeat Interval Features[END_REF] ANN 94.4 [START_REF] Khazaee | Classification of electrocardiogram signals with support vector machines and genetic algorithms using power spectral features[END_REF] SVMGA 98.91 [START_REF] Krasteva | QRS Template Matching for Recognition of Ventricular Ectopic Beats[END_REF] ANN 98.4 Our algorithm with NR -98,2 Our algorithm for M C c -99

Conclusion

In this paper, we proposed an automated heartbeat classifier, using Support Vector Machines with an embedded reject option. The proposed system accomplishes preprocessing, feature extraction and recognition tasks for ectopic heartbeats detection.

A set of features including frequency informations, RR intervals, QRS morphology and AC power of QRS detail coefficients is exploited to characterize each beat. The resultant feature vector representing each ECG heartbeat is used as input of the modified SVM classifier.

For this purpose, a cost-sensitive reject rule for SVMs is set. The main feature of our system is that the decision rule uses dynamic reject thresholds following the cost of rejecting a sample and the cost of misclassifying a sample in order to optimize the classification cost. This approach is important in medical applications when rejection may even be more preferable than risking misclassifications.

Our results shown above illustrate a good error reject tradeoff and indicate an error rate of 2.8% with no rejection and less than 1.2% for the minimal classification cost. These results are competitive to other published studies.

Figure 1 :

 1 Figure 1: Different formes of ECG beats.

Figure 2 :

 2 Figure 2: a: Example of R peak detection and QRS segmentation. b: Zoom on QS segments of PVC beat (left) and normal beat (right).

Figure 3 :

 3 Figure 3: RR histogram for N, L, R and V beats

Figure 4 :

 4 Figure 4: Power spectrum of different beats (a) and the relation between frequencies and power levels (b). For better representation, we limited the frequency axis to f s /4.

Figure 6 :

 6 Figure 6: Expected losses against posterior probabilities P (x).

  d

Figure 7 :

 7 Figure 7: scatter plot showing the reject region induced by the reject thresholds in correspondence to the costs of misclassifying and rejecting samples. Positive cases are represented by black circles and negative cases by red triangles. The pink region corresponds to classification Y = +1, and white regions to classification Y = -1. The blue and the light blue regions correspond to no decision. The lines +0.8473 and -0.8473 correspond respectively to f + and f -and the line 0 corresponds to f (x) = 0 or P (Y = 1 | X = x) = 0.5.

Figure 8 :

 8 Figure 8: Error-reject tradeoff (a) and the classification cost against the reject rate (b).

Where x ij is the j th component of the i th feature vector . x j and σ x j are the mean and the standard deviation of the j th component of the feature vectors, and tansig(.) the hyperbolic tangent sigmoid transfer function. The objective is to make each j th component to be normal distributed with zero-mean and unit standard deviation.

Classifier Model

Bayes rule

A discriminant function f : X → R classifies an observation x ∈ X into one of two classes, labeled +1 or -1 . Viewing f (x) as a proxy value of the conditional probability P (x) = P(Y = 1|X = x), one is less confident for small values of | f (x) | corresponding for P (x) around 1/2. The strategy used in this work is to report sgn(f (x)) = +1 or -1 if |f (x)| exceeds a threshold δ corresponding to P (x) = P o and no decision otherwise. Assuming that the cost of making a wrong decision is c and that of withholding a decision is r, the smallest risk R for a decision d is: R(d) = min{cP (x), c(1 -P (x)), r}.

(3)


Ion Candel
email: ion.candel@edf.fr

Cornel Ioana
email: cornel.ioana@gipsa-lab.inpg.fr

David Lemon
email: dlemon@aslenv.com

Bertrand Reeb
email: bertrand.reeb@edf.fr

USING WIDE BAND SIGNALS FOR OBSTACLE PATH CORRECTION IN ACOUSTIC SCINTILLATION FLOW METERS

Keywords: flow rate measurement, scintillation, CHIRP, wide band, matching filtering

This paper deals with the challenges of water flow measurement in hydroelectric plants using acoustic scintillation. One of the drawbacks of ultrasonic signal propagation through water is that measurements require relatively clean water conditions. In water with significant amount of impurities, the level of transmitted signals may not be sufficient for an accurate flow measurement. In terms of signal processing, the effect of impurities is represented by a random loss or fading of received signals. In order to improve the measurement, the wide band signals can be used in order to improve the quality factor of the signals by adaptive filtering using a static mode reference. Results show that, in real configurations, the interference due to the obstacles can be reduced.

INTRODUCTION

Acoustic scintillation flow metering promises to offer an elegant and efficient non-intrusive solution to flow metering by, in general: mounting the transducers on the outer walls of the pipe, transmitting a signal to be affected by the flow, and analyzing the received signal in order to obtain the flow parameters [START_REF] Lynnworth | Ultrasonic Measurements for Process Control: Theory, Techniques, Applications[END_REF]. However, a series of tests in water containing impurities showed that the resulting insufficient signal levels may cause the system to stop the measurement.

ACOUSTIC SCINTILLATION DRIFT

The acoustic scintillation drift technique measures the water flow velocity perpendicular to a number of acoustic paths established across an intake to a turbine [START_REF] Lemon | Improvements to the Accuracy of Discharge Measurements by Acoustic Scintillation resulting from Revisions to Data Processing Procedures[END_REF].

Short ultrasound pulses are sent from the transmitting arrays, on one side, to the receiving arrays on the other, at a rate of approximately 250pulses/second [START_REF] Asl Aqflow | ASFM Operations Manual[END_REF]. Turbulence embedded in the flow causes amplitude fluctuations in the acoustic signals (Fig. 1).

USING LARGE BAND SIGNALS ON SCINTILLATION

Amplitude variations found at the receiving end of the scintillation acoustic path can be calculated using a different technique, based on wide band signals, such as linear frequency modulations or CHIRPs. Instead of using a pulse package of a single carrier frequency, the frequency within the burst is changed (swept) throughout the duration of the transmission, from one frequency to another [START_REF] Skolnik | Introduction to radar systems[END_REF]. In our application, we have constructed a wide band signal in order to replace the impulse package found in the emission stage. A matched filter is obtained by correlating a known signal, or reference, with an unknown signal to detect the presence of the reference in the unknown signal. This is equivalent to convolving the unknown signal with a conjugated time-reversed version of the template.

Due to the high spectral content of a large band signal, the autocorrelation between two large band signals will always have a significant maximum for the main lobe.

Considering a CHIRP signal received when an obstacle blocks the acoustic path, the spectral content of the transmitted signal is still found when performing a time frequency analysis. Figure 2 illustrates the spectrogram for a received chirp signal in the case of an obstacle: For the same measurement, it is possible to view also the spectral content of a CHIRP signal received with no obstacle blocking the acoustic path, as shown in figure 3. Using this property, we select a reference signal which will be the first term of the correlation. This reference is the received CHIRP signal which has traveled along the acoustic path during a clear water test. Next, we conduct an obstacle test that attenuates the acoustic path. For each received pulse we then calculate the corresponding correlation and extracted a maximum from each one. The extraction of maximums of correlations provides a new time series, which can be compared with the one based on envelope detection.

Figure 4 shows the difference between the time series calculated using the correlation and envelope detection: the influence due to the impurities/obstacles presented in the acoustic path is greatly reduced. The two time -frequency representations in figures 2 and 3 show the potential of wide band signals in dealing with impurities/obstacles encountered across the acoustic path. It should be noted that these tests were carried out using signals with low levels of power, hence the low signal to noise ratio for the signal in figure 2.

Another advantage of using large band signals arises from the "personalized" spectrum content: the algorithm can be rendered robust to outside perturbation and the complexity of the equipment can be reduced.

Figure 1 .

 1 Figure 1. Diagram of the acoustic scintillation principle. If the two paths are sufficiently close (Δx), the turbulence effect on the signals remains virtually the same, except for a time delay, Δt. This time delay corresponds to the peak in the time-lagged cross-correlation function calculated for Signal 1 and Signal 2 . The mean velocity perpendicular to the acoustic paths becomes v mean = Δx/ Δt.

Figure 2 .

 2 Figure 2. Spectrogram of a received CHIRP signal in the case of an obstacle.

Figure 3 .

 3 Figure 3. Spectrogram of a received CHIRP signal for no obstacle.

Figure 4 .

 4 Figure 4. Time series computed using the correlation (green) and envelope detection (blue).

  This paper deals with an improvement of a flow rate measurement method in the context of impurities/obstacles, as impurities/obstacles inherent in the water flow pose a threat to the overall accuracy of the measurement results.Signal processing techniques can eliminate this drawback by constructing a new type of transmitting signal and by using a new method to quantify the turbulence embedded in the water flow. The method employs wide-band signals (CHIRPS) and the matching filtering. Worst case scenario tests have shown encouraging improvements in time series amplitudes.Future work will consist in adding impurities of different sizes in the water flow and quantifying their influence on the time series, as well as on the water flow measurement accuracy.


Bertrand Reeb
email: bertrand.reeb@edf.frion.candel@edf.fr

Ion Candel

Gilles Proulx
email: gilles.proulx@hydro.qc.cadlemon@aslenv.com

Eng. PhD. Eng David Lemon

Investigation of the Performance of Acoustic Scintillation Flow Meter when Turbulence Levels are Low

Electricité de France (EDF) is funding a 3-year PhD work on acoustic scintillation flow metering, in association with Hydro-Québec and ASL AQFlow. The PhD project aims at improving the discharge estimation when hydraulic conditions are not ideal for acoustic scintillation measurements, as well as providing a better understanding of the effects of sources of interference which can be encountered during the measurement process and can cause inaccuracies in the velocity estimation. In order to achieve these objectives, a fast and portable data acquisition system was set up, which relies on high speed acquisition cards. Each four channel acquisition card can be connected in parallel with up to three others, thus forming a high speed multi channel data acquisition (DAQ) system. A first test of this DAQ system was performed at one of Hydro-Québec's hydroelectric plant. As the ASFM replaced the existing trash rack elements, special equipment was designed and built to create turbulence in the flow necessary for the acoustic scintillation to operate. HPP performance tests by Hydro Quebec were duplicated with the EDF measurement system which recorded the acoustic scintillation raw signals. Using the resources available at the time, the acquired time series have a lower resolution than those provided by the Acoustic Scintillation Flow Meter (ASFM). However, valid velocity estimates are possible even at this lower resolution.

I.

Introduction

Water flow estimation in hydro power plant intakes provides important data used to evaluate turbine efficiency and to assure compliance with regulatory release requirements. A flow metering system, which is robust, adaptable and precise, is therefore needed.

As far as low head HPPs, with small converging intakes, are concerned, acoustic scintillation flow metering seems an interesting an innovative solution to the discharge measurement issue.

EDF's initiative to use the Acoustic Scintillation Flow Meter (ASFM) is part of a continuing process of improving the performances of the hydroelectric plants. Previous measurement sessions performed with the ASFM revealed that flow effects such as transducer support vibrations, grid-induced wakes and low turbulence levels could negatively influence the results. The EDF PhD thesis project was started in February 2011 with the objective to improve flow metering by the ASFM, and to provide solutions for existing problems using advanced signal processing techniques.

This paper deals with the challenges that arose in the first months of the thesis, as well as with the difficulties experienced during the measurements at one of Hydro Quebec's HPP.

II. The Acoustic Scintillation Flow Meter

The Acoustic Scintillation Flow Meter (ASFM) uses a technique called acoustic scintillation drift to measure the flow velocity perpendicular to a number of acoustic paths established across the intake to the turbine. Short (16 µsec) pulses of high-frequency sound (307 kHz) are sent from transmitting arrays on one side to receiving arrays on the other, at a rate of approximately 250 pings/second [1]. Fluctuations in the amplitude of those acoustic pulses result from turbulence carried along by the flow. The ASFM measures those fluctuations (known as scintillations) and from them computes the lateral average (i.e. along the acoustic path) of the velocity perpendicular to each path. In its simplest form, two transmitters are placed on one side of the measurement section, two receivers at the other (Fig. 1). The signal amplitude at the receivers varies randomly as the turbulence along the propagation paths changes with time and the flow.

If the two paths are sufficiently close (Δx), the turbulence remains embedded in the flow, and the pattern of these amplitude variations at the downstream receiver will be nearly identical to that at the upstream receiver, except for a time delay, Δt. This time delay corresponds to the peak in the time-lagged cross-correlation function calculated for Signal 1 and Signal 2. The mean velocity perpendicular to the acoustic paths is then Δx/Δt.

Using three transmitters and three receivers at each measurement level allows both the magnitude and inclination of the velocity to be measured. The ASFM computes the discharge through each bay of the intake by integrating the horizontal component of the velocity over the cross-sectional area of the intake. In a multi-bay intake, the discharges through each bay are summed to compute the total discharge. The ASFM instrument consists of five major components: transmitting (Tx) and receiving (Rx) transducers and underwater cabling, switching canisters, surface connection cables, a data acquisition Surface Unit, and a PC computer with the user interface for controlling and operating the ASFM (ASFM Link software). The ASFM typically is divided into 3 groups for operation in multi-bay intakes; each group consists of a set of transducers, two switching canisters, and cabling. A single ASFM system can have up to 30 paths with 3 groups (typically 10 paths per bay in a 3-bay intake -Fig. 2). One path from each group may be sampled simultaneously, provided the path lengths are equal.

III. The multi-channel acquisition system

The scientific investigations done throughout the thesis are based on signals recorded from the ASFM and the first task was to identify the appropriate points of measurement inside the ASFM. In order to study closely the operation mode of the ASFM, as well as to replicate the velocity algorithm inside the ASFM, several signals had to be viewed at the same time, which exceeded the capacities of conventional oscilloscopes. A multi-channel acquisition system had to be designed.

In addition to this problem, a new challenge arose from the fact that the scintillation pulses are being transmitted with a frequency of approximately 250Hz for a period varying from 33 seconds to 120 seconds for each level [6]. Therefore, the acquisition system has to record these pulses for long periods of time and in a format that can be read easily. Furthermore, the acquisition system must be compact as the measurements often take place in zones with too little space to deploy a large and complicated acquisition system.

The solution has been found in the form of a portable computer and as many as four HS4DIFF data acquisition units provided by TiePie Engineering®, as illustrated in Fig. 3.

Each HS4DIFF unit connects to the computer via an USB connection and a synchronizing cable is used to connect two or more units together. A computer application controls the acquisition process by creating a combined instrument from two or more HS4DIFF units.

The combined instrument can support up to five HS4DIFF units. Several options are available for recording the data such as Data Loggers, DiskWriter and AutoDisk. The preferred choice is DiskWriter because it records data easily readable with Matlab® [2]. After extensive testing in the laboratory, the system was deployed at one of Hydro Québec's HPP. Two HS4DIFF units were used to record the analog form of the signals, both the transmitted signals and the signal found at the input of the analog to digital converter. The first advantage of such a system is the ability to visualize in the same time the scintillation signals of the ASFM unit and the signals at different stages along the signal path in the electronics, including the envelopes of the received scintillation pulses located at the input of the analog to digital converter (ADC). Most acquisition systems that provide a number of acquisition channels higher than four have low sampling frequencies (up to several hundred kilohertz) or low resolution ADCs.

However, even with such extended capabilities, the system was unable to record all the 250 scintillation pulses which are transmitted in one second. This inconvenient means that the calculated time series have a shorter length compared to the ones provided by the ASFM software and that velocity estimations are more difficult to obtain.

IV.

The analysis of the data recorded at Hydro Québec's hydro power plant

The most important part of the data recorded were the received pulses from the three elements of the scintillation transducers (Fig. 5). These pulses are used to construct the time series which characterize the turbulence generated by the water flow and thus calculate the flow velocity. In Fig. 5, each subplot represents the signal recorded for each element of the acoustic transducer. Each receiving element receives, at the expected moment in time, the corresponding pulse from its counterpart at the transmission level. Since the elements are close, the angles of the beams emitted allow each element to receive the signals transmitted for the other two, hence the succession of the three pulses for each element (the fourth pulse present at the beginning of the signals corresponds to the previous transmission). This phenomenon does not introduce errors in the velocity estimation algorithm because the pulses are separated in time and the necessary identification subroutine is built into the ASFM software.

The signals shown in Fig. 5 display an interesting aspect: on all subplots there are three smaller pulses in the (1 ÷ 2.5) msec. time interval. These pulses (called echoes) are actually reflections from previous pings which are reflected initially at the receiving side (Rx) and then travel to the emission side (Tx) to be once again reflected back to Rx.

These echoes could provide an alternative method to determine the travel time of the signals between the transducers. When the emitted pulse for ping "n" reaches Rx, a small fraction of the signal's energy is reflected back at Tx. After reaching Tx, the signal is reflected again and travels to Rx, after ping "n+1", as shown in Fig. 6: frequency is 10MHz, but it can be easily increased to 50MHz, which would translate to calculating a time of flight between the transducers with an error of 20 nanoseconds.

0 0.5 1 1.5 2

Velocity estimation using the recorded data

These on-site measurement sessions were the first in which real data was recorded with the DAQ system and velocity estimations were possible. The ASFM used for data recording had to manage up to 10 levels. Continuous recording of data ensured that for each level an almost equal amount of pings was recorded. Due to the time needed for the data to be written on the disk, a large number of received pings were not recorded, thus making velocity estimations more complicated [3]. This is due to the recording process performed by the oscilloscopes: data blocks are copied into the oscilloscope's buffer memory and then sent to the computer via the USB port. During this time, the acquisition is inactive and the subsequent scintillation pulses are not recorded. The ASFM unit emits 250 pulses / second, therefore, we can say that the turbulence embedded in the flow is digitized 250 times /second, thus a sampling frequency of 250Hz. For the signals recorded, the actual time series had an average sampling frequency of 15Hz. The raw time series obtained from these recordings had to be processed using an adaptive interpolation method and a filtering technique in order to correct the spectral modifications introduced by the interpolation.

Figure 8. Time series for the UPSTREAM element provided by the ASFM and the recorded data.

The interpolation technique used is a piecewise polynomial form of the cubic "spline interpolant" to the data values [START_REF]A Practical Guide to Splines[END_REF]. This technique has been used successfully to calculate the envelopes of the received pings and gives very good results with interpolating under-sampled data. The recorded signals were digitized with a frequency of 10 MHz, so no interpolation was needed in order to calculate the envelopes of the received pulses. The time series resulted from the envelope detection had to be interpolated in order to be able to estimate the time delays. The interpolation inserts more points between the points of the under-sampled data, thus enhancing the time resolution. However, too much interpolation can have undesirable effects, such as enhancing unwanted spectral components. In order to avoid such drawbacks, the interpolation has been adapted to the spectral content of each time series in such way that the spectral changes are being kept to a minimum. Reducing the negative effects of the interpolation is done by comparing the resulting spectrum of the time series after interpolation (more precisely the high frequency content) with the spectrum obtained for the time series provided by the ASFM software. The latter have a spectral content which stretches in theory up to 125 Hz (according to the Nyquist -Shannon sampling theorem), so any spectral components found above this value on the interpolated time series spectrums are filtered. Moreover, some spectral components which do not exceed the 125Hz value can dwarf the other components, rendering the correlation to output a zero delay between the time series. These components were filtered using a Butterworth band-pass filter. Throughout the entire process, the ASFM time series provided a reference for the results.

Velocity estimations were performed both on the interpolated time series and the ASFM time series using an algorithm derived from a diagram provided by AQFlow.

Fig. 9 shows the results of the velocity estimation using the ASFM time series and the interpolated data. The two traces in Fig. 9 show that velocity estimations are possible to extract from the recorded data and the nature of the profile can be the same. However, there are points on the profile (especially for levels 6 and 9) where the interpolation cannot provide an exact value for the velocity. This effect is due to the relatively low number of pings recorded (for these levels the actual sampling rate was lower than 15Hz), but also due to the turbulence levels not producing enough variation in the amplitudes.

From a signal processing point of view [START_REF]Optimum Signal Processing. An Introduction[END_REF], the time series calculated present statistic parameters (mean and variance) which resemble the ones of a Gaussian noise (Fig. 10). Practically, there is no perfect noise waveform and therefore correlations are possible, but the actual correlation peak (which corresponds to the time delay between the two time series) can be comparable to the adjacent peaks, making the delay estimation unsure. Fig. 10 shows an example of two correlations from time series with low and high turbulence levels: In the first subplot of Fig. 10 the two most predominant peaks correspond to two values of time delay (and therefore velocity) and may create confusion. The ideal situation is to obtain a correlation waveform which resembles the one displayed in the second subplot of Fig. 10.

V. On-site ASFM measurements

A first set of tests were carried out on unit 21 of the HPP under test with the trash racks removed and replaced with the ASFM transducer frame supports. It was evident that the turbulence level was very low. In order to improve the quality of the measurement, it was decided to install a device to generate more turbulence. A series of chains were attached to a barge, forming a chain curtain. The chains were maintained vertically using a heavy weight attached to the lower end. The spacing between the chains was 30 cm for the entire width of the intake bay (Fig. 11). As it can be seen on the ASFM time series in Fig. 12, and the spectral comparison in Fig. 13, the turbulence level has been significantly increased. However, even with the higher turbulence level generated by the chains, the spectral content was yet low compared to a good measurement condition and therefore the system can benefit from further improvements.

Fig. 12. Turbulence level data of an acoustic path with and without chain curtain

.

Fig. 13. Comparison of spectral amplitudes of time series with and without chain curtain

Another way to view the effect of the turbulence level is to look at the velocity profile calculated from measurements with and without the chain curtain. Fig. 14 shows that the velocity profile without the chain curtain is inconsistent, especially in the middle of the measurement section. On the other hand, by using the chain curtain, the velocity profile becomes smoother.

VI. Conclusions

A high speed data acquisition system has been designed and tested in order to acquire raw acoustic scintillation signals used in ASFM flow meter. The data proved useful for direct visualization of the transmitted signals, but velocity estimation proved more difficult to obtain. Signal processing techniques were used to compensate the amount of data that was not collected and in several cases the results matched with the ones provided by the ASFM data series. Concerning the high speed acquisition, the system has been later improved and configured in such a way that has allowed a complete acquisition of the entire scintillation transmission sequence, providing better results of velocity computation.

Future work will concentrate on simulating various causes of ASFM unfavorable situations, such as signal interference, and on applying new signal processing techniques (time-phasefrequency, wavelets decompositions, etc.) to verify if they can improve the actual acoustic scintillation flow metering concept. In the case of low turbulence time series, two possible techniques will be considered: the use of the "Singular Value Decomposition" method to extract relevant correlation data from the time series and the use of the ambiguity function for a more precise localization in time of the correlation peaks.

Fig. 1 :

 1 Fig. 1: Schematic representation of acoustic scintillation drift.

Figure 2 :

 2 Figure 2: Components of a typical multi-bay ASFM.

Figure 3 .

 3 Figure 3. The HS4DIFF recording unit.

Figure 4 .

 4 Figure 4. Data is being recorded from the ASFM using two HS4DIFF units (on the right).

Figure 5 .

 5 Figure 5. Recorded scintillation pulses.

Figure 6 .

 6 Figure 6. Echo forming between two scintillation pings.

Figure 7 .

 7 Figure 7. Pings and echoes used for water temperature estimation.

Figure 9 .

 9 Figure 9. Velocity profiles obtained from the data recorded and ASFM.

  Figure 10. Comparison between correlation of ASFM time series with low and high turbulence levels. (plus grand)

Fig. 11 .

 11 Fig. 11. Chain curtain used for creating turbulence in the intake bay.

  windows belonging to high turbulence time series

Fig. 14 .

 14 Fig. 14. Velocity profiles calculated with and without chain curtain.

References 1 .

 1 "Improvements To The Accuracy Of Discharge Measurements By Acoustic Scintillation Resulting From Revisions To Data Processing Procedures", D.D. LEMON, ASL AQFlow Inc., Victoria, BC, Canada, D. R. TOPHAM, Topham Scientific Consulting, Victoria, BC, Canada and D. BILLENESS, ASL AQFlow Inc., Victoria, BC, Canada. 2. Tech report -EDF DTG, Ion CANDEL, September 2011.


Richard J Nowakowski

Gabriel Renault
email: gabriel.renault@labri.fr

Emily Lamoureux

Stephanie Mellon

Timothy Miller

The Game of timber!

We analyse timber, a game played on a graphs. We find the P positions for both normal and misère play on paths and show how to win the game. In passing, we also show a correspondence with Dyck paths, the Catalan and Fine numbers. We present an algorithm for winning the Normal Play game on trees.

Introduction

The game of timber arises out of the coming together (loosely) of the work of C. Mynhardt's work on "altitude" and "depression" on graphs, see [START_REF] Mynhardt | Trees with depression three[END_REF][START_REF] Mynhardt | Altitude of regular graphs with girth at least five[END_REF] for examples, and an offshoot of the game toppling dominoes. The game timber is played on a directed graph, with a domino on each edge. The edges will be denoted by (x, y) with the arc starting at x and ending at y. There are two players, herein called Alice and Bob, moving alternately. On a move, the player chooses a domino on some edge, say, (x, y) and topples it in the direction of y. (This is the only time the direction of the edge is important.) The domino then topples the dominoes on the edges incident with y, independent of whether the edge is directed into or away from y, and the process of toppling the dominoes continues until no more dominoes topple. The toppled dominoes and corresponding edges are removed from the graph. More formally, let G be a graph, a timber position based on G will be denoted by G t which will be G with its edges directed. Note that for a given G there will be many possible G t . Given a graph G, a corresponding G t and an arc (x, y), let Y be the set of all vertices that can be reached by a path in G starting at y and not having x as the second vertex. Choosing to topple the domino on the arc (x, y) removes Y and all edges incident with Y from G and G t . For example, see Figure 1. We invite the reader to find the unique winning move in this position.

x toppling x leaves We will attempt to make the paper self-contained but the reader can consult [START_REF] Albert | Lessons in Play[END_REF][START_REF] Berlekamp | Winning ways for your mathematical plays[END_REF] for more details on game theory. Under Normal play, the first player who cannot move is the loser; under misère play this player would be the winner. timber is an impartial game, since both players have the same moves available, i.e., Alice and Bob can start with any domino. In impartial games there are just two outcome classes: N -positions in which the Next player to move can force a win; and P-positions in which the next player cannot force a win but the Previous player can. In Normal play, the end position is a P-position (since the next player has no move at all and therefore has no way to force a win) and it would be an N -position in misère play. Note that the outcome classes can be determined recursively starting with the end of the game. From a position G t , if a player can move to a P-position then G t is an N -position, otherwise G t is a P-position. We write, o + (G t ) for the outcome class of the position G t under Normal play and o -(G t ) under misère play.

Since timber could be played on a set of disconnected components we should be using the Sprague-Grundy theory for impartial games, however, we found the game was too complicated for us to solve at this level of detail, even on a path. We do completely characterize the outcome class for all G t , where G is any connected graph, under Normal play and find the misère outcome classes where G is a path or a 2-connected graph.

Our first result shows that, for most graphs, the outcome class is very easy to calculate.

Theorem 1.1 If G is a connected graph with a 2-connected subgraph then for any G t , o + (G t ) = N . If G is a 2-connected graph then for any G t , o -(G t ) = P.
Proof. Let G be a connected graph and G t be some timber position of G. Let xy be an edge of G in the 2-connected subgraph and (x, y) the corresponding arc in G t . Toppling the domino on (x, y) causes all the dominoes to be toppled. In Normal play, this is a winning move. In misère play it is a losing move and if G is 2-connected the next player only has losing moves. Thus, in Normal Play, only graphs which are trees need be considered. The paper is divided into two main sections. In Section 2, timber played on a path is analyzed. Although the algorithm for trees could be used to analyze the game on paths, the path analysis uncovers relationships to Dyck paths, Catalan and Fine numbers and is worthy of independent investigation. The Normal play version on paths was the original game analyzed by Lamoureux, Mellon, and Miller as a first year class project at Dalhousie University, supervised by the first author. The game was called toppling peaks. See Section 2 and Figure 2 to see how the game got its name and why a possible relationship to 'altitude' and 'depressions' of graphs could be inferred. Nowakowski and Renault then found the misère analysis and generalized the game to graphs in general. In Section 3, we give a reduction algorithm to determine the outcome class of a tree. Surprisingly, nimbers are required even though we are only investigating outcome classes. In the last section we raise some questions. timber on paths, also called toppling peaks, arose as an impartial variant of the partizan game toppling dominoes [START_REF] Fink | Toppling Conjectures[END_REF] which is played with the same set-up: a line of dominoes, each domino marked L or R where in toppling peaks an L can only be toppled to the left and an R to the right. See Figure 2.

After some preliminary analysis of timber, using CGSuite [START_REF]Combinatorial Game Suite (CGSuite[END_REF], we quickly determined the number of P-positions of length n = 1, 2, . . . , 10 is 0, 1, 0, 2, 0, 5, 0, 14, 0, 42, respectively. Ignoring the 0s, the On-line Encyclopedia of Integer Sequences [START_REF] Sloane | The On-Line Encyclopedia of Integer Sequences[END_REF] suggests that the resulting sequence is the Catalan numbers, 1 n+1 2n n . The first two references are [START_REF] Deutsch | A bijection of Dyck paths and its consequences[END_REF][START_REF] Guy | The nesting and roosting habits of the laddered parenthesis[END_REF]. The secondnesting and roosting habits of the laddered parenthesis-suggests a pairing approach that is very useful in actual play. The first reference-A bijection of Dyck paths and its consequences-gives rise to the name of the game and hints at the solutions.

A position can be represented visually on a 2-dimensional graph: start at (0, 0) and let an L be a line joining the lattice points (x, y) and (x+ 1, y + 1) and R be the line joining (x, y) and (x + 1, y -1). We'll sometime refer to this as the peak representation. For example, all the representations of LLRL are given in Figure 2. A Dyck path of length 2n is one of these paths that also ends at (2n, 0) and which never goes below the x-axis. See Theorems 2.3 and 2.4.

For normal play, Theorem 2.3 shows that P-positions are exactly those positions with Dyck path representations. For misère play, there are two types of P-positions. Theorem 2.4 shows that one kind, a subset of the Normal play P-positions, correspond directly to the Dyck paths with no 'peaks' of height 1, i.e., a Dyck path with no sub-path (j, 0), (j + 1, 1), (j + 2, 0) for any j. The misère P-positions of length 2n + 1 are formed from these P-positions of length 2n by the addition of an R at the start or an L at the end. The extra P-positions of length 2n + 2 are obtained by adding both the L and the R dominoes. In general, the number of misère P-positions of length 2n is the number of Dyck paths of length 2n with no peaks of height 2 before the path returns to 0, but this is a counting argument since the representation of a length 2n P-position does not have to be a Dyck path and can have a peak of height 2.

The Analysis.

We will refer to the two players as 'Alice' and 'Bob' usually with Alice to go first in a given position.

The main tool for characterizing outcome classes are the following theorems. Then A is the set of P-positions and B is the set of N -positions.

We now translate the intuitive 'peak' approach into one where it is easier to give proofs and more useful when we consider trees.

Given the alphabet {L, R}, for a word w, let |w| L be the number of Ls in w, |w| R the number of Rs in w and w [i,j] the subword w i w i+1 • • • w j . Let W P be the set of words w such that for any i,

|w [0,i] | L |w [0,i] | R and |w| L = |w| R .
Theorem 2.3 In Normal play, the P-positions of toppling peaks are exactly those which correspond to word of W P .

Proof. Let w ∈ W P be a position and n its length. First let's assume that Alice plays on an R domino, say w i . Then Alice has moved to the position

w [0,i-1] . As w ∈ W P , |w [0,i] | L |w [0,i] | R . Hence |w [0,i-1] | L = |w [0,i] | L |w [0,i] | R = |w [0,i-1] | R + 1, and |w [0,i-1] | L = |w [0,i-1] | R . So w [0,i-1] / ∈ W P .
If Alice plays on an L domino, instead, say w i . Then Alice has moved to the position

w [i+1,n-1] . As w ∈ W P , |w [0,i-1] | L |w [0,i-1] | R . Hence |w [i+1,n-1] | L = n 2 -|w [0,i] | L = n 2 -|w [0,i-1] | L -1 n 2 -|w [0,i-1] | R -1 = |w [i+1,n-1] | R + 1, and |w [i+1,n-1] | L = |w [i+1,n-1] | R . So w [0,i-1] / ∈ W P . Let w / ∈ W P be a position. Assume |w| L = |w| R . Let i = min{0 k n -1 | |w [0,k] | L < |w [0,k] | R }. w k is an R domino, and |w [0,k] | L = |w [0,k] | R -1. By toppling w k , Alice moves to w [0,k-1] ∈ W P . Assume now |w| L = |w| R . Without loss of generality, we can assume |w| L < |w| R . Let i = min{0 k n -1 | |w [0,k] | L < |w [0,k] | R }. w k is an R domino, and |w [0,k] | L = |w [0,k] | R -1. By toppling w k , Alice moves to w [0,k-1] ∈ W P .
Note that w ∈ W P iff the corresponding position produces a Dyck path.

We now consider Misère play. Let SW P be the set of words w such that w ∈ W P and ∀w 1 , w 2 ∈ W P , w = w 1 LRw 2 . We define X = (SW P \{∅}) ∪ {Rw | w ∈ SW P } ∪ {wL | w ∈ SW P } ∪ {RwL | w ∈ SW P }. We note w the word obtained from w after removing the first character if it is an R and the last one if it is an L. Theorem 2.4 In misère play, the P-positions of toppling peaks are exactly those which correspond to word of X.

Proof. Let w ∈ X be a position. Assume w ∈ (SW P \{∅}). From the normal play analysis, we know that Alice cannot move to a position in SW P ⊂ W P . Assume Alice can move to a position Rw 0 with w 0 ∈ SW P . Then it follows that w = w 1 LRw 0 for some w 1 . As w, w 0 ∈ W P then w 1 ∈ W P , which is not possible since w ∈ SW P . Similarly, we can prove Alice has no move to a position of the form w 0 L or Rw 0 L with w 0 ∈ SW P . Similarly, we can prove Alice has no move to a position in X from a position in X.

Let w /

∈ X ∪ {∅}. Assume w ∈ W P . Then there exists w 1 , w 2 ∈ W P such that w = w 1 LRw 2 , and we can choose them such that w 2 ∈ SW P . From w, Alice can move to Rw 2 ∈ X. Similarly, we can prove Alice has a move to a position in X from a position in ({Rw

| w ∈ W P } ∪ {wL | w ∈ W P } ∪ {RwL | w ∈ W P })\X. Now assume w [0,1] = RR. Alice can move to R ∈ X.
Now assume w is none of the above forms. w starts with an L and ends with a R, and is not in W P , so Alice has a move from w to a position w 0 ∈ W P \{∅}. Without loss of generality, we can assume it is by toppling a domino leftward. If w 0 ∈ SW P , the same move from w leaves the position w 0 ∈ X or w 0 L ∈ X. Otherwise, there exists w 1 , w 2 ∈ W P such that w 0 = w 1 LRw 2 and we can choose w 2 ∈ SW P . Alice can then move from w to Rw 2 ∈ X or Rw 2 L ∈ X.

Counting positions

Surprisingly, there are few games for which the number of P-positions are known. In fact, the authors only know of [START_REF] Hetyei | Enumeration by kernel positions[END_REF] where the game was developed to relate the number of P-positions to Bernoulli numbers of the second kind. Even for the basic game of nim with n pieces, it is not presently known how many P-positions there are (see [START_REF] Sloane | The On-Line Encyclopedia of Integer Sequences[END_REF] sequence A048833).

Let W P (n), SW P (n) and X(n) be the number of P-positions of length n respectively: (1) in normal play; (2) in misère play where the position is a Dyck path; and (3) in misère play in general.

A peak at height k on a Dyck path is a point (j, k) of the path that is preceded by the point (j -1, k -1) and followed by the point (j + 1, k -1). The Fine number, F n , is the number of Dyck paths of length 2n without peaks at height 1.

Theorem 2.5 Let n be a nonnegative integer, then 1. W P (n) is the number of Dyck paths of length n, specifically, W P (2n) = (2n)! n!(n + 1)! ; W P (2n + 1) = 0.

SW P (n)

is the number of Dyck paths without a peak of height 1; specifically

SW P (2n) = 1 2 n-2 i=0 (-1) i c n-i ( 1 2
) i ; SW P (2n + 1) = 0.

X(2n

) is the number of Dyck paths of length 2n without peaks at height 2 before the first return to height 0; specifically X(0) = 0 and ; W P (2n + 1) = 0.

X(2n) = F n-1 + F n for n > 0; X(2n + 1) = 2SW P (2n) = n-2 i=0 (-1) i c n-i ( 1
(2): Again, from the definition, it is also clear SW P is the set of Toppling peaks position that has a representation that is a Dyck path without peaks at height 1. That is SW P (2n) = F n . From [START_REF] Peart | Dyck paths with no peaks at height k[END_REF], Theorem 2, we have the recursive formulas

F n = n-1 i=1 c i * F n-i-1 and F n = 1 2 n-2 i=0 c n-i ( - 1 2 ) i .
(3): The number of Toppling peaks P-position of length 2n in the misère version is F n + F n-1 , where the second term counts the number of positions in SW P (2n -2) which have had an R appended at the beginning and an L at the end. The number of Toppling peaks P-position of length 2n + 1 in the misère version is 2 * F n .

Let E n be the number of Dyck paths of length 2n without peaks at height 2 before the first return to height 0. These can be counted by starting the sequence with an L then a member of SW P (2i) (a sequence without a peak of height 1) then an R then any member of W P (2n -2i -2) (that is a Dyck path of length 2n -2i -2). This gives the first line in:

E n = n-1 i=0 F i * c n-i-1 = F n-1 + n-1 i=1 c i * F n-i-1 = F n-1 + F n = X(2n)
3 timber on trees.

When timber is played on a tree, a move is to choose an arc (x, y) and remove from the tree the connected component containing y but not x. There may be an immediate move that finishes the game. If there isn't we present a reduction algorithm each part of which preserves the outcome class of the tree. One step removes 'peaks'. The second merges two 'peak'free paths incident with a common vertex.

The one move to finish the game is characterized first. Lemma 3.1 Let T be a tree, v a leaf of T and x the vertex adjacent to v. Let T t be a timber position which contains the arc (v, x), then o + (T t ) = N , that is T t is a next-player win position.

Proof. Alice wins by toppling the domino on the arc (v, x).

The next lemma shows that the two edges forming a 'peak' can be eliminated.

Lemma 3.2 Let T 1 t , T 2 t be two timber positions. Choose y ∈ V (T 1 t ), z ∈ V (T 2 t
) and let x be a vertex disjoint from T 1 and T 2 . Let T t be the position with vertex set V (T ) = V (T 1 ) ∪ {x} ∪ V (T 2 ) and arcs E(T t ) = E(T 1 t ) ∪ {(x, y), (x, z)} ∪ E(T 1 t ). Let T ′ t be the position with vertex set V (T ′ ) = V (T 1 ) ∪ V (T 2 ) where y and z are identified, and the arcs

E(T ′ ) = E(T 1 t ) ∪ E(T 2 t ). Then o + (T t ) = o + (T ′ t ).
Proof. We show it by induction on the number of vertices of T ′ . If V (T ′ ) = {y}, then there is no move in T ′ t and T t is the Toppling Peaks position LR. Hence o + (T ) = P = o + (T ′ ). Assume now |V (T ′ )| > 1. Assume Alice has a winning move in T t . It cannot be by choosing (x, y) or (x, z) because Bob would choose the other and win. If the chosen arc removes x from the game, choosing the same arc in T ′ leaves the same position. Otherwise, choosing the same arc in T ′ t leaves a position which has the same outcome by induction. Hence Alice has a winning move in T ′ t . The proof that Alice has a winning move in T t if she has one in T ′ t is similar.

If the tree has been reduced to two paths directed away from a common vertex x, then no move in either path affects the other. In either path any number from 1 to all the edges can be removed. This is the game of nim played with 2 heaps. Suppose the size of the paths/heaps are m and n. The Sprague-Grundy theory of Impartial games (see [START_REF] Albert | Lessons in Play[END_REF][START_REF] Berlekamp | Winning ways for your mathematical plays[END_REF] or numerous other books and papers) shows that the two heaps are equivalent to playing with one heap of size m ⊕ n where ⊕ is the XOR of integers. (Or write the numbers in binary and add without carrying. For example, 14 ⊕ 5 = 1110 2 ⊕ 101 2 = 1011 2 = 11.) The next result shows that even if there is more of the tree affixed at x, we can still make the replacement. Every impartial game is equivalent to a playing nim with single heap. The proof of the next result actually shows that when this merging of paths is done, the size of the equivalent nim heap doesn't change.

Lemma 3.3 Let T 0 be a tree, w ∈ V (T 0 ) a vertex, and n, m ∈ N two integers. Let T t be the position with vertex set V (T t ) = V (T ) ∪ {y i } 1 i n ∪ {z i } 1 i m and arcs

E(T t ) = E(T ) ∪ {(y i , y i+1 )} 1 i n-1 ∪{(z i , z i+1 )} 1 i m-1 ∪ {(w, y 1 ), (w, z 1 )}.
Let T ′ t be the position with vertex set

V (T ′ ) = V (T ) ∪ {x i } 1 i n⊕m
and the arcs

E(T ′ ) = E(T ) ∪ {(x i , x i+1 )} 1 i (n⊕m)-1 ∪ {(w, x 1 )}. Then o + (T t ) = o + (T ′ t ).
Proof. We prove it by induction on |V (T 0 )| + n + m and show that o + (T t + T ′ t ) = P which shows that o + (T t ) = o + (T ′ t ). If n + m = 0, T = T 0 = T ′ . Assume now |V (T 0 )| + n + m > 0. Any arc of T 0 is in both T t and T ′ t , thus if Alice chooses such an edge in one of T t or T ′ t then Bob can choose the corresponding arc in T ′ t or T t , which leaves a P position (either by induction or because the two remaining positions are the same). Assume Alice chooses the arc (y i , y i+1 ) (or (w, y 1 ) = (y 0 , y 1 )). If (i ⊕ m) < (n ⊕ m), Bob can choose the arc (x i⊕m , x (i⊕m)+1 ) (or (w, x 1 ) if i ⊕ m = 0) which leaves a P position by induction. Otherwise, there exists j < m such that (i ⊕ j = n ⊕ m), and Bob can choose the arc (z j , z j+1 ) which leaves a P position by induction. Similarly, we can prove that Bob has a winning answer to any move of the type (x i , x i+1 ) or (z i , z i+1 ).

A position which is neither of the form T for Lemma 3.2 or Lemma 3.3 is called minimal. A leaf-path is a path from some vertex x to a leaf y, x = y, consisting only of vertices, other than y and possibly x, of degree 2. Lemma 3.4 ∅ is the only minimal position which is P.

Proof. Let T t be a minimal position with at least one arc. If it has exactly one arc, it is obviously in N , so we can assume T t has at least two arcs. Then there exist a vertex w at which there are two leaf paths {x i } 1 i n and {y i } 1 i m . If (x n , x n-1 ) or (y m , y m-1 ) is an arc, Alice can choose it and win. Now assume both (x n-1 , x n ) and (y m-1 , y m ) are arcs. As T is minimal, it cannot be of the form from Lemma 3.2, so all (x i , x i+1 ), (y i , y i+1 ), (w, x 1 ) and (w, y 1 ) are arcs. But then it is of the form of Lemma 3.3, which is a contradiction.

The next algorithm shows how to apply the previous results to find the winning move, if there is one. Algorithm 3.5 Let T 0 be a tree and T 0 t a timber position.

1. Check each leaf x ∈ V (T 0 ); if (x, y) ∈ E(T 0 ) then toppling the domino on (x, y) wins, by Lemma 3.1, and o + (T 0 t ) = N and the algorithm stops.

If no such edge exists then let n = 0.

2. Flatten: Choose a vertex x ∈ V (T n ) at which there is a leaf path P to a leaf y and with at least one edge directed toward x. If there is u, v, w ∈ V (P ) such that (v, u), (v, w) ∈ E(T n t ) then T n+1 and T n+1 t are the tree and position obtained by applying Lemma 3.2 by removing v and identifying u and w. Let n := n + 1.

If there is a new edge, say (u, v) at the 'leaf ' end of the path and it is directed away from the leaf vertex, then by Lemma 3.1, and o + (T n t ) = o ( T 0 t ) = N . If no 'inward' arc is produced, eventually all the leaf-paths from x only have edges which are directed away from x.

Merge: If all vertices which have leaf-paths have all edges directed

toward the respective leaves then there is a vertex x which has two leafpaths say, P 1 , P 2 . Form T n+1 t by applying Lemma 3.3. Let n := n + 1 and repeat step 2.

The merging of paths in

Step 3 makes it awkward to describe how to find the winning move but with a little practice it becomes an easy matter for a player to discover it. The algorithm runs in O(n 2 ) time. Each edge is visited most once before at least one edge is deleted or a merging (Step 3) takes place. An edge on a merged path needs only be looked once, that in Step 2 when it will be immediately deleted.

As an example, we solve the position given in the Introduction. That the position is an N -position is shown by toppling the edge labelled '!' in the reduced position. This corresponds to toppling the edge labelled y. (After toppling y the merge of the two paths results in the empty path.) Note that the Flatten step can be applied at any appropriate degree 2 vertex at anytime. Given G and G t , suppose we regard the dominoes as physical. If a domino on (x, y) is toppled and there is another domino on (u, v) where the distance (in G) from y to u is the same as distance from y to v then this domino would have dominoes trying to topple it from both sides and so would remain standing! Call this game phyminoes. Then phyminoes is the same as timber on bipartite graphs. Problem 4.3 Is there an efficient algorithm to find the outcome class of a phyminoes on a non-bipartite graph?

Merge

Figure 1 :

 1 Figure 1: Example and Problem

Figure 2 :

 2 Figure 2: Three representations of one position.

Theorem 2 . 1 [

 21 Partition Theorem][1][page 41] Suppose the positions of a finite, impartial, normal-play game can be partitioned into mutually exclusive sets A and B with the properties: every option of a position in A is in B; and every position in B has at least one option in A. Then A is the set of P-positions and B is the set of N -positions. Theorem 2.2 [Misère Partition Theorem] Suppose the positions of a finite, impartial, misère-play game can be partitioned into mutually exclusive sets A and B with the properties: every option of a position in A is in B; and every position in B that has an option has at least one option in A. every position with no option is in B.

Figure 3 : 2 . 4 . 2

 3242 Figure 3: Example and Problem

  From the definition, it is clear that W P is the set of Toppling peaks positions that has a representation that is a Dyck path. Dyck paths only have even length. It is well known that the number of Dyck path of length 2n is the n th Catalan number c n =

		2	) i .
	Proof. (1): (2n)! n!(n+1)! , for example see [4]
	p.194. Thus
	W P (2n) =	(2n)! n!(n + 1)!