Nalin Srun

Sotheara Leang

Ye Kyaw

Sethserey Sam

Convolutional Time Delay Neural Network for Khmer Automatic Speech Recognition

Keywords: Khmer ASR, Time Delay Neural Network, Convolutional Neural Network, Lowresource Language

Convolutional Neural Networks have been proven to successfully capture spatial aspects of the speech signal and eliminate spectral variations across speakers for Automatic Speech Recognition. In this study, we investigate the Convolutional Neural Network with Time Delay Neural Network for an acoustic model to deal with large vocabulary continuous speech recognition for Khmer. Our idea is to use Convolutional Neural Networks to extract local features of the speech signal, whereas Time Delay Neural Networks capture long temporal correlations between acoustic events. The experimental results show that the suggested network outperforms the Time Delay Neural Network and achieves an average relative improvement of 14% across test sets.

Introduction

Google Assistant, Siri, Alexa, and Cortana are transforming the way users interact with their devices, homes, and automobiles. These voice assistants simplify our lives, and many individuals become dependent on them. A key aspect of that smart merchandise is Automatic Speech Recognition (ASR). It permits you to talk into your computer or tool, and it translates your voice into textual information automatically. At present, a variety of approaches including conventional models [START_REF] Lawrence R Rabiner | A tutorial on hidden markov models and selected applications in speech recognition[END_REF][START_REF] Stolcke | Mllr transforms as features in speaker recognition[END_REF], hybrid models [START_REF] Hinton | A fast learning algorithm for deep belief nets[END_REF][START_REF] Mohamed | Deep belief networks for phone recognition[END_REF] and end-to-end neural models [START_REF] Graves | Speech recognition with deep recurrent neural networks[END_REF] are carried out in almost ASR systems to achieve the state of the art speech recognition accuracy. The rapid advancement of ASR is heavily dependent on the massive amounts of audio and annotated transcripts. Recognizing the majority of languages having a rich resource gives a very mature performance [START_REF] Stenman | Automatic speech recognition an evaluation of google speech[END_REF] when compared to low-resource languages, which have so far achieved limited accuracy [START_REF] Besacier | Asr and translation for under-resourced languages[END_REF].

In a recent study [START_REF] Soky | Khmer automatic speech recognition system based on dnn models[END_REF], Khmer ASR was studied with Deep Neural Network (DNN) versus Gaussian Mixture Model (GMM), and the result showed that the DNN outperformed by approximately 3.65% and 1.10% in the open test and close test, respectively. However, when working in a language with limited resources, DNN results in overfitting [START_REF] Li | Dnn online adaptation for automatic speech recognition[END_REF]. Furthermore, the DNN can only acquire a limited amount of contextual information from the speech, limiting its ability to deal with long-range correlations in the speech signal [START_REF] Li | Lstm time and frequency recurrence for automatic speech recognition[END_REF][START_REF] Yu | Recent progresses in deep learning based acoustic models[END_REF]. In recent years, various novel networks have been studied, and it has been discovered that Time Delay Neural Network (TDNN) is one of the potential neural networks for modeling the long temporal information of the speech [START_REF] Peddinti | A time delay neural network architecture for efficient modeling of long temporal contexts[END_REF][START_REF] Povey | Semi-orthogonal low-rank matrix factorization for deep neural networks[END_REF]. The TDNN is a stacked architecture that employs a modular and progressive design to build a bigger network from the sub-components. In the lower layers, it learns the narrow context of the input features, whereas in the top layers, it learns the larger context. Furthermore, it is quicker to train TDNN and requires much less training data compared to Recurrent Neural Network (RNN) and Long Short-term Memory (LSTM) [START_REF] Shewalkar | Performance evaluation of deep neural networks applied to speech recognition: Rnn, lstm and gru[END_REF]. Speech can be regarded as an image because the spectral properties of speech consist of two dimensions, such as time and frequency. As a result, several researches on Convolutional Neural Network (CNN) has been conducted for speech recognition [START_REF] Musaev | Image approach to speech recognition on cnn[END_REF]. The three critical components of CNN, including local receptive field, weight sharing, and sub-sampling, are known to be capable of minimizing the variance of translation, scaling, and distortion of the speech signal [START_REF] Guiming | Speech recognition based on convolutional neural networks[END_REF]. Therefore, CNN has additional advantages for speech processing when dealing with diverse pronunciations for different speakers. In this research, we propose to investigate the benefits of combining CNN and TDNN in dealing with large vocabulary continuous speech recognition for Khmer. Our goal is to utilize the CNN as an extra function to the TDNN to learn local features of the speech signal. This research is structured as follows: the proposed method is given in Section 2. Section 3 introduces the experiments including the datasets, language model, lexicon, and baseline model. The results and discussions are presented in Section 4. Finally, Section 5 presents the conclusion and future work.

Proposed Method

In this investigation, the CNN-TDNN of Kaldi's Mini-LibriSpeech recipe [START_REF]Kaldi mini-librispeech recipe[END_REF] was adopted. The network consisted of 6 convolutional neural network blocks, followed by 9 delay neural network blocks, and a fully connected block (Fig. 1). It took the 40-dimensional Mel-frequency Cepstral Coefficients (MFCCs) and 100-dimensional i-vector as the input features to estimate 2,976 states of the triphone Hidden Markov Model (HMM). The network contains 4.7 million parameters and can gather ±30 contextual inputs from a particular acoustic frame.

Feature Extraction

The speech signals were first split into short-time frames of a length of 25ms with an overlap of 10ms. A Hamming window, a pre-emphasis factor of 0.97, a cepstral liftering coefficient of 22 and cepstral mean normalization (CMN) were used. Finally, 40dimensional high resolution MFCCs were computed and 100-dimensional i-vectors were extracted on top of PCA-reduced spliced-MFCC features for speaker adaptation.

Preprocessing Block

The input features were preprocessed before being fed to the network as shown in Fig. 2. First, mel-filterbanks were computed from 40-dimensional MFCCs through inverse discrete cosine function (MFCCs are more compressible, thus it is preferred to dump them to disk rather than mel-filterbanks). Second, batch normalization and SpecAugment were applied on the mel-filterbanks. In Kaldi, the SpecAugment is implemented as a neural layer that does time and frequency masking on-the-fly during the training. Third, 100-dimensional i-vectors were projected to 200-dimensional vectors through a linear transformation, followed by a batch normalization. Finally, 40-dimensional mel-filterbanks and 200-dimensional i-vectors were combined to produce 40x6 input features.

Convolutional Neural Network Block

Our network employed 6 CNN blocks, each with a convolutional layer followed by a rectified linear unit (ReLU) activation function and a batch normalization, to learn local spatial information and to eliminate variations in inter-speaker spectra (Fig. 3). Each convolutional layer took three consecutive outputs from the preceding layer of time steps (t-1, t, t+1), allowing the CNN blocks to gather information from 12 neighboring frames. During each convolution operation, 3x3 filters with a strip size of 1 and a zeropadding size of 1 were employed. The first three convolutional layers use 48 filters; the next two lay-ers use 64 filters; and the last layer uses 128 filters. Layers 3, 5, and 6 use a subsampling size of 2 to lower the height dimension of the input feature. The final outputs (1x5x128) were transformed into 640dimensional vectors and utilized as inputs to TDNN sequence blocks.

Time Delay Neural Network Block

Our network was made up of 9 TDNN blocks, and thus, it can extract information from the long temporal context of the inputs. Each block begins with a TDNN layer, followed by a ReLU activation function and a batch normalization function (Fig. 4). The Factored Time Delay Neural Network (TDNNF) was used in the network. It uses Single Value Decomposition (SVD) to improve matrix multiplication, resulting in fewer parameters and much less computation. The first TDNNF layer has a dimension of 768 while SVD has a dimension of 192. Time stride was not applied in this layer. The subsequent TDNNF layers consist of 768 dimensions and 96 dimensions for SVD. Each layer took 6 consecutive outputs from ±3 contexts and concatenated them with the current output from the previous layer. Therefore, the final output from the TDNN blocks can capture information from the previous and the prior 24 inputs. RestNet connections were applied in TDNN layers to capture information from the input from the previous layer. Except for the first layer, the output from each layer was summed with the output from the previous layer with a scaling factor of 0.66.

Training Recipe

A context-independent (CI) model was first trained on 13-dimensional MFCCs with their first and second-order derivatives. MFCCs were computed using 23 Mel-filterbanks on overlapping frames of 25ms with 10ms shift. A Hamming window, a pre-emphasis factor of 0.97, a cepstral liftering coefficient of 22 and CMN were used. The first context-dependent (CD) Figure 4: Time delay neural network block. The number n denotes the size of adjacent inputs to be combined with the current input t.

was trained over 40 iterations using 2.5k states and 25k Gaussian distributions. Second CD training using 5k states and 50k Gaussians was performed by concatenating ±4 contexts of 13-dimensional MFCCs and projecting into a 40-dimensional subspace using linear discriminant analysis (LDA), and the features were decorrelated using maximum likelihood linear transform (MLLT). Speaker-adaptive training (SAT) was performed with the same configuration of 5k states and 50k Gaussians to normalize the features using feature space maximum likelihood linear regression (fMLLR). Note that the configurations of CD models are based on [START_REF] Soky | Khmer automatic speech recognition system based on dnn models[END_REF].

Three data augmentations were used to expand the amount of the training data and provide greater signal variability in speech. First, speed perturbation with factors of 0.9, 1.0, and 1.1 was carried out, producing 3 times the training data. Second, volume perturbation was used to make the acoustic model resilient to volume fluctuations by using random factors sampled from a uniform distribution of [1/8, 2]. Finally, SpecAugment [START_REF] Daniel S Park | Specaugment: A simple data augmentation method for automatic speech recognition[END_REF] was conducted on the fly in the neural layer of the network by masking the frequency channels of at most 50% of the spectrum and masking the time steps of the frames of at most 20%.

The 40-dimensional high-resolution MFCCs were extracted from the augmented training data, and half of the data was used to train the i-vector extractor [START_REF]Kaldi ivector recipe[END_REF], which was used to compute the 100-dimensional i-vectors from speech signals. Finally, CNN-TDNN was trained on the aligned training data from the SAT model across 20 epochs with starting and final learning rates of 0.002 and 0.0002, respectively, using lattice-free maximum mutual information (LF-MMI) [START_REF] Povey | Purely sequence-trained neural networks for asr based on lattice-free mmi[END_REF].

Datasets

In this study, the training set includes portions of the Basic Travel Expression Corpus (BTEC) [START_REF] Kikui | Creating corpora for speech-to-speech translation[END_REF] consisting of 47 speakers and 48,369 utterances. It is a speech corpus that was produced from the previous work [START_REF] Soky | Khmer automatic speech recognition system based on dnn models[END_REF] and designed to cover the expressions in the traveling domain. More data was collected from the internet (online websites and YouTube) and combined into the datasets to examine the performance of Khmer ASR in several domains. Finally, the new speech corpus, which is mostly in natural speech and background noise, and contains 119 speakers, 73,660 utterances, and around 153 hours, was used during the learning phase. Furthermore, we proposed to investigate the performance of CNN-TDNN in five domains, including News, Law, Health, Daily life and General. The specification of the training and testing sets are given in Table 1.

Language Model

The language model is an essential component of the ASR system. By estimating the likelihood of word sequences, it is able to model how words are put together to form sentences. It can predict which words will follow particular words and with what probability. A 3-grams language model was built using SRILM [START_REF]The sri language modeling toolkit[END_REF]. A Witten-Bell Smoothing technique was used to prevent assigning zero probability to the unknown words and a pruning factor of 1e-8 was applied to reduce the size of the language model. Using an in-house toolkit [START_REF] Chea | Khmer word segmentation using conditional random fields[END_REF], the text of each dataset was segmented into morphemes and subjected to punctuation removal and numbers to text conversion. Table 2 shows specification of the training and testing datasets as well as their perplexities.

Lexicon

The lexicon is the pronunciation dictionary that connects sub-word units (phonemes) to words. The acoustic model uses it to learn mapping from sequence of phonemes to words. The corpus from [START_REF] Soky | Building wfst based grapheme to phoneme conversion for khmer[END_REF] was used to build a grapheme-to-phoneme (G2P) model using Phonetisaurus G2P [START_REF]Phonetisaurus g2p[END_REF]. The dataset, which included 57 phonemes and 34k words, was created using the International Phonetic Alphabet (IPA) standard. Finally, a new lexicon was constructed using the vocabulary of the language model, which consists of around 100k words.

Baseline Model

In this study, we proposed using TDNN as the baseline model. The suggested network design is given in Fig. 5. It is made up of 15 TDNN blocks and is capable of capturing 30 adjacent contexts from a partic-ular input feature. The network is similar to CNN-TDNN except for the first six TDNN blocks. Our idea is to investigate the performance of the acoustic model when CNNs are utilized in the network while maintaining the same configuration. Each TDNN accepts 3 successive input features (t-1, t, t+1), and thus, the TDNN blocks can obtain information from 6 neighboring input features.

Evaluation Method

Word Error Rate (WER) was used as the metric to evaluate the accuracy of TDNN and CNN-TDNN. It is a common measure for evaluating the performance of an ASR system. The following is the formula:

W ER = S + D + I S + D + C ( 1
)
Where S is the number of substitutions, D is the number of deletions, I is the number of insertions, and C is the number of correct words. WER is between zero and one. The lower the WER, the better performance of the acoustic model. In addition, error analysis was conducted to examine the inaccuracy of the model in generating the hypothesis by using SCLITE [START_REF]Sctk, the nist scoring toolkit[END_REF]. We discovered that the majority of the problems were due to phonetic and encoding issues. A phonetic mistake occurs when the model estimates a word that sounds the same as another word but has a different meaning (homophones). When the model generates the same word but with a different encoding, this is referred to as an encoding error. We contend that these issues are caused by the quality of the corpus used to build the lexicon and the language model. Tables 4 and5 The first example of phonetic errors, the demonstrative adjective "ហ្ន ឹ ង" ("this" in English), was predicted as the conjunction "ន ិ ង" ("and" in English) when they In encoding error instances, coeng "ដ" (178A) was predicted as coeng "ត" (178F). As a result, the word "ឧត្ដ ម", and "សេម្ដ ច" were synthesized as "ឧ + ត + ◌្ + ត + ម" and "ស + ម + ◌្ + ត + េ◌ + ច" instead of "ឧ + ត + ◌្ + ដ + ម" and "ស + ម + ◌្ + ដ + េ◌ + ច", respectively.

Results and Discussions

Conclusion

In this study, we reported the state-of-the-art for dealing with Khmer ASR using CNN and TDNN.

The results reveal that integrating the CNNs into the network of TDNN significantly improves the performance of the model across all test sets. In future work, a new acoustic modeling approach will be investigated to improve the overall performance of the model in terms of accuracy and decoding speed. Furthermore, due to the constrained variety of datasets in comparison to different well-resourced languages, a larger dataset is desired for future development of Khmer ASR. More data needs to be acquired and used for various data augmentation methods in order to examine more environmental factors. Khmer encoding issues and discrepancies will also be studied to eliminate encoding errors while building a corpus, and it is widely assumed that different misspelling errors are no longer a problem, as long as there is a huge amount of correct data for speech and text corpora.

Figure 1 :

 1 Figure 1: Adopted architecture for CNN-TDNN

Figure 2 :

 2 Figure 2: Preprocessing block

Figure 3 :

 3 Figure 3: Convolutional neural network block

Figure 5 :

 5 Figure 5: Proposed architecture for baseline model

  show the most common examples of phonetic and encoding errors, respectively. Phonetic error examples: Scores: (#C #S #D #I) 20 1 0 0 REF: េដាយសារ ែត មុ ខ របរ ហ្ន ឹ ង មាន ការ រ ី ក ចេ្រម ី ន ខា ្ល ំ ង នាង ខ្ញ ុ ំ ក៏ ព្រងី ក សាខា េទ លក់ េន ក្ន ុ ង ខណ ្ឌ ែសន (Due to the rapid growth of this business, I expanded the branch to sell in Sen district) HYP: េដាយសារ ែត មុ ខ របរ ន ិ ង មាន ការ រ ី ក ចេ្រម ី ន ខា ្ល ំ ង នាង ខ្ញ ុ ំ ក៏ ព្រងី ក សាខា េទ លក់ េន ក្ន ុ ង ខណ ្ឌ ែសន (Due to the business and rapid growth, I expanded the branch to sell in Sen district) Scores: (#C #S #D #I) 19 1 0 0 REF: កាត់ បន្ថ យ នូ វ ការ ជួ ប ជុ ំ គា ្ន ែដល បច្ច ុ ប្បន្ន េយ ី ង កំ ពុ ង ែត ជួ ប នូ វ បញ ្ហ ជម្ង ឺ ឆ្ល ង កូ វ ីដ ដប់ ្របាំ បួ ន (Reducing the gatherings when we are currently facing with Covid 19) HYP: កាត់ បន្ថ យ នូ វ ការ ជួ ប ជុ ំ គា ្ន ែដរ បច្ច ុ ប្បន្ន េយ ី ង កំ ពុ ង ែត ជួ ប នូ វ បញ ្ហ ជម្ង ឺ ឆ្ល ង កូ វ ីដ ដប់ ្របាំ បួ ន (Reducing the gatherings as well when we are currently facing with Covid 19)

Table 3

 3

compares the performance of CNN-TDNN with TDNN across the five test sets. TDNN performs

Table 1 :

 1 Dataset Preparation

	Dataset	Domain	Hour Utterance	Speaker
	Train Set Multi-domain ∼ 153	73,660	57 males, 62 females
		News	∼ 4	970	2 males, 3 females
		Law	∼ 3	1,035	2 males, 3 females
	Test Set	Health	∼ 4	1,399	2 males, 3 females
		Daily Life	∼ 5	1,840	2 males, 3 females
		General	∼ 4	1,474	2 males, 3 females

Table 2 :

 2 Dataset Preparation for Language Model

	Dataset	Domain	Vocabulary Sentence Average Word Perplexity
	Train Set Multi-domain	∼ 100K	∼ 3M	53	036.42
		News	∼ 3K	∼ 3K	28	090.38
		Law	∼ 2K	∼ 3K	20	158.29
	Test Set	Health	∼ 3K	∼ 3K	20	158.29
		Daily Life	∼ 3K	∼ 4K	24	238.94
		General	∼ 3K	∼ 3K	22	169.29
	best in the News, with a WER of 11.53%, followed by
	17.57% in the Health and Daily Life, and a WER of
	roughly 19% in the Law and General. CNN-TDNN,
	on the other hand, outperforms by a significant mar-
	gin of 14.4% on average throughout the test sets. It
	achieves the best WER of 11.53% in the News, fol-
	lowed by 14.83% in the Daily Life and around 15.7%
	in the Health and General. According to the exper-
	imental results, introducing CNNs into the network
	significantly improves performance when compared
	to TDNN with the same configuration and capability
	in capturing long temporal dependency of 30 adja-
	cent contexts. Thus, CNNs enable the network to
	efficiently learn representation of the speech signal.

Table 3 :

 3 Experimental results in WER(%) on each test set

	Model	News	Law	Health Daily Life General
	TDNN	11.53	19.18	17.57	17.57	19.17
	CNN-TDNN	09.94 (↓13%)	16.02 (↓16%)	15.70 (↓10%)	14.83 (↓15%)	15.68 (↓18%)
	Table 4: Phonetic Error
	REF ⇒ HYP Freq
	ហ្ន ឹ ង ⇒ ន ិ ង អ្ន ក⇒នាក់	233 38
	ែដល⇒ែដរ	28
	are both homophones. Similarly, the particle "ែដល" of adverb "បច្ច ុ ប្បន្ន " in second example, was predicted as adverb "ែដរ" and these two words have the same pronunciation.

Table 5 :

 5 Encoding Error

	REF ⇒ HYP Freq
	សេម្ដ ច⇒សេម្ត ច ឧត្ដ ម⇒ឧត្ត ម	32 30
	Encoding error examples:
	Scores: (#C #S #D #I) 13 1 0 0
	REF: ការ សេ�ងា ្គ ះ ជាត ិ ជា ្របវត្ត ិ សា�ស្ត របស់ សេម្ដ ច អគ្គ មហា េសនា បត ី េតេជា ហ៊ ុ ន ែសន (National salvation is the history of Samdech Akka Moha Sena Padei Techo
	Hun Sen) (សេម្ដ ច ⇒ 179F 1798 17D2 178A 17C1 1785 ) HYP: ការ សេ�ងា ្គ ះ ជាត ិ ជា ្របវត្ត ិ សា�ស្ត របស់ សេម្ត ច អគ្គ មហា េសនា បត ី េតេជា ហ៊ ុ ន ែសន (សេម្ត ច ⇒ 179F 1798 17D2 178F 17C1 1785)
	Scores: (#C #S #D #I) 9 1 0 0
	REF: រដ្ឋ ម�ន្ត ី ្រកសួ ង ការពារ ជាត ិ រួ ម ជាមួ យ ន ឹ ង ឯក ឧត្ដ ម (Minister of National Defense together with His
	Excellency) (ឧត្ត ម ⇒ 17A7 178F 17D2 178A 1798) HYP: រដ្ឋ ម�ន្ត ី ្រកសួ ង ការពារ ជាត ិ រួ ម ជាមួ យ ន ឹ ង ឯក ឧត្ត ម (ឧត្ត ម ⇒ 17A7 178F 17D2 178F 1798)

Acknowledgment

We would like to express our gratitude to CADT students, particularly Mr. Rithborey Chorn and Mr. Samnang Pov, for their support in collecting and preparing the speech corpus.


P Roldan

S Trilles

X Serena

A Tajdine

Xavier Serena Alos

Novel Composite Clock Algorithm for the Generation of Galileo Robust Timescale

Pedro Roldan works on time and frequency systems and their application to GNSS. He holds a master degree in Aerospace Engineering and in Meteorology and Geophysics, and is currently following PhD studies. He has participated in the design, development, monitoring and operation of the Galileo system, focused on the application of algorithms of time transfer, generation of timescales and steering. Currently, he is in charge of time and frequency activities in Thales Alenia Space.

Sébastien Trilles is an expert in navigation algorithms and performances. He received his Ph.D. degree in Pure Mathematics from the Paul Sabatier University and an Advanced M.S. in Space Technology from SUPAERO. He heads the Performance and Processing department where high precise algorithms are designed as the navigation, integrity and ionosphere modelling algorithms of SBAS systems, the orbitography solutions and the composite time reference generation for Galileo.

INTRODUCTION

Both master clock and composite clock approaches are currently used by GNSS systems, being both able to comply with the very demanding stability and availability levels required by GNSS applications. In the master clock concept, this is granted by the use of very stable master clocks based on hydrogen technologies and the definition of the necessary redundancies to cope with anomalies in the clocks and maintenance and operational activities, while in the paper time concept it is ensured by using composite clock algorithms capable of detecting anomalies in the contributing clocks and adapting the list of participants and their characterization when anomalies are detected. Algorithms of composite clock allow to combine different clocks to obtain a timescale more stable than any of the input clocks and insensitive to anomalies and gaps in them.

Most algorithms of composite clocks are based on estimating the state of the clocks with a Kalman filter, considering the noises associated to the measurements and the clock states. The time given by the individual clocks is not directly accessible for the composite clock generation, but only the time difference between pairs of clocks. The use of time differences as input measurements to the filter implies that the number of measurements is always smaller than the number of states. If N clocks are considered, at least N different states should be estimated with only N-1 independent measurements, obtained as differences -݅ܭܮܥ .݆ܭܮܥ In case this observability problem is not addressed, the covariance associated to these N states would indefinitely increase and a drift could appear in the estimated states. This can be solved by limiting the covariance to the part associated to the independent measurements, as done in the algorithms of [START_REF] Brown | The Theory of the GPS Composite Clocks[END_REF] and [START_REF] Greenhall | A Kalman filter clock ensemble algorithm that admits measurement noise[END_REF][START_REF] Greenhall | A Kalman filter clock ensemble algorithm that admits measurement noise: corrections and update[END_REF], or by introducing constraints as additional observations, making the number of observations equal to the number of states, as done by the algorithms of IGS, version 1 [START_REF] Senior | Developing an IGS Time Scale[END_REF] and version 2 [START_REF] Senior | The Next Generation GPS Time, NAVIGATION[END_REF].

The differences of clocks used as measurements can be considered in phase, as done by the algorithms of Brown, Greenhall and IGS version 2, or in frequency, as done by IGS version 1. For modelling the clocks, different states can be considered. Generally, the three states of phase, frequency and drift are considered when processing measurements in phase, and two states of frequency and drift are considered when processing measurements in frequency. For the case of IGSv2, an additional state for the total phase, including a white noise, is also included.

In a general problem, the measurements would be those of Equation ( 1), being phase or frequency depending on the algorithm, while the estimated states would be those of Equation ( 2), with some of the states potentially set to zero and with the possibility of adding other states depending on the modelling of the clocks.

ܼ = ܭܮܥ൫ -ܭܮܥ ൯ / (1) ܺ = ܭܮܥ( -)ܥܥ ܭܮܥ( -)ܥܥ ܭܮܥ( -)ܥܥ ௗ (2)
This implies that certain states are directly observable with the measurements, while some others only appear when propagating the observable states. Focusing on the case with processing of phase measurements and estimation of three states, the measurements are directly related to the phase states, while the frequency and drift states are indirectly observed through the propagation step of the Kalman filter. The presence of frequency and drift states increases the stability of the solution in the long-term, since this indirect observability makes the estimation of these states varying more slowly than the estimation of the phase states. Despite this benefit related to the increase of stability, the slower estimation of frequency and drift states may generate an anomalous behavior in the composite timescale in case the frequency and drift of some of the input clocks suddenly change, since it would not be able to properly reflect this change in the estimated states.

This limitation when facing sudden changes in the clocks has been identified and analyzed by [START_REF] Galleani | Time and the Kalman Filter[END_REF], together with other limitations of a Kalman-based estimation of clocks, like the behavior in case of addition and removal of clocks. The estimation of not directly observed states may also produce unexpected behaviors during the convergence of the Kalman filter. As shown by Wright (2007a[START_REF] Wright | Composite clocks with 3-state models[END_REF], accumulated errors in the estimation of frequency and drift states during the convergence can lead to an arbitrary offset in frequency and drift in the resulting timescale.

For the application as a GNSS system time, the robustness against operational events shall be the main design driver of a composite clock algorithm. For this reason, the algorithm described in this paper has been designed to mitigate as much as possible these limitations. The non-observable states have been removed with the addition of constraints linking the estimated states to the predicted ones and allowing a change from epoch to epoch in the level of the noise of each particular clock, following the approach of IGS version 2 [START_REF] Senior | The Next Generation GPS Time, NAVIGATION[END_REF]. The indirectly observable states have been limited to those cases in which they provide a clear benefit for the modelling of the clocks. To do that, the possibility of adding steered clocks contributing to the composite has been considered, being the state vector for these clocks only composed by phase states, directly observed with the clock differences used as measurements in the filter. To avoid the impact of reconfigurations of the ensemble, all the input clocks are aligned before being injected in the filter, so that the phase, frequency and drift offset between them is never too large to avoid a smooth transition in the composite timescale. To avoid any dependency with a particular clock, the approach of injecting differences with respect to a single reference clock used in the algorithm of IGS has been here replaced by the injection of a set of differences without any privileged clock. Among the available clock differences, most of them would depend on the others and would not really add new information to the filter. A set of independent differences could be selected instead, according to particular criteria like the residuals of the estimation.

ROBUST COMPOSITE CLOCK

The composite clock algorithm considered in this paper contains the steps of prediction and update of the Kalman filter, together with additional functions for preprocessing of measurements and anomaly detection based on the residuals of the estimation. A high level architecture of this algorithm is included in Figure 1.

Figure 1 Architecture of robust composite clock

The preprocessing function analyses the input differences to detect and correct clock anomalies and operational events like outliers, phase and frequency jumps. To do that, it compares the clock differences for one epoch with the clock differences predicted from the previous epochs. In case an event is detected, two strategies would be possible: rejecting the clock differences or correcting the event. If the accuracy is the main driver, the strategy of measurement rejection would be very likely the best one, since the correction of events like phase or frequency jumps could be not accurate enough to allow for an accurate estimation of the clocks in the filter. However, this strategy could degrade the availability, especially in scenarios with a large amount of events. If the availability is the main driver or there are too frequent events, the option of correcting the events and include the corrected measurements in the filter could provide better results.

Even if the algorithm considers constraints to limit the change in the composite timescale from one epoch to the next one, in case very large phase, frequency or drift offsets exists among the different clocks, a significant discontinuity may appear in case clocks are added or removed from the ensemble. To mitigate this problem, an alignment of the input clocks is also performed in the preprocessing step, being its goal not to accurately model the input clocks, but only to ensure that the differences between them do not exceed what can be absorbed by the Kalman filter constraints.

In the prediction step, the state vector of Equation ( 2) and the associated covariance matrix ܲ are propagated a time step ߬according to the Equation (3). In this equation, the transition matrix ܨ is obtained by composing the transition matrices of each clock ܨ , as per Equation ( 4), and the process noise matrix ܳ is generated with the matrices of each clock ܳ , following Equation ( 5).

ܺ = ܺܨ ܲ = ܨܲܨ ᇱ + ܳ (3) ܨ = 1 ߬ ߬ ଶ 2 ⁄ 0 1 ߬ 0 0 1 (4) ܳ = ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎡ 1ݍ • ߬+ 2ݍ • ߬ ଷ 3 + 3ݍ • ߬ ହ 20 2ݍ • ߬ ଶ 2 + 3ݍ • ߬ ସ 8 3ݍ • ߬ ଷ 6 2ݍ • ߬ ଶ 2 + 3ݍ • ߬ ସ 8 2ݍ • ߬+ 3ݍ • ߬ ଷ 3 3ݍ • ߬ ଶ 2 3ݍ • ߬ ଷ 6 3ݍ • ߬ ଶ 2 3ݍ • ߬ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎤ (5)
The noises 1ݍ , 2ݍ and 3ݍ respectively represent the contributions of Random Walk PM, Random Walk FM and Random Walk DM to the clock noise. These noises can be obtained based on the Hadamard variance of the clocks, following the approach in [START_REF] Hutsell | Relating the Hadamard Variance to MCS Kalman Filter Clock Estimation[END_REF].

With the predicted state, the pre-fit residuals ܻ can be computed as in Equation ( 6). In this equation, the measurements ܼ are compared against the reconstructed measurements, obtained by applying the design matrix ܪ to the state vector ܺ. To solve the observability problem derived from the estimation of N clocks with N-1 clock differences, the estimated state is constrained to the predicted one, following the approach in [START_REF] Senior | The Next Generation GPS Time, NAVIGATION[END_REF]. Taking this into account, the vector of measurements ܼ is composed by the clock differences themselves ܼ and additional rows for the constraints ܼ ௦௧ , as per Equation ( 7). In the same way, the design matrix ܪ contains one row ܪ for each clock difference included in the filter and additional rows ܪ ௦௧ for the constraints, as shown in Equation ( 8).

ܻ = ܼ -ܺܪ (6) ܼ = ܭܮܥ൫ -ܭܮܥ ൯ ܼ ௦௧ = ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎡ ݓ , • ܭܮܥ( -)ܥܥ ,ௗ ୀ ݓ , • ܭܮܥ( -)ܥܥ ,ௗ ୀ ݓ ,ௗ • ܭܮܥ( -)ܥܥ ௗ,ௗ ୀ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎤ (7) ܪ = [… 1 0 0 … -1 0 0 …] ܪ ௦௧ = ݓ ଵ, 0 0 … ݓ ே , 0 0 0 ݓ ଵ, 0 … 0 ݓ ே , 0 0 0 ݓ ଵ, ௗ … 0 0 ݓ ே ,ௗ (8)
These constraints allow to define the evolution of the composite timescale, by constraining a weighted average of the estimated states ܭܮܥ( -)ܥܥ , ܭܮܥ( -)ܥܥ and ܭܮܥ( -)ܥܥ ௗ to a weighted average of the predicted states ܭܮܥ( -)ܥܥ ,ௗ , ܭܮܥ( -)ܥܥ ,ௗ and ܭܮܥ( -)ܥܥ ௗ,ௗ . Those clocks with a larger weight in the constraint will have a larger contribution to the evolution of the composite timescale. This approach, even if strongly dependent on the selection of weights, allows to easily remove the contribution of a particular clock to the composite timescale by reducing its weight. It is therefore a suitable approach for an operational application like a GNSS system time, in which it shall be possible to totally remove the contribution of clocks impacted by anomalies or involved in operational and maintenance activities.

The residuals are computed for all the available clock differencesܭܮܥ -ܭܮܥ , but the processing of all these clock differences in the composite clock would not be useful to improve the estimation of the clock states, since they provide redundant information. For example, if the measurements ܼ , ܼ and ܼ are included in the Kalman filter, only two of them would efficiently contribute to the estimation, since the third one could be obtained from the others as ܼ = ܼ -ܼ , potentially with a different measurement noise but with the same link to the estimated parameters. The residuals of the clock differences could be used then as a criterion to select those differences that better allow to characterize the behavior of the clocks, and only include these differences in the update step of the Kalman filter.

In the update step, the Kalman gain ܭ is obtained based on the predicted covariance ܲ and the measurement noise ܴ. The residuals ܻ of a set of independent measurements are used then to update the state vector ܺ and the associated covariance, according to Equation ( 9).

ܵ = ܪܲܪ ᇱ + ܴ ܭ = ܪܲ ᇱ ܵ ିଵ ܺ = ܺ + ܻܭ ܲ = -ܫ( -ܫ(ܲ)ܪܭ )ܪܭ ᇱ + ′ܭܴܭ (9)
The measurement noise matrix ܴ can be obtained by composing the measurement noise of each measurement ܴ and the measurement noises associated to the constraints ܴ ௦௧ . Since the constraints are based on a weighted average of the predicted states, a weighted average of the associated covariance could be a good choice, as described in Equation ( 10).

ܴ = 0ݍ + 0ݍ ܴ ௦௧ = ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎡ ݓ , • ܲ , ୀ 0 0 0 ݓ , • ܲ , ୀ 0 0 0 ݓ ,ௗ • ܲ ,ௗ ୀ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎤ (10)
The noise of the measurements ܴ could be modelled as a white noise 0ݍ and 0ݍ , even if additional contributions may be needed, depending on the source of the measurements, to take into account the noise introduced by the measuring system.

The residuals can be also used as a second barrier to detect anomalies that have not been detected in the preprocessing. For this purpose, the residual of the measurement ܻ is normalized with the covariance of the clocks ܲ , and ܲ , and the measurement noise 0ݍ + 0ݍ to obtain the normalized residual ܰ , as per Equation ( 11).

ܰ = ܻ ඥ ܲ , + ܲ , + 0ݍ + 0ݍ (11)
Those clock differences showing a large normalized residual are rejected, and in case several differences are rejected for the same clock, the clock is excluded from the filter. In a first step, the clock is de-weighted, by progressively reducing its weight in the constraints of Equations ( 7), ( 8) and ( 10). If the clock is not recovered, it is de-correlated from the other clocks by removing the crossed terms of the covariance matrix and using the clock differences involving that clock only to estimate its own state. Like this, the estimation of the state of the other clocks is not contaminated by the differences between them and the anomalous clock. If the situation is not automatically recovered, the state vector and the covariance matrix for the anomalous clock are finally reinitialized, allowing for a new estimation from scratch.

Apart from the presence of anomalies, other factors like the noise and the availability of data also contribute to define the weight of the clock, as shown in Equation ( 12). The weights are obtained as the inverse of the process noise for each state, combined with multiplicative factors to reduce the weight in case of data gaps or anomalies.

ݓ , = ݅ܽݒܣ ݈ • ݉݊ܣ 1ݍ • ߬+ 2ݍ • ߬ ଷ 3 + 3ݍ • ߬ ହ 20 ݓ , = ݅ܽݒܣ ݈ • ݉݊ܣ 2ݍ • ߬+ 3ݍ • ߬ ଷ 3 ݓ ,ௗ = ݅ܽݒܣ ݈ • ݉݊ܣ 3ݍ • ߬ (12)
The factors ݅ܽݒܣ ݈ and ݉݊ܣ would be 1 in case the clock has a nominal behaviour, and would smoothly reduce in case no more data is available for a particular clock, or in case an anomaly is detected when analysing its residuals.

PROPOSED INTEGRATION IN GALILEO

The use of a composite clock as a GNSS system time can be put in place by setting it as the reference time of the Orbit Determination and Time Synchronization (ODTS) process that generates the predicted orbits and clocks included in the navigation message of the system. To do that, the composite clock could be integrated with the ODTS, in a tightly coupled approach, or could be done separately from the ODTS, in a loosely coupled approach. The tightly coupled approach would be optimized in terms of performance, since both the composite clock and the ODTS are included in the same estimator. However, this solution may be not optimal from a robustness point of view. The loosely couple approach would allow for a more modular design, that could prevent the propagation of errors from the ODTS to the composite clock and from the composite clock to the ODTS. For this reason, a loosely coupled approach, in which the ODTS and the composite clock are executed in totally independent modules, have been selected. With this approach, the ODTS can be also independently configured to use as reference time the composite clock, or to switch to a master clock, for example in case the composite clock needs to be reinitialized.

If the composite clock is executed outside the ODTS, an alternative source of clock differences would be needed. The use of predicted clocks from the ODTS could be considered, but to avoid a degradation of the composite clock performance the use of predictions at a very short time would be needed, compromising the separation of composite clock and ODTS modules. An additional module of time transfer is therefore considered, that allows to perform a preliminary estimation of the clock difference between satellite and station clocks independent from the ODTS. This module uses a GNSS time transfer with Galileo satellites, based on the same GNSS observations used in the ODTS and fixing the orbits to those predicted by the ODTS. Even if the use of predicted orbits from the ODTS introduces a dependency between composite clock and ODTS, the validity of the predicted orbits is much longer than that of the predicted clock, so in case the quality of the new predicted orbits cannot be ensured, the previous ones could be kept.

This logics is summarized in Figure 2. The observations from a station network and the predicted orbits of the ODTS are used to perform time transfer and obtain the differences ܭܮܥ -ܭܮܥ , including differences between satellite clocks, between station clocks and between satellite and station clocks. These differences are injected as measurements in the composite clock, that estimates models for the differences ܭܮܥ -.ܥܥ These model are used by the ODTS to refer the estimated orbits and clocks to the composite timescale. Finally, the predicted orbits and clocks from the ODTS are used to generate the navigation message broadcast to the users.

Figure 2 Proposed integration of composite clock in Galileo

In a master clock mode, the reference time in the ODTS is defined by setting one of the estimated clocks ܭܮܥ ୫ ୟୱ୲ୣ୰ െ ܵܶܦܱ ୰ୣ to zero, as per Equation ( 13). In a composite clock mode, the reference time would be defined by applying constraints to all the estimated clocks, as per Equation ( 14), so that the estimation of ܭܮܥ െ ܵܶܦܱ ୰ୣ performed in the ODTS is similar to the values of ܭܮܥ െ ܥܥ provided by the composite clock.

ܭܮܥ ୫ ୟୱ୲ୣ୰ -ܵܶܦܱ ୰ୣ = 0 (13) ܭܮܥ( -ܵܶܦܱ ୰ୣ ) -ܭܮܥ( -)ܥܥ = 0 (14)

EXPERIMENTATION RESULTS

The composite clock algorithm proposed in Section 2 and the integration within a GNSS system proposed in Section 3 have been tested with data from a station network of 16 stations for the whole year 2020. Even if the observations from 16 stations have been considered for the time transfer, the ensemble of the composite clock only contained the Galileo satellite clocks and two stations, equipped with hydrogen masers and steered to UTC.

As discussed in Section 1, to reduce as much as possible the indirectly observed states, the station clocks have been modelled only with a phase state, being their frequency and frequency drift removed by the steering. For the satellite clocks, not steered, the three states of phase, frequency and frequency drift have been considered.

This scenario includes several operational events, like the outages in the satellite E19 in March (NAGU 2020006) and June (NAGU 2020012), in the satellite E24 in September (NAGU 2020015), and in the satellite E11 in November (NAGU 2020017). For the station clocks, a jump has been considered in one of the stations (STA1) in February, while the other station (STA2) was impacted by a gap in November and a jump in December.

Figure 3 shows the Allan Deviation of the clocks resulting from this experimentation. It can be observed that the stability of stations is well in line with the stability of the hydrogen masers, showing that the composite clock is not degrading the stability of the differences ܭܮܥ െ .ܥܥ

Figure 3 Allan Deviation of the clock differences CLKi-CC

The stability of the station clocks is much better than that of satellite clocks, leading to a larger weight for the stations, as shown in Figure 4. The weight associated to the satellite clocks is only relevant in case of gap or jump in one of the stations, like in February, November and December. In those cases, the weight of the other station was limited to 0.7, to prevent that a single clock drives the composite timescale.

Figure 4 Operational events during 2020 and weights of station and satellite clocks

This strategy of weighting allowed to obtain a continuous and stable composite timescale, despite the events in the station and the satellite clocks. Figure 5 shows that the composite timescale remained aligned to the station clocks within 5ns during the whole year.

For the periods when both stations were available, the composite timescale remained close to the average of the station clocks. In November, when one of the stations was impacted by a large gap, the composite clock smoothly transitioned to a different phase, in line with the remaining station.

Figure 5 Time series of station clocks with respect to CC

The jumps in the station clocks were properly detected and corrected, and no impact on the composite timescale was observed.

Similarly, the exclusion of some satellites like the E19 and the E11, and the frequency jumps in the satellite clocks, as those observed in Figure 6, did not lead to any alteration of the composite timescale.

CONCLUSIONS

An algorithm of composite clock based on a Kalman filter can be a good solution for generating a GNSS system time. However, mechanisms to ensure the robustness of the algorithm in case of operational events, like gaps, outliers or jumps in the contributing clocks, are required. The proposed solution, based on:  the processing of a set of independent clock differences without any privileged clock,  the addition of constraints in the filter following the IGS approach to remove the non-observable states,  the use of indirectly observable states only for those clocks that require them,  the detection and correction of outliers and jumps in a preprocessing step,  the alignment of input clocks to allow for smooth reconfigurations of the composite clock, and  the anomaly detection based on the residuals of the estimation, is able to provide good results in a long-term scenario representative of the real operation of a GNSS system, using the satellite clocks and a limited set of station clocks.

The integration with the other functions of a GNSS system in a modular way allows to isolate as much as possible any potential problem, increasing the robustness of the system. It also allows to keep both the composite clock and the master clock modes, and to be able to switch from one to the other only by changing the constraints in the ODTS.

From a performance point of view, the stability of the composite timescale can be better than that of the most stable clocks, and in any case much better than that of the satellite clocks. This implies that in case the composite timescale is used as reference time for the broadcast satellite clocks, the predictability of the offsets would be driven by the predictability of the satellite clocks, minimizing the ranging error at user level.

Figure 6

 6 Figure 6 Time series of Galileo satellite clocks with respect to CC, after removing a quadratic model

engineering activities for the Galileo Second Generation, always focusing on the navigation and performance aspects. Currently, he is managing the engineering of Galileo Mission Segment in Thales Alenia Space.


Abdelkader Gouaïch
email: gouaich@lirmm.fr

Nadia Hocine
email: hocine@lirmm.fr

Liesjet Van Dokkum
email: van-dokkum@univ-montp1.fr

Denis Mottet
email: mottet@univ-montp1.fr

Digital-Pheromone Based Difficulty Adaptation in Post-Stroke Therapeutic Games

Keywords: I.2.1 [Applications and Expert Systems]: Games Difficulty adaptation, post-stroke rehabilitation, serious game

In this paper, we propose a dynamic difficulty adaptation approach for serious games dedicated to upper-limb rehabilitation after stroke. The proposed approach aims to provide a personalized rehabilitation session in which the training intensity and challenges can be adapted to patient"s abilities and training needs. The objective is to increase the rehabilitation volume by getting the patient engaged in the therapy session. This approach has been implemented and tested through a point and click game. It has been also experimented on healthy people in order to explain how this approach will be integrated to post-stroke therapeutic games.

INTRODUCTION

Stroke is a medical emergency that can cause permanent neurological and functional damages. It can cause amongst others hemineglect, the inability to move one or more limbs on the controlateral side of the body and to understand or formulate speech. The specific abilities that can be affected by stroke depend on the location, type and size of the lesion. Each patient is characterized by a specific combination of deficits.

Most rehabilitation techniques are founded on the principles of motor learning and skills" acquisition established for the healthy nervous system. Several studies suggest that intensive training (many repetitions) while giving feedback and motivating patients can have an important impact on their motor skill recovery [3] [10]. During rehabilitation sessions, therapists try to improve patients" functional skills by defining a set of exercises to reach or grasp objects in a 2D work plan. Therapists try to train patients by adapting the activities to the patient"s abilities and current health condition.

Therapeutic games or "theragames" may be considered as useful rehabilitation tools as they allow personalizing the rehabilitation session according to the patient"s abilities and training needs. In fact, theragames provide an environment in which the training intensity, difficulty, duration and frequency can be manipulated and enhanced. They also offer the opportunity for therapists to monitor and analyze rehabilitation sessions. In addition, theragames have to fulfill game design principles and rules as well, in order to be entertaining and engaging while serving their pedagogical/therapeutic purposes.

Our objective is to propose a dynamic difficulty adaptation (DDA) technique dedicated to a family of post-stroke therapeutic games. This technique focuses on a digital-pheromone based algorithm that is used to assess and track the patient"s movement capabilities.

In the following section, we will discuss theragames key challenges that have been identified by our research activities with therapists and researchers from different hospitals in Montpellier and the Movement to Health laboratory. Then, we provide readers with the state the art of difficulty adaptation in therapeutic games. In the fourth, section we describe the proposed adaptation technique as well as the digital-pheromone based algorithm. After that we will explain how this approach has been implemented by a simple training game. The six section aims to describe the pilot experiment we have conduced to control for the effectiveness of the adaptation technique in terms of motivating players. Finally, we conclude the paper by discussing the results of the pilot experiment and explaining our contribution regarding the related work.

POST-STROKE THERAGAMES CHALLENGES AND REQUIREMENTS

Theragames may be considered as a promising rehabilitation tool, especially in post-stroke rehabilitation context. The usefulness of such tool can be can be established by overcoming adequately the challenges identified by: (i) the medical aspects of the therapy and (ii) the technical, psychological and creative aspects of the game.

From a therapeutic point of view, the theragame should provide a motivating rehabilitation environment allowing for the manipulation of training duration, intensity and frequency. In addition, a theragame should take into account patients" abilities and training needs in order to provide them with a personalized game experience. Actually, each patient can be characterized by a specific combination of deficits and has a particular recovery rate. Therefore, theragames should be adaptive to the patient"s observed abilities, performances and her or his health condition.

Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. To copy otherwise, or republish, to post on servers or to redistribute to lists, requires prior specific permission and/or a fee. IHI'12, January 28-30, 2012, Miami, Florida, USA.

Another challenge identified, is the accessibility of the theragame. In fact, stroke patients would need different devices depending on their stage of recovery and rehabilitation. For instance, during the first weeks of their rehabilitation program, patients often cannot grasp an object due to their motor deficits. Thus, the theragames devices that should be anticipated for rehabilitation are for instance: motion capture based devices, IR sensors, graphic tablet, wii-board, kinect, etc.

From the game design point of view, it is important to create a game without being too much constrained by therapeutic requirements. In fact, creativity and imagination are very important in the game design to provide meaning to both gaming and therapeutic tasks. Allowing thus to create an engaging experience for the patient. Indeed, as a stroke causes long-term impairments, it is important to build several games (and consequently try to reduce the development costs) or imagine game scenarios that can sustain several months of rehabilitation. It is necessary to mention the importance of providing game designers with means to meet therapeutic requirements without being overloaded with medical details.

In order to meet these requirements, an adaptive approach could be useful in the early stages of the therapeutic game design. We focus in this paper on the dynamic difficulty adjustment of therapeutic tasks that takes into account the patient"s abilities, motivation and training needs. In addition, the proposed approach allows building low cost theragames without the necessity of having a prior knowledge about the medical details.

RELATED WORK

We can find in the literature different works that deal with poststroke rehabilitation based on a serious gaming strategy.

For instance, Heuser et al [START_REF] Heuser | Tele-Rehabilitation using the Rutgers Master II glove following Carpal Tunnel Release surgery[END_REF] suggest five therapeutic gaming exercises aiming to train impaired hand movements based on reaching and grasping targets. The adaptation technique consists in providing real-time feedbacks to the patients. Each difficulty level differs from the others by changing the color, speed and the size of virtual objects. Furthermore, Chen Y. et al [START_REF] Chen | A dynamic decision network framework for online media adaptation in stroke rehabilitation[END_REF] propose a framework for media adaptation in task-oriented neuromotor rehabilitation based on biofeedback [START_REF] Huang | Novel design of interactive multimodal biofeedback system for neurorehabilitation[END_REF]. This technique is based on the adaptation of players' feedback according to their performances. Scenery images and their clarity, sounds variety and music instruments are used to reflect the arm movement. The adaptation technique depends on the choice of the musical instrument, the variation in tempo and the task difficulty, which is annotated by the therapist. Like the previous work, the difficulty level is designed in a static manner. This can lead to a situation in which the patient is confronted with continuous failure when the task is too difficult or continuous boredom when it is too easy. Both situations can decrease motivation to continue the therapy.

Ma et al, [START_REF] Ma | Adaptive virtual reality games for rehabilitation of motor disorders[END_REF] propose an adaptive post-stroke rehabilitation system based on a virtual reality game. This therapeutic game aims to improve (i) functional arm capabilities, (ii) visual discrimination and (iii) selective attention using three difficulty levels: beginner, intermediate and expert. The dynamic adjustment of the system focuses on an in-game difficulty adaptation based on patient"s progress. This is used to select tasks and configure the initial difficulty level of the game. Actually, the system suggests a difficulty level to the patient according to her/his model. This model includes mainly information about the patient's motor skills by using metrics such as the standard medical motricity index [START_REF] Demeurisse | Motor evaluation in vascular hemiplegia[END_REF]. The adaptation technique is based on a specification matrix that determines the relationship between the patient"s profile parameters and the difficulty levels.

Another interesting example of dynamic difficulty adaptation is proposed by Cameirao et al [START_REF] Cameirao | Neurorehabilitation using the virtual reality based Rehabilitation Gaming System: methodology, design, psychometrics, usability and validation[END_REF]. The authors propose a personalized training module based on a rehabilitation scenario, the so-called Spheroids. The adaptation consists in providing an appropriate task difficulty by capturing specific features of the movements of the arm. The adaptation technique focuses on an assessment step which represents the system calibration. The player is asked to move her hands to numbered dots positioned at specific locations on the tabletop. In each trial, hand and target position is randomly defined by the system.

In the above presented work, the game difficulty is often annotated by therapists in the beginning of the therapy session. Thus, when the patient fails to continue to play when the game difficulty is beyond her/his capabilities, engagement in the therapy and motivation to continue decrease. In addition, we can mention the lack of techniques allowing the reusability of game elements in the game framework. This is important as it allows building low cost theragames, while introducing variability in therapeutic activities that can also improve patients" motivation.

PROPOSITION

This section describes the proposed adaptation technique that aims to meet both therapy and game design requirements. The overall difficulty adaptation process is presented in figure 1.

A virtual rehabilitation session based on therapeutic serious game may provide an appropriate and helpful environment for both patients and therapists. The serious game can be helpful as it consists of a training tool, which can be used during: (i) the first months of in-centre stroke rehabilitation and (ii) when patients return home to help them to remotely continue their rehabilitation program in order to improve their quality of life.

Like in actual rehabilitation sessions, the game proposes tasks to the patient that allow her/him to perform movements in a 2D plan (a table) the so-called work plan. Due to the stroke, the patient finds it more or less difficult to reach certain regions of the work plan. In the experiment conduced, only pointing tasks, whereby the patient is asked to reach towards a target, are considered. Indeed, for post-stroke rehabilitation reaching represents an important part of upper-limb rehabilitation. When the patient is able to achieve pointing tasks, the therapists can focus on other types of tasks such grasping.

The first step of the adaptation process consists of collecting initial data of the patient"s functional abilities. This allows constructing a model about patient"s movement abilities on the work plan called "ability matrix".

The actual game needs information in order to adapt the difficulty of game tasks allowing movements towards certain regions of the work plan. In our adaptation model, this information is provided by a structure called "interaction matrix". The latter represents the work plan that is modelled by a matrix in which cells contain a Boolean value stating whether or not it is recommended to situate a target at this place.

The interaction matrix generator is defined by two parameters: the patient"s ability matrix and a recommendation from the motivation module to either increase decrease or maintain the difficulty level of tasks.

Figure 1: Overall process of dynamic difficulty adaptation

We describe in following sub-sections the dynamic difficulty adaptation components and their integration within a theragame application.

Player's Profile and the Abilities Assessment Step

In order to make decisions on the therapeutic objectives, therapists use generally clinical test scores, such as Fugl-Meyer score [START_REF] Fugl-Meyer | The post-stroke hemiplegic patient. 1. a method for evaluation of physical performance[END_REF], or a kinematic based evaluation in virtual reality therapy. As for serious games, a kinematic based evaluation can be used and especially in our approach, the assessment step is included in the game. The game asks the player to reach for targets placed at different positions in the work plan. All information on execution of this initial task is stored in the patient"s specific profile, alongside more general information like age, gender, dexterity, and impaired side, date of stroke, initial fugl-meyer score and neglect. The ability data is modelled by a structure so-called "ability zone". The player"s profile is then used to choose the appropriate adaptation strategy according to the objective of the proposed game.

Ability Zone

Definition: An ability zone 𝐴(𝑚, 𝑛) is a matrix of a dimension 𝑚 × 𝑛 where m represents the number of rows and n the number of columns. Each cell 𝐴 𝑚, 𝑛 ∈ ℝ value represents an easiness score for the patient to reach the mapped area of the work plan.

The ability zone represents the patient"s movement capabilities model within a 2D work plan. The goal of constructing such model is to build an instantaneous image of the patient"s movement capabilities.

While the definition of the ability zone provides only its structure and intended semantics, it does not provide any means to construct this structure and different approaches can be used. For instance, in our previous work, a statistical approach is used in which each cell 𝐴(𝑖, 𝑗) contains a probability of success of the patient when asked to reach the mapped area [START_REF] Hocine | Therapeutic games' difficulty adaptation: An approach based on player's ability and motivation[END_REF].

In this paper, we suggest using a bio-inspired approach to build the ability zone using a virtual ant.

Digital-pheromone based algorithm

The ant algorithm was originally introduced by Dorigo [START_REF] Dorigo | Ant Colony Optimization[END_REF] as a metaheuritic to solve problems that can be reduced to an optimal path finding when exploring a graph. The principle of this approach can be summarized as follows:

1. The tackled problem is modelled as a graph consisting of a set of nodes and links between these nodes.

2. Virtual ants are then introduced to explore nodes of the graph following links with the highest level of pheromones.

3. In fact, when exploring the graph, ants can deposit "pheromones" to indicate an interesting "finding". Following this behaviour, all interesting paths of the graph will be marked with high levels of pheromones and all other links will have low levels because pheromones are volatile.

In our context, the principles of ant algorithms are used to construct the patient"s ability zone as follows:

1. A virtual ant is placed on the hand of the patient and follows its movements;

2. By moving her/his hand, the patient makes the ant explore the work plan;

3. The ant deposits pheromones following the movement path of the patient: the pheromones represent a signal that can be propagated to the work plan following a propagation law. Furthermore, the pheromone can also evaporate following an evaporation law. All areas that have been reached many times by patients are considered as easy and the others as difficult.

Propagation law:

The propagation law is an important feature since it allows predicting the difficulty of nearest cells without requiring the patient to actually explore them. The predicted difficulty is less important than the actual one and the pheromone signal has to be accordingly decreased. In our approach we focused on a linear propagation function. This function takes as parameter a range w (defined according to the wanted precision of the movement) and propagates the pheromone signal to all cells within the w-range of the current position. The level of the pheromone that reaches the cell, c, when the ant is located at cell, s, is calculated as follows:

𝑙𝑒𝑣𝑒𝑙 𝑠 𝑐 = -𝜑 𝑤 𝑑𝑖𝑠𝑡 𝑠, 𝑐 + 𝜑 , if 𝑑𝑖𝑠𝑡(𝑠, 𝑐) ≤ 𝑤 0, otherwise
𝜑 represents the initial (or nominal) level of the pheromone deposited at the location of the ant and 𝑑𝑖𝑠𝑡 𝑠, 𝑐 indicates the distance between the cell c and the location of the ant s.

Evaporation law:

The evaporation law is as important as the propagation law since it allows forgetting regions that the patient has reached by chance. In fact, since the patient"s movements are not totally controlled, the patient may reach involuntarily some areas of the work plan. The volatile nature of pheromones allows ignoring this information, which may disturb the decision making on difficulty. The memory of the evaporation law is defined by task parameters according to level of details chosen by the therapist.

Improving the method accuracy:

To be accurate, proposed method has to take into account the movement as well as its goal. In fact, successful movements where the player reaches the target with "a good quality", have to be given more importance than unsuccessful movements where the player does not reach the target or reaches it with inaccurate and uncontrolled movement path. For this reason, the nominal pheromone signal deposited by the ant is modulated by a movement indicator that expresses its quality. Within the context of this study we have used the following model to modulate the level of the pheromone.

Knowing that the patient has reached the target moving along the path p, the pheromone model can be described as follow:

𝜑 ′ =

Dynamic difficulty adjustment technique

General framework

In [START_REF] Hocine | Therapeutic games' difficulty adaptation: An approach based on player's ability and motivation[END_REF] we have introduced a motivation model based on job satisfaction and activation theory to adapt task difficulty during rehabilitation sessions. In this approach motivation is considered as a regulation process (see figure 3). When the reality (outcome of a therapeutic game) is consistent with the expectations of the patient (perception of own motor skills) then this is considered as a stabilized satisfaction state. The patient attempts to provide effort to maintain this steady state. In the case of a minor disruption (constructive dissatisfaction), inducing challenge, the patient may produce efforts to regulate her/his satisfaction state.

The difficulty adaptation process aims to produce a constructive dissatisfaction situation to keep the patient engaged in the therapy session, which can increase the rehabilitation volume (i.e total time and number of therapy session)

Depending on successes or failures of the patient the motivation module provides one of the following recommendations:

 Decrease task difficulty: this happens when the player is confronted with a situation of consecutive failures; it is important to decrease difficulty level in order to avoid her/his frustration.



Increase task difficulty: when the player is in a situation of consecutive successes and finds the game very easy, it is important to increase the difficulty to maintain her/his engagement and avoid boringness.



Maintain the level task difficulty: the player is reaching a steady state that should be maintained.

In order to get more details about the motivation module to make these decisions see [START_REF] Hocine | Therapeutic games' difficulty adaptation: An approach based on player's ability and motivation[END_REF]

Interaction Matrix

Definition: An interaction matrix of a dimension 𝑚 × 𝑛 is defined as 𝐼(𝑚, 𝑛) where m represents the number of rows and n the number of columns. Each cell 𝐼 𝑚, 𝑛 ∈ {0,1} contains a binary value indicating whether or not an interaction with the player is desired within this cell.

The interaction matrix is a structure that models the 2D work plan used by patients to perform their rehabilitation tasks. Depending on the pointing device used (IR emitter, graphics tablet, mouse, Kinect) and the specific movement capability, the actual working area becomes more or less important. For this reason, the actual size of the work plan can be estimated to 1.5m width and 1.5m height.

The work plan is mapped to the action matrix by defining a scale used to translate each element (𝑖, 𝑗) of the action matrix to a rectangular area defined by the following indices, left-bottom and right-top points: (i * cell_height, j * cell_width)

and ((i + 1) * cell_heigh, ( j + 1) * cell_width)

The interaction matrix can be used to either: (i) challenge the patient, (ii) support her/him or (iii) maintain the current state. The interaction matrix will influence positions of game objects in the virtual word. Depending on the motivation module and the patient"s ability matrix, the decision making difficulty levels can influence the interaction matrices.

Case 1: Increase difficulty

When the motivation module suggests increasing the difficulty, the player is succeeding too easily, which may produce boringness and decreases her/his engagement. The goal is to explore more difficult cells. However, as said previously, it is important to create a constructive dissatisfaction without selecting very difficult areas of the work plan that are represented by the outside frontiers of the ability zone. In order to take this decision on difficulty, a simple gradient based algorithm is used, in which the gradient of each cell indicates direction of increase or decrease of the pheromone:

𝑔𝑟𝑎𝑑(𝐴 𝑖,𝑗 ) = (𝐴 𝑖,𝑗 -𝐴 𝑘,𝑙 ) 𝑘,𝑙∈𝑖-1..𝑖+1,𝑗 -1..𝑗 +1
The interaction matrix is constructed by the following challenge function that simply keeps cells with negative gradient:

Challenge(𝑋) 𝑖,𝑗 = 1, 𝑔𝑟𝑎𝑑 𝑋 𝑖,𝑗 > 0 0, otherwise

Case 2: Decrease difficulty

When the player is failing too often, as the exploration of the work plan is too large in regards to the player"s current capabilities, the motivation module suggests decreasing the difficulty. The solution is then to come back to known areas where the player was succeeding. This can be established by selecting the inner frontier of the current ability zone.

In this case the interaction matrix is constructed by a support function that keeps cells with positive gradient:

Support(𝑋) 𝑖,𝑗 =
1, 𝑔𝑟𝑎𝑑 𝑋 𝑖,𝑗 < 0 0, otherwise

Case 3: Maintain difficulty

In this case the motivation module recommends maintaining the current task"s difficulty level. Consequently, the previous interaction matrix is maintained.

Example:

Let suppose that a player has performed a movement path when she/he was asked to reach a target positioned at (3,3) (see Figure 5). When a virtual ant follows this path (with nominal pheromone level=1, diffusion parameter W=2 and evaporation=0.5) constructs the ability zone presented in Figure 5 (a).

In a challenging mode, the interaction matrix, which represents cells with negative gradients, is shown in figure 5 (b). In a supporting mode, the interaction matrix, which represents cells with positive gradients, is shown in figure 5 (c).

IMPLEMENTATION

The technical architecture used for integrating the adaptation mechanism through theragames is presented in Figure 6. The proposed game, used as a case study, has been developed using Unity3D. This game acts as a client and communicates with the dynamic difficulty adjustment (DDA) server using REST protocol. All information is stored using a NoSql server (CouchDb)

On one hand, the game traces patient"s movement which allows the difficulty adjustment service to update the ability and interaction matrices. On the other hand, the interaction matrix is used by the theragame to determine targets" positions. In fact, they are placed as close as possible to recommended cells.

The DDA server makes it possible to play the theragame and dynamically adjust the difficulty of pointing tasks based on the observed patient"s capabilities.

It is important to note that the game designer has no prior knowledge about therapy requirements. She or he has simply to take the recommendations of the interaction matrix into account, when building the scene of a game stage. These results represent a loose coupling between the therapeutic level and game level. Finally, the DDA server allows the therapist to access a set of monitoring services in order to consult patient"s profile, performances and the evolution of its ability zone during the game session.

The presented approach is generic and can be applied to any pointand-click game (based on pointing tasks) adapted to stroke patients. 6. PILOT EXPERIMENT

Protocol and participants

To experiment with the presented approach we have executed a pilot study on healthy players. In fact, before testing on patients and disturb their classical rehabilitation program (planned for five weeks at least), it was necessary to experiment on healthy participants. Nonetheless, the experiment scheme with healthy players has to simulate difficulty of pointing tasks. 2) To reach a target (1) she/he has to extend the rubber by continually hitting the space key on the keyboard using her/his non-dominant hand. More the frequency of hit is high, more the virtual rubber will extend; if the hit frequency diminishes then the rubber length diminishes to its minimal size. The mouse only provides the direction of the virtual rubber and is controlled by the dominant hand.

The experiment follows a repeated-measures design where each player plays with a random strategy of difficulty adaptation and our proposition for difficulty adaptation. For all players the experiment was conducted in a similar way. An initial assessment allows collecting the maximum hit frequency of each user. The player then plays an assessment game by cleaning the screen using the virtual rubber. This builds the initial ability zone. Then, the player is asked to reach a sequence of targets in two sessions: session "A" that uses a random algorithm and a session "B" using our proposition to place the target. The order of these sessions has been alternated between players. The player was given a rest between these two sessions to avoid fatigue. Finally, the player reports on her/his own perceived difficulty, frustration and fatigue using the DP-15 scale [START_REF] Deligneieres | La perception de la difficulte dans les taches motrices[END_REF] for both sessions A and B. It is worth noting that both sessions have the same duration of 2min.

The experimental group was composed by n=10 participants with the following characteristics: We have introduced the following experimental hypotheses:

 I.H0 There is no difference between session A and session B concerning the success of pointing tasks  II.H0 There is no difference between session A and session B concerning the level of perceived difficulty, frustration and fatigue.

Results

Results are reported as (Mean ± Standard Deviation). All statistical analyses were performed using R (http://www.rproject.org) version 2.12.0. A repeated measure t-test was used to reject hypothesizes.

IH0 was rejected since our proposition of difficulty adaptation has increased the number of successful pointing tasks by an average of M=4.3 with SD=3.5. The increase was statistically significant: t(9)=3.85, p<0.01 and r2=.97 II.H0 has not been rejected. This means that there is no statistically significant difference between session A and B concerning the perceived difficulty, frustration and fatigue.

Discussion

I.H0 hypothesis stating that there is no difference concerning the number of successful pointing tasks between session A and B has been rejected. We reason that the adaptation strategy has increased the number of successes by placing targets on the edge of the ability zone. It is important to note that tasks of session B have not been perceived as less difficult than tasks of session A by users. In fact, II.H0 shows that users have perceived no difference in difficulty of tasks between session A and B. This can be explained by the fact that the adaptation strategy places targets at the edge of ability zones maintaining consequently the challenge. During qualitative interviews, all users have reported fatigue and frustration when controlling the virtual rubber by taping on the keyboard using the non-dominant member. This shows that the pointing task was challenging.

CONCLUSION

In this paper, we described a theragame framework aiming to meet both therapy and game requirements. The usefulness of such rehabilitation tool consists in providing a personalized and motivating rehabilitation session in which the training intensity and challenge are adapted to patient"s abilities and training needs. We focused in this work on the adaptation of therapeutic game difficulty. The proposed difficulty adaptation technique is generic and dedicated to a family of therapeutic games for upper limb rehabilitation. The proposed technique is based on both the player"s profile and her/his performances to dynamically adapt the game experience. In order to test the effectiveness of such approach we initially executed a pilot experiment on healthy players using a random strategy and the proposed adaptation strategy to place targets.

Results showed that the proposed adaptation strategy increased the total number of successful pointing the same amount of time, compared to a naive random algorithm.

Despite the large effect size of 79%, results of this experiment have to be carefully handled due to the small sample size (n=10) Further experiments with stroke patients have to be executed to confirm these initial findings.

Figure 2

 2 Figure 2 presents an example where the ant has followed patient"s movements. Grey scales are used to indicate the different level of the pheromone (the light grey cells represent cells were successfully reached). As shown in (b) the volatile nature of pheromone helps reducing noises coming from uncontrolled movements. This image provides a simple indication of the patient"s ability zone.

Figure 2 :Figure 3 :

 23 Figure 2: Example of pheromones deposited by an ant following the patient's hand (a); the application of evaporation law reduces noises due to uncontrolled movements (b)

Figure 5 : matrices of example figure 4 Figure 6 :

 546 Figure 6: Technical architecture DDA server with the game.

Figure 4 :

 4 Figure 4: Example of a movement path.

Figure 6 :Figure 6

 66 Figure 6: Screenshot of the pointing game used experiment

Figure 7 :

 7 Figure 7: Assessment, the user is asked to erase as quickly as possible the screen by moving the rubber band

Table 1 . Pilot experiment participants

 1

	Gender	Age	Dominant hand
	7M 3F	(25.9 ± 4)	8 right; 2 left