Gabriel Vasile
email: gabriel.vasile@gipsa-lab.grenoble-inp.fr

Rayna Guy D'urso
email: guy.durso@edf.fr

Charlatchka

Elena Lungu
email: elungu@sigintec.fr

Calibration of an Active Ultrasound Bedload Monitoring System for Underwater Environments

This paper introduces the concept of coupling ultrasound transducers with piezoelectric accelerometers for active ultrasound bedload monitoring in underwater environments. A general system calibration procedure is proposed which allows the enduser to establish a real-time link between the energy of the received signals and the bedload transport phenomena.

I. INTRODUCTION

Monitoring the bedload transport in underwater environments is crucial for understanding stability issues such as reservoir silting or channel self-cleaning. In this context, sediment traps are widely used to derive the sediment balances. The Birkbeck sampler has become one of the preferred method for in situ bedload measurements [START_REF] Laronne | The continuous monitoring of bedload flux in various fluvial environments[END_REF]. Being an integral measure of the transported sediment volume in a certain time interval, this direct method suffers from one main drawback: it has a finite capacity. Without the use of devices such as a sludge pump to empty the accumulating sediment, the sampler fills and eventually ceases to yield data.

An alternative nonintrusive bedload monitoring instrument is the buried geophone station [START_REF] Rickenmann | Bedload transport measurements at the Erlenbach stream with geophones and automated basket samplers[END_REF], [START_REF] Rickenmann | Bedload transport measurements with impact plate geophones: comparison of sensor calibration in different gravel-bed streams[END_REF], [START_REF] Hilldale | Installation of impact plates to continuously measure bed load: Elwha river, Washington, USA[END_REF]. Under protection of steel plates, several geophones can provide continuous and automatic measurements, even during large floods. Every stone passing the steel geophone equipped plate generates an impulse recorded by the signal acquisition board. Using a calibrated voltage thresholding scheme the grain impacts are recorded and counted. One can notice that this instrument is highly sensitive to the environmental noise, since it is operating at low frequencies (< 1KHz).

In [START_REF] Bogen | Bedload measurements with a new passive acoustic sensor[END_REF], Bogen and Moen used piezoelectric acoustic transducers for bedload monitoring stations in Norway. The reported results present a clear correspondence with the sediment flow, but no estimation scheme is proposed.

In this paper, the potential of ultrasound (US) transducers coupled with piezoelectric accelerometers for bedload monitoring stations is further explored in a controlled laboratory environment. The main objectives are:

• to increase the operating frequency of the acoustic transducers, • to provide consistent calibration and processing,

• to evaluate the uncertainty of the derived measurements.

To pursue these objectives, an experimental platform, illustrated in Fig. 1-(

II. CALIBRATION AND MEASUREMENT

As illustrated in Fig. 1-(b), the following sensors have been employed: 4 ultrasound contact transducers PAC R50 (EM, R2, R3, R4), 2 calibrated piezoelectric accelerometers Endevco 233E (A1, A1) and 1 geophone R.T. Clark (G). The EM US transducer has been excited using the ENI 2100L power amplifier and the Picotest G5100A arbitrary signal generator. The received signals from A1, A2, R3, R4 have been conditioned using the Nexus low noise amplifier and recorded using the PXIe-1082 NI system. For the active configuration, the multi-sensor platform has been calibrated both in static and dynamic regimes.

The static calibration is a direct relationship between the output and the input of the system determined experimentally. Static calibration requires a reference standard, to which the output is reported. For our system of sensors, the calibration process determines a relationship between the input voltage controlling the EM US transducer and the output characteristic, called static characteristic curve or calibration. Fig. 2 shows the static characteristic in water, measured by transmitting a sinusoid signal and looking at the reception of the peak value of the Fourier transform. The value obtained supports the idea that our system is linear in all cases. In circumstances where the assumption of stationary is unsustainable, linear dynamic models need to be replaced by more sophisticated models. The linear theory provides a point of departure for such developments. The dynamic characteristic of each sensor is related to the performance when the input EM voltage signal is a function of time. The speed of response of a system refers to the ability to respond to sudden changes in the amplitude of the input signal. Thus, the response rate is calculated based on the parameters describing the system performance under the transition: constant measuring of the time delay (latency measurement), setup time (settling time), the dead time and range dynamics.

For a second order system, the transfer function can generally be expressed as:

G(s) = Y (s) X(s) = b 0 a 2 s 2 + a 1 s + a 0 , (1)
where the X(s) and Y (s) are the Laplace transforms of the input and output signals.

The received signal characterizing the two transducers types are illustrated in Fig. 3 One can notice that the system settling time of the two accelerometers is greater than the settling time of the two ultrasound transducers. As a consequence, the triggering signal is set to be the signal received from the two accelerometers.

To obtain the performance parameters of a second order system, one can apply the Hilbert transform in order to retrieve the envelope of the received signal. According to Fig. 4, the important parameters are:

1) the payback period (T d ) -time between two consecutive peaks after exceeding 2) the overshoot (OS) -difference between the maximum signal value and the final value

OS = y max -y f in = (y f in -y 0 )e -ξπ √ 1-ξ 2 , (2)
3) the depreciation rate (ξ) -describes how fast the signal is stabilized ξ

1 -ξ 2 = 1 π ln y f in -y 0 OS , (3)
4) the natural damping frequency (ω n )

ω d = 2π T d = ω n 1 -ξ 2 , (4)
5) the static sensitivity (K) -difference between the amplitudes of the received signal before and after the settling time

K = y f in -y 0 C -C 0 . ( 5
)
The obtained dynamic characteristic are illustrated in Fig. 5 and match the previously derived static characteristic. Finally, the transfer function of an acoustic transducer can be computed by positioning two transducer elements in transmission configuration:

1) by continuous excitation of one transducer with wide band white Gaussian noise, the other responds in its nominal frequency band. 2) by considering the performance parameters illustrated in Fig. 4-(b), one can define an impulsive excitation for estimating the system transfer function. Accordingly, by taking the impulse duration (short sine waveform) to be smaller than the system stabilization time, the transducer is sensing this calibration signal as an impulse and not as a step function. In this paper, the second method for estimating the transfer function is employed and the results are illustrated on Fig. 6.

III. RESULTS

Four configurations have been tested in the controlled water tank facility :

• passive configuration in the nominal bandwidth of the impact (< 14 KHz), • passive configuration outside the nominal bandwidth of impact (20 to 40 KHz), • active configuration with continuous sine excitation at 21 KHz, • active configuration with chirp excitation at (20 to 40 KHz). After using controlled steel balls for calibration, two types of sediments were investigated and confidence intervals for the mass measurements have been obtained. Notice that the maximum speed in water have been reached before hitting the multi-sensor bedload monitoring platform in all cases. An Figs. 8,9, 10 and 11 illustrate the obtained results. One can observe that the active ultrasound bedload monitoring platform can provide similar results (in terms of measurement error and confidence intervals) as conventional passive bedload monitoring in controlled noise free environment. This opens a new opportunity to increase the robustness of bedload monitoring system in underwater environments, where the frequency bands below 15 KHz are quite sensitive to the environmental noise.

IV. CONCLUSION

The general system calibration procedure presented in this paper offers an alternative option to monitor bedload transport phenomena in underwater environments. The main advantage of the proposed experimental setup is better noise robustness at lower acoustic frequency bands. Future studies include the full performance analysis in real underwater conditions.

Fig. 1 .

 1 Fig. 1. Multi-sensor bedload monitoring platform: (a) general view, (b) setup of the acoustic sensors on the top steel plate.

Fig. 2 .

 2 Fig. 2. Multi-sensor bedload monitoring platform: static characteristics.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Multi-sensor bedload monitoring platform: step response.

Fig. 5 .

 5 Fig. 5. Multi-sensor bedload monitoring platform: dynamic characteristics.

Fig. 6 .

 6 Fig. 6. Transducer transfer function: (a) accelerometer, (b) ultrasound transducer.

Fig. 7 .

 7 Fig. 7. Multi-sensor bedload monitoring platform: spectrogram of the received signals for steel ball calibration.

Fig. 8 .

 8 Fig. 8. Bedload mass measurements with passive configuration in the nominal bandwidth of the impact (< 14 KHz).

Fig. 9 .

 9 Fig. 9. Bedload mass measurements with passive configuration outside the nominal bandwidth of impact (20 to 40 KHz).

ACKNOWLEDGMENT

The authors would like to thank the Électricité de France company for funding this research.


Chihebeddine Ammar -Kais Haddar
email: ammarchihebeddine@hotmail.com

Patents editor in order to automatically construct terminological databases

Keywords: Normalisation, modèle de document de brevets, Bases de données terminologique Normalization, Document patent model, Terminological databases

Les demandes de brevet ont une structure similaire dans le monde entier. Elles comprennent une page de couverture, un mémoire descriptif, des revendications, des dessins (si nécessaire) et un abrégé. Au monde Arabe, il n'y a pas de collection numérique uniforme de document de brevets et donc pas de collection XML. Dans ce contexte, nous présentons la création d'un modèle de document de brevet standard pour les brevets Arabes et le développement d'un générateur de collection XML de brevets ayant une structure uniforme est simple à utiliser. Cette collection va nous être utile pour la construction d'une base de données terminologique pour les domaines scientifiques et techniques.

ABSTRACT. Patent applications are similarly structured worldwide. They consist of a cover page, a specification, claims, drawings (if necessary) and an abstract. In the Arabic world, there is no uniform digital collection of patent documents and therefore no XML collections. In this context, we aim to create a standardized document model for Arabic scientific patents and develop a generator of XML patent collection having a uniform and easy to use structure. This collection will be useful for us to build a terminological database for scientific and technical fields.

Introduction

One of the very rich in terminology work streams are the scientific patents. They are similar, for example, to a scale repository. They also cover several scientific and technical fields, while offering rich interdisciplinary relations. That is why we will need several terminological databases, one for each field. Indeed, standardized modeling patent allows us to maintain a standard for the representation of texts in digital form, so that we protect patents data by bringing them in digital databases. It will provide a single common data model for all terminological data regardless of the data's language, source, field, etc. Also, we will be able to build collections of uniform patents which facilitate the extraction and the exploitation of patents data and the extraction of links between valid terms. Standardized modeling patent ensure also interoperability between applications. Finally, it will allow us to easily enrich other terminological databases. Patent information, in the Arabic world, remains almost the preserve of patent agents or lawyers versed in research needed to be done before any patent application filing or preparation a lawsuit. That's why, the development of computerized databases on patents, will open access to all categories of users: businessmen, economists, researchers, etc., and make them aware of the potential value of the information that patents contain. Patents are available in different formats: Full text, PDF document, set of images, XML, etc. They have heterogeneous components that require different modeling. Also, patents have linguistic structures like text and titles, and nonlinguistic structures like figures, citations, tables and formulas. In fields such as mechanics, automatic extraction based only on the text will fail. In addition to the text, figures and citations information, all patent publications contain a relatively rich set of well-defined metadata. These metadata are often found in the cover page of patents and titles of figures and tables. To cope with the large volume of data and metadata, we will develop a patents terminological editor to generate terminological databases. This will allow us to develop heuristics, based on metadata such as the applicant(s) name(s), the inventor(s) name(s) or priority documents, etc., for finding interesting documents. The structure of the XML documents may be used for the processing performed to differentiate various elements according to their semantic. Thus, a section title, a summary, bibliographic data, or examples can be used to identify different aspects of the text. Indeed, scientific patents can be easily processed as XML documents. So we can treat their structures1 as a source of information. The work presented in this paper is the continuation of studies [START_REF] Ammar | A standard TMF modeling for Arabic patents[END_REF] on the standardization for Arabic patents. We propose a standardized model for Arabic patents and create a patents editor in order to generate a patent collection having similar structure. It is an original idea because nobody treated terminology in Arabic patents in previous works. This article is organized as follows. In section 2, we present previous works. We present, in section 3, our Arabic standard patent document editor. Section 4 is devoted to the evaluation and discussion and we conclude and enunciate some perspectives in section 5.

Previous works

Previous works on patents [START_REF] Lopez | Multiple retrieval models and regression models for prior art search[END_REF]Lopez and Romary, 2010b;[START_REF] Magdy | Exploring Standard IR Techniques on Patent Retrieval[END_REF] were mainly based on purely statistical approaches. They used standard techniques of information retrieval and data extraction. Some of the previous works use machine learning tools to extract header metadata from patents using support vector machines (SVM) [START_REF] Hoang | Extracting and Matching Authors and Affiliations in Scholarly Documents[END_REF], hidden Markov models (HMM) [START_REF] Binge | Scientific literature metadata extraction based on HMM[END_REF], or conditional random fields (CRF) [START_REF] Lopez | GROBID: Combining Automatic Bibliographic Data Recognition and Term Extraction for Scholarship Publications[END_REF]. Others use machine learning tools to extract metadata of citations [START_REF] Hetzner | A simple method for citation metadata extraction using hidden markov models[END_REF], tables [START_REF] Liu | Tableseer: automatic table metadata extraction and searching in digital libraries[END_REF], figures [START_REF] Sagnik | Figure Metadata Extraction from Digital Documents[END_REF] or to identify concepts [START_REF] Pattabhi Rk Rao | Automatic Identification of Concepts and Conceptual relations from Patents Using Machine Learning Methods[END_REF]. All these approaches rely on previous training and natural language processing. In (Lopez et al., 2010a), the authors developed multilingual terminological database called GRISP covering multiple technical and scientific fields from various open resources. The European Patent Office, EPO2, offers inventors a uniform procedure of application, and a register of multilingual patents (English, French, and German). The MAREC (MAtrixware REsearch Collection) database is formed of patent documents from the European Patent Office. It is a standard corpus of patent data available for research purposes. It consists of 19 million of patent documents in different languages (English, French, and German) in a standardized XML schema highly specialized. The ePCT 3 is a WIPO (World Intellectual Property Office) online service that provides secure electronic access to the files of international applications filed under the international patent system as maintained by the International Office. It is also possible to file international applications using ePCT-Filing. In the literature, works on Arabic patents are missing. This is caused by the fact that in the Arabic world, there is no North African or Arabian Intellectual Property Office and therefore no uniform collections of Arabic patents. In Tunisia, for example, the INNORPI4 (National Institute for Standardization and Industrial Property) does not propose a digital collection of patent documents and therefore no XML collections. As a result, Arabic patents have no unique structure. For the Tunisian patents, as illustrated in Table 1, the cover page doesn't have abstracts and patent documents could be in one of the three languages (Arabic, English or French). In the regional office5 for the Gulf Cooperation Council (GCC Patent Office6 ), there are only Arabic patents and there is an Arabic abstract in the cover page. The layout of the description part varies also from place to place. For example, the summary and the background of the invention could not exist in some patent descriptions. The Tunisian patents themselves have no unique structure in that some of them have no abstract, have missing bibliographic data and even no cover page. For these reasons, a normalization phase for Arabic patents is necessary.

Tunisian patents GCCPO patents Language Arabic, French, English Arabic

Digital document

No Yes

International patent classification

No Yes

Search report

No No

Citations

No Yes

International publication Yes No

Table 1. Comparison between Tunisian and GCCPO Patents.

Arabic standard patent document editor

Patent applications are similarly structured worldwide. They consist of a cover page, a description, claims (Hong, 2013), drawings (if necessary) and an abstract. The cover page of a published patent document usually contains bibliographic data such as the title of the invention, the filing date, the priority date, the names and addresses of the applicant(s) and the inventor(s). It also has an abstract, which briefly summarizes the invention, and a representative drawing. Bibliographic data are extremely useful for identifying, locating and retrieving patent documents. The patent description must describe the claimed invention and give technical information. The claims determine the patentability and define the scope of the claimed invention.

Patent documents are often difficult to understand and have a variety of structures. So, we propose an Arabic patent document model and develop a patent editor which automatically generates a collection of XML patent documents having a similar structure. It will facilitate the task of terms and keywords extraction.

In the following, we will present our UML based patent document model for bibliographic and application data. The structure of the patent can be divided into two parts: bibliographic data taken from the cover page and application data from the rest of the patent document. The Figure 2 shows the class diagram of the patent application data in which all associations are also strong composition associations, because, if a composite is removed, all of its component parts will be removed with it. It presents the association of the Invention class with one or more Abstract in different languages, one Claim and Description parts and zero or one Drawing part. The two above presented diagrams allow us to introduce a DTD (Document Type Definition) for scientific Arabic patents, as shown on Figure 3. The role of the DTD is to precisely define the unique structure of Arabic patent documents, no matter the patent offices.

Evaluation and Discussion

We did not have a collection of document in digital form because it is not the official in Tunisia for example. So we created our small collection of multilingual patents from various fields using our patents editor.

To cope with the large volume of patents data and metadata, we developed a patents terminological editor to generate TMF (ISO 16642:2003 -Terminological Markup Framework) terminological databases. This will enable us to facilitate the extraction and information retrieval tasks from the cover pages (metadata), and the other parts (data) of patents. Our terminological database contains terms of different technical and scientific fields and various patents with different structures. We can distinguish two categories of terms: the scientific and technical terms and the other terms. Scientific and technical terms in their turn were divided according to their technical and scientific fields.

The results of our terminological database are presented in

Conclusion

Our main obstacle is that the structure of patents differs from an intellectual property office or institute to another in the Arabic world. The cover page of a Tunisian patent differs from the Egyptian or Moroccan patent cover page. We conducted a standardized modeling for Arabic patents based on the forms of patents published in the Arabic world. It provides us a single common patent document model. We developed a patent document editor to create a collection of Arabic patents having unified structure. In future, we plan to enlarge our patents collection and then our terminological database. We will merge several terminology databases of patents. We aim to better extract information from a collection of multilingual scientific patents and to combine onomasiological and semasiological models. We are also developing a new annotation procedure, to annotate our learning and test collections. We aim also to evaluate the relevance of our terminological database with an information retrieval system for patent documents, and to realize a prior art search for the collection of patents.

Information retrieval technics in multilingual patents are not lacking in previous works, we will test whether the results of this works remain valid if one expands the collection by documents into other languages (Arabic, for example), and if they will be affected by changing the type of the document collection, calculating noise, redundancy, cost, precision, recall, silence, etc.

Figure 1

 1 Figure 1 shows the class diagram of the patent bibliographic data in which all associations are strong composition associations. It contains Bibliographic Data class which includes the Filing Number and Date, the Publication Date and Language and Classification of the patent. Bibliographic Data object is associated with one or more Title of Invention in different languages, zero or more Priority patent applications, one or more Inventor(s) and Applicant(s), zero or one Representative (agent) and zero or more International Publications (PCT).

Figure 1 .

 1 Figure 1. The class diagram of patent bibliographic data.

Figure 2 .

 2 Figure 2. The class diagram of patent application data.

Figure 3 .

 3 Figure 3. A DTD for scientific Arabic patents

Table 2 .

 2 Table 2 (only technical and scientific terms are counted). It concerns Tunisian and Gulf Arabic patents and it can be easily merged with other terminological databases. We hope that our terminology database will improve patent search. Over view of the number of technical and scientific terms in our terminological database.

	Collection Number of patents	Number of terms Full text Cover page Abstract
	INNORPI	28	236	25	63
	GCCPO	30	302	312	96

Remind that an XML document is structured as a tree consisting of hierarchical elements which may have one or more attributes, the leaf nodes have information.

EPO: European Patent Office, http://www.epo.org/

https://pct.wipo.int/LoginForms/epct.jsp

INNORPI: National Institute for Standardization and Industrial Property, http://www.innorpi.tn

Certificates of Patents granted by the GCC Patent Office secure legal protection of the inventor's rights in all Member States.

GCC Patent Office: Gulf Cooperation Council Patent Office, http://www.gccpo.org/

<?xml version="1.0" encoding="UTF-8"?> <!ELEMENT Invention (Bibliographic_Data, Description, Claim, Drawing?, Abstract+) > <!ATTLIST Invention ID_Inv CDATA #REQUIRED Inv_Lang (AR|FR|EN) "AR"

Drawing (Figure+) > <!ATTLIST Drawing ID_D CDATA #REQUIRED > <!ELEMENT Abstract (p+) > <!ATTLIST Abstract ID_Ab CDATA #REQUIRED Num_Ab CDATA #REQUIRED >


Eric Guy D'urso
email: guy.durso@edf.fr

Jérôme De Oliveira

Guillet

Elena Lungu
email: elungu@sigintec.fr

This paper introduces the concept of automatic reference selection for active ultrasound wild salmon monitoring systems in turbulent underwater environments. A general in situ calibration procedure is proposed which allows grate improvements in terms of fish detection, identification and tracking capabilities.

I. INTRODUCTION

Generally, in dams or weirs (natural or artificial), hydro power plants have a certain impact on the water ecosystems especially for fish (migratory species) and aquatic organisms. Many fish species migrate on more or less extended areas, such as the wild salmon and eels which migrate several thousand of kilometers during their life cycle. The ecological continuity concept is defined as the possibility of free movement of animals (especially fish) and a free sediment transport. Considering the European regulation (European Water Frame Directive, EEL regulation), ecological continuity is a key criterion for fulfillment of the objectives of the regulation in rivers.

Fish passes and downstream bypass are of increasing importance for the restoration of free passage for fish and other aquatic species in rivers. Such devices are often the only way to make it possible for aquatic fauna to pass obstacles that block their up and down river journey. Efficiency or passability of these new and old structures must be quantified in situ. For this, various methods have been developed and used, such as trapping, telemetry or video counting. All these methods require manual and heavy processes (fish manipulation) or are dependent of the light and water quality.

Firstly, Spampinato et al. introduced in [1] an underwater imaging system for detecting fishes in Ken-Ding, Taiwan. The proposed imaging systems has three subcomponents: texture analysis, fish detection and fish tracking. The main inconvenient of such a system is the fact that it can operate during the day or with artificial illumination, only. Such a constraint is not applicable for wild salmons in monotonous regions since the smolts are mainly migrating during the night and any artificial illumination source will perturb them.

A second method for wild salmon counting is based on resistive measurements [START_REF] Smith | Development of a "weir-less" electrode array for a resistivity fish counter: results of field trials in 1993[END_REF]. Using two electrodes, this monitoring system is measuring the conductivity fluctuations (slight decrease in amplitude) when a fish is passing through. Such techniques are well suited for calm underwater environments such as the fish elevators where no considerable turbulent follow is present. This is not well adapted for downstream smolts migrations in Alpine rivers.

Finally, An et al. recently introduced ultrasound for fish monitoring [START_REF] An | Ultrasound backscatter from free-swimming fish at 1 MHz for fish identification[END_REF]: by using an ultrasound emitter, the acoustic signal scattered by the fish is analyzed using three underwater cameras. In this context, an active ultrasound fish monitoring system can efficiently contribute to detect, quantify and identify in the water column the fish species continuously. This paper illustrates an example of active ultrasound wild salmon monitoring system which can be installed on mountainous rivers. The current system is composed of several acoustic barriers operating at multiple frequencies. According to Fig. 1-(b), each salmon passage will obstruct the acoustic propagation channel between the emitter and the receiver, which makes possible wild salmon monitoring. Each acoustic barrier acts as a conventional immersed wideband ultrasound system in bistatic configuration: the emitted waveform is a linear Hamming windowed chirp. However, after passing through the underwater media, the power and the low noise amplifiers, the emitted waveform is deformed. This happens also due to the frequency varying electric impedance of the ultrasound transducers. As a consequence, the results, after the match filter detection at the reception, are less accurate. In this paper, we introduce an automatic method for reference selection in order to provide a more robust fish detection result.

II. REFERENCE SELECTION PROCEDURE

In many signal processing applications the analyzed signals are nonstationary in the sense that the frequency changes over time and are comprised of multiple components each having its particular instantaneous frequency law (IFL). Most of the times the phase function of an analytic non-stationary component can be described by a parametric model. If there is no a priori information of the signals nature, a widely employed model is the polynomial phase signal (PPS) model based on the Weierstrass approximation theorem.

The proposed approach is based on time warping [START_REF] Shen | Wideband timevarying interference suppression using matched signal transforms[END_REF] the given signal with different parametric functions and identifying for which parameters the warped signal has the best spectral concentration [START_REF] Stankovic | A measure of some time-frequency distributions concentration[END_REF]. The cost function used for evaluating the spectral concentration is derived from the measures designed to evaluate time-frequency distributions concentration. Full detail of the range acoustic autofocusing algorithm [START_REF] Anghel | Modelbased parameter estimation of non-stationary signals using time warping and a measure of spectral concentration[END_REF] is illustrated on Fig. 2. The acoustic autofocusing approach has the main advantage that it does not need a particular reference target to estimate the correction nonlinearities of the acoustic channel and it works in the presence of multiple reflective targets such as fishes. However, the number of model parameters is limited from an operational point of view: it is not practical to minimize a cost function of too many variables. In our case three is enough.

For the ultrasound wild salmon monitoring system under study, the received frequency-voltage characteristic can be approximated by a third order polynomial expression. Applying the proposed algorithm on the received signal, the nonlinearity coefficients are determined and the range profile computed by match filtering with the signal warped in the optimal basis as reference.

The only question remains how distinguish between unperturbed reference signal and the revived signal affected by the fish passages. This is done by computing histograms (empirical probability density function) with the LFI coefficients obtained from different realizations of the received signal for each acoustic barrier. By considering the fish passages as sparse events, we obtain a bi-modal histogram shape on which conventional histogram thresholding methods such as [START_REF] Otsu | A threshold selection method from gray-level histograms[END_REF] can be applied.

III. RESULTS

The experimental setup is composed by one acoustic barrier equipped with Ultran GS-30 ultrasound transducers. The transducer transfer function is shown on Often, multiple acoustic echoes can be encountered in the same range profile. Fig. 6-(a) illustrates a full range profile which presents at least four indirect paths. One question is rising, namely whether we should use these reflections in order to improve the detection result after match filtering or not (the improvement is marginal). The proposed acoustic autofocusing algorithm is used to take into account the first two reflected paths. Indeed, the reflections can be considered as harmonics of the direct path in the proposed approach. The obtained result is presented on Fig. 6-(b). One can observe that the practical detection improvement by considering the first two indirect paths is quite marginal for our study. The smolt detection results are illustrated on Fig. 7. One can notice that the smolts detection can be directly performed after match filtering. Indeed, when a smolt is passing through the acoustic barrier, a decrease appears in the focused received signal. The principle is graphically explained on Fig. 8-(a),(b).

Finally, an interesting phenomena is occurring when the smolt is passing in the proximity of the acoustic barrier (e.g. as in Fig. 8-(c)): due to indirect path reflexion, the focused received signal level is increased.

IV. CONCLUSION

The general automatic reference selection procedure presented in this paper offers an alternative option to monitor wild salmon migration in turbulent underwater environments. The main advantages of the proposed algorithm are improved detection capabilities at the reception, the improved resolution and the opportunity for automatic in situ acoustic channel equalization. Future studies include developing full identification and tracking capabilities of different fish types.

Fig. 1 .

 1 Fig. 1. The principle of the ultrasound active system.

Fig. 2 .

 2 Fig. 2. Range acoustic autofocusing algorithm diagram for a 2-parameter nonlinearity model.

Fig. 3 .

 3 One can notice three different frequency bands with linear phase: around 36 KHz, around 110 KHz and around 460 KHz. The emitted waveform is a linear chirp in the frequency band 400 KHz to 500 KHz. The measured directivity diagram and the static characteristics in all three operating frequency bands of the ultrasound transducers are illustrated on Fig. 4.

Fig. 3 .

 3 Fig. 3. Ultran GS-30 ultrasound transducers transfer function.

Fig. 5 -

 5 Fig. 5-(a) illustrates the obtained histograms of the first and the second dictionary coefficients. One can notice a bimodal empirical probability density function in both cases. Fig. 5-(b) shows the comparison between the range profile obtained by match filtering with the emission (in red) and with the estimated reference (in blue): the detection in this range profile is improved by 30 dB using the proposed reference selection scheme.Often, multiple acoustic echoes can be encountered in the same range profile. Fig.6-(a) illustrates a full range profile which presents at least four indirect paths. One question is rising, namely whether we should use these reflections in order to improve the detection result after match filtering or not (the improvement is marginal). The proposed acoustic autofocusing algorithm is used to take into account the first two reflected paths. Indeed, the reflections can be considered as harmonics of the direct path in the proposed approach. The obtained result is presented on Fig.6-(b). One can observe that the practical detection improvement by considering the first two indirect paths is quite marginal for our study.

Fig. 4 .

 4 Fig. 4. Ultran GS-30 ultrasound transducers: (a) measured directivity diagram, (b) static characteristics in all three frequency bands.

Fig. 6 .

 6 Fig. 6. (a) obtained range profile -direct path and indirect paths, (b) range acoustic autofocusing using the direct path and the first two indirect paths.

Fig. 7 .

 7 Fig. 7. An example of smolt detection result.

Fig. 8 .

 8 Fig. 8. (a) acoustic channel obstruction when a smolt is passing through the barrier, (b) impact of the smolt passage on the focused received signal, (c) acoustic channel perturbation (by reflexion) when a smolt is passing in the proximity of the barrier.


Jean Chandezon
email: jean.chandezon@univ-bpclermont.fr

Jean-Pierre Plumey

New approach of diffraction of electromagnetic waves by a rough surface

We consider the diffraction of a plane wave by a rough surface. In the sake of simplicity the study is restricted to the case of perfectly conducting surfaces. Solving this problem is only possible by limiting the infinite rough surface to a window of width D. We show that we can obtain the diffraction pattern at infinity in the Fraunhofer zone from the modeling of diffraction by a grating with period D whose elementary pattern coincides with the rough surface in the window D. We give some numerical results for triangular profiles or rectified cosine. We show that for small heights we find that the most widely Kirchhoff approximation is very well checked. This modeling can be applied to Fraunhofer diffraction problem by a non-planar metal strip and the complementary problem of diffraction by a perfectly conducting screen, infinitely thin and with a slit of one or more periods.

I. INTRODUCTION

The problem of diffraction of electromagnetic waves by rough surfaces occurs in many areas: optics, radar detection, telecommunications and more recently in computer graphics for the characterization of textures of materials. If diffractive surfaces are perfectly conducting, the equations of the problem are dimensionless and therefore apply regardless of the frequency used.

According to the geometric nature of the surface methods used to solve are of different types depending on the value of the radius of curvature with respect to the wavelength. They appeal to geometrical optics for very large radii of curvature, physical optics or geometrical theory of diffraction for large curvatures and electromagnetism to the curvatures of the order of magnitude of the length wave. Scalar Kirchhoff approximation theory (KA) is the most widely used theory in the wave scattering from rough surfaces [START_REF] Beckmann | The scattering of elecromagnetics waves from rough surfaces[END_REF][START_REF] Born | Principles of optics[END_REF][START_REF] Ogilvy | Theory of wave scattering from random rough surfaces[END_REF]. The reason is that it is easily understandable physical basis and leads to relatively simple analytical expressions for scattered field amplitudes (in some important limits). The basic feature of the theory is the assumption that the wave field on the surface of a scatterer is approximated as follows: each point of the surface is considered as a point of an infinite plane surface parallel to the local tangent surface.The interaction of the wave field with the surface is treated as the interaction of that field with that plane surface. The KA method gives good results in a rather large domain, but when the roughness becomes large, it appears multiple reflections; then it's necessary to use rigorous electromagnetic methods. In this paper, we propose a new approach to this problem in the resonant domain based on the resolution of the diffraction gratings, surfaces for which there are rigorous numerical models and computer gratings codes (GC) very effective. For the sake of simplicity we have chosen to limit our study to perfectly conducting surfaces.

Our results concern the diffraction of a plane wave by a window of width D of a rough surface. They also apply to the case of diffraction by a metal ribbon of width D and the complementary problem (Babinet's theorem) of the diffraction by an infinitely thin metallic grating with a slit of width D. In the case of surfaces with low roughness height we show that one finds well the results given by KA.

In the first part we recall the basic principles of the study of diffraction by gratings and in the second, those of the aperiodic diffraction surfaces. In the third and in the fourth part we propose a new approach to numerical modeling of the diffraction of a plane wave by a rough surface. In the fifth part we give some numerical results which show that the proposed method, although much more general, is in perfect agreement with (KA).

II. DIFFRACTION BY A GRATING

Diffraction of a plane wave by a grating

Let us consider a grating whose surface coincides with a cylindrical surface, described in rectangular Cartesian coordinate system Oxyz by a periodic profile function y = ha(x) with period d (Fig. 1). The grating is illuminated from the vacuum by an incident monochromatic plane wave with wavelength λ under the incident angle θ 0 . The time dependance is ∂/∂t → iω where ω is the angular frequency [START_REF] Born | Principles of optics[END_REF][START_REF] Petit | Electromagnetic Theory of Gratings[END_REF][START_REF] Neviere | Light Propagation in Periodic Media[END_REF]. Although the method is applicable to any grating, for the sake of simplicity, we consider only the case of perfectly conducting grating.

For 2D problems (∂/∂z → 0) there are two kinds of solutions TE and TM. In both cases the electromagnetic field is completely determined by the knowledge of its component along Oz; we note that component F with F = E z for TE mode and F = H z for TM mode. In the classic problem the grating is illuminated by one normalized plane wave, under incidence θ 0 : F -(x, y) = e -ikα0x e -ikβ0y . Above the grating, y > max(ha(x)), the diffracted wave has the form of a Rayleigh development which consists of a sum of plane and evanescent waves

F + (x, y) = +∞ n=-∞ A + n e -ikαnx e -ikβny with α n = α 0 + n λ d . (1)
If |α n | ≤ 1 then A + n correspond to outgoing propagating plane wave with angle θ n so that α n = sin θ n +nλ/d with n ∈ U :

U = n ∈ Z| sin θ 0 + n λ d < 1 . (2)
If |α n | > 1 then A + n correspond to an evanescent wave. With our notations the propagation constant β n according to Oy is given by:

β n = 1 -α 2 n = cos θ n if |α n | ≤ 1 -i α 2 n -1 if |α n | > 1 . (3)
Far from the grating, in the Fresnel zone, evanescent waves vanished, the diffracted field is only composed of a finite set of plane waves with angles θ n

F + (x, y) = n∈U A + n e -ik sin θnx e -ik cos θny . (4)
The boundary conditions on the grating surface are used to calculate the diffracted amplitudes A + n in the directions θ n and therefore efficiencies in these directions:

E n = |A + n | 2 cos θ n / cos θ 0 .
For numerical modeling we use the C-method or Chandezon method [START_REF] Chandezon | A new theoretical method for diffraction gratings and its numerical application[END_REF][START_REF] Chandezon | Multicoated gratings: a differential formalism applicable in the entire optical region[END_REF][START_REF] Li | Improvement of the coordinate transformation method for surface-relief gratings with sharp edges[END_REF][START_REF] Li | A rigorous and efficient grating analysis method made easy for optical engineers[END_REF] developed in our Laboratory from the late 70s. This is a differential method of solving Maxwell's equations in a curvilinear translation coordinates system (u, v) adapted to the geometry of the surface such that the surface v = 0 coincides with the grating surface.

The extended grating problem

In fact, numerical modeling gives the solution of a more general problem called here extended grating problem (EGP), where the incident wave is not a single plane wave, but a discrete sum of plane propagating waves or evanescent waves according α n correlated with each other, with complex amplitudes A - n , whose α n propagation constants obey to grating equation (Eq. ( 2)) . For y > max(ha(x)) the total field is written

F (x, y) = +∞ n=-∞ A - n e -ikαnx e +ikβny + +∞ n=-∞ A + n e -ikαnx e -ikβny . (5)
The first sum represents the incoming waves and the second outgoing waves. In the numerical approximation of order M with -M ≤ n ≤ +M we keep only 2M +1 terms in this development. The writing of boundary conditions on the surface y = ha(x) allows to determine the linear relationship between the 2M + 1 amplitudes of the incoming waves A - n assumed to be known, and those of the 2M + 1 amplitudes of outgoing waves

A + m A + n = +M m=-M S m,n A - m . (6)
The S matrix with elements S m,n is the diffraction matrix associated to the profile for incident direction θ 0 .

If we are only interested by the far field this matrix is restricted to real orders, it is denoted S ∞ with elements S m,n where m, n ∈ U . We can then write

A + n = m∈U S m,n A - m , with: m, n ∈ U. ( 7
)
the knowledge of S ∞ is sufficient to determine the far field diffracted by the grating if we know the incident far field. It is not possible to calculate directly S ∞ but only S. After normalization, efficiencies or energy distribution between the different directions are given by

E n,m = A + n,m A + n,m (8)
Quantity E n,m corresponds to the diffracted energy in the direction θ n when the grating is illuminated by a plane wave in the direction θ m .

Diffraction of a beam by a grating

In the most general problem only a finite portion of the grating is illuminated by the incident wave which is then a beam consisting of a continuous sum of plane waves

F -(x, y) = ˆ+1 α=-1 A -(α)e -ikαx e +ikβy dα. (9)
A plane wave with incidence θ 0 is a particular case of beam for which A -(α) = δ(sin θ 0 -α) where δ(x) is a Dirac function.

The diffracted wave is also a beam including propagating and evanescent waves

F + (x, y) = ˆ+∞ α=-∞ A + (α)e -ikαx e -ikβy dα. ( 10
)
Solving the problem is to determineA + (α) fromA -(α) assumed known. This can be done by discretizing α from an origin α 0 with a regular spacing ∆α where ∆α = ∆(sin θ) is the chosen angular resolution. Calculations are then made with the GC for the EGP for all associated α 0 + n∆α incidences. This rigorous resolution takes into account diffraction of all the grating patterns.

In the following we propose a method of resolution to consider only a single pattern of the grating.

Problem formulation

We consider an aperiodic infinite rough surface Σ (Fig. 2), perfectly conducting, with profile y = ha(x) -y 0 on which apply the boundary conditions between the incident and the diffracted field. The dimensionless function a(x) is normalized such that max(a(x)) -min(a(x)) = 1. The parameter h is the maximum amplitude of the modulation.

The plane surface origin y = 0 is chosen so as to be immediately above the rough surface then y 0 = h max(a(x)).

We are looking for the diffracted field F + (x, y) by Σ illuminated by the incident wave F -(x, y). In Cartesian coordinates the most general solution of Maxwell's equations in the half space y > 0 above Σ is a beam which is written in the integral form

F ± (x, y) = ˆ+∞ α=-∞ A ± (α)e -ikαx e ∓ikβy dα, (11)
with α 2 + β 2 = 1 and where F ± (x, y) is the component of the electromagnetic field along Oz. The ± sign corresponds respectively to outgoing or incoming waves.

If |α| ≤ 1 then A ± (α) corresponds to an outgoing or incoming plane propagating wave according to the θ angle with α = sin θ and if |α| > 1 to an evanescent or anti-evanescent wave whose propagation constant β along Oy is given by Eq. ( 2) . At infinity it remains only the propagating waves for which |α| < 1.

In the surface y = 0, the field is

F ± (x) = ˆ+∞ α=-∞ A ± (α)e -ikαx dα. ( 12
)
Amplitudes A ± (α) are Fourier transform of the field for y = 0

A ± (α) = 1 λ ˆ+∞ x=-∞ F ± (x)e +ikαx dx. (13)
The energy radiated by each of the plane waves with angle θ which compose the diffracted beam can be written

E (θ) = A + (θ)A + (θ). (14)
Solving the problem of diffraction by the rough surface is, writing the boundary conditions, to determine A + (α) knowing A -(α).

Classical methods

Most of the time, the problem of diffraction by an aperiodic surface is modeled by the methods of physical optics like KA [START_REF] Beckmann | The scattering of elecromagnetics waves from rough surfaces[END_REF][START_REF] Born | Principles of optics[END_REF] which is to assume that the surface can, at any point, be represented by its plane tangent surface and calculating the wave diffracted by the plan from the Fresnel coefficients. KA gives very good results for the roughness of low amplitude and low slope ie when there are no multiple reflections.

In the case where the dimension of the roughness is of the order of the wavelength, and where there are multiple reflections, it is necessary to use a much more expensive electromagnetic theory in calculation time. In this domain we find many works [START_REF] Voronovich | Wave Scattering from Rough Surfaces[END_REF][START_REF] Sailllard | Scattering from metallic and dielectric rough surfaces[END_REF][START_REF] Edee | A curvilinear coordinate based hybrid method for scattering of plane wave from rough surfaces[END_REF][START_REF] Dusséaux | Scattering of plane wave by 1dimensional rough surface study in a nonorthogonal coordinate system[END_REF][START_REF] Soto-Crespo | Electromagnetic scattering from very rough random surfaces and deep reflection gratings[END_REF].

Maystre et al [START_REF] Maystre | Scattering by one or two dimensional randomly rough surfaces[END_REF] showed that in their opinion, there are three main methods to solve this problem:

1) to consider that an aperiodic surface is the limit, when the period tends to infinity, of a periodic surface leading to periodize the surface and then solve the problem like a diffraction grating. The grating program gives only a finite and discrete set of diffraction directions angularly separated by ∆ sin θ = λ/D . Therefore, we are led to select a large period to limit the interpolation problems between calculation points. 2) to illuminate by a plane wave over a window of width L of the surface and then solve the problem with a GC considering the incident truncated plane wave as a beam. 3) to illuminate the window L with a Gaussian beam for example, for having an illumination which tends to zero at the edges of the window. In fact this method is a generalization of method 2 to arbitrary beams.

To these three cases, we can add the canonical problem of diffraction of a plane wave by a plan locally deformed where, most of the time, the proposed solutions are type 1. There is an additional difficulty which is the response of the deformation is embedded in the specular reflection of the plane mirror. Then it is necessary to subtract the response of plane mirror to the found solution. The numerical experiment shows that this method is rapidly becoming unusable when the height h increases.

Finally, in all cases, the problem is solved by a periodization then, in addition to the desired solution, this method introduces the influence of multiple diffractions between all patterns of the grating associated with Σ. The influence of these multiple diffractions increases very rapidly with h which limits the field of application of these methods. One way to eliminate the influence of multiple diffractions is to introduce between the patterns an area where we place a fictitious material (PML) eliminating the couplings between the patterns [START_REF] Edee | Complex coordinate implementation in the curvilinear coordinate method: application to plane-wave diffraction by nonperiodic rough surfaces[END_REF].

As for us we propose an interpolation method which easily suprime the multiple diffractions with an optimal modeling of the problem. In the limiting case where the period tends to infinity we find method 1.

IV. NEW APPROACH TO DIFFRACTION BY AN APERIODIC SURFACE

Discrete approximation of the electromagnetic field

Except for trivial case h = 0, where the surface is a plane mirror there is no analytical solution; then we are led to replace the original continuous problem in α by a discrete problem by discretizing α and then amplitudes A ± (α) of incident and diffracted fields. By assumption, samples A n = A(α n ) are taken regularly according to the variable α posing α n = α 0 + nλ/D where D is a parameter that has the dimension of a length.

Around an arbitrary origin according to x, define a spatial window D (Fig. 2 ) such as D = {x ∈ R | -D/2 < x < +D/2}. Consider the electromagnetic field F 0 (x) = F (x, y = 0) in the plan y = 0, the discrete representation of the field transforms the aperiodic function in a periodic function of period D

F 0 (x) = λ D +∞ n=-∞ B n e -ikαnx , x ∈ R. ( 15
)
On the other hand we put

F D (x) = Π(x/D)F 0 (x), (16)
where Π(x) is the rectangular function

Π(x) = 1 |x| ≤ 1 0 |x| > 1 . (17)
F D (x) is a function that is zero outside of x ∈ D written F 0 (x) = F D (x) + (1 -Π(x/D)) F 0 (x). (18)
The field F D (x), can be represented by a Fourier series

F D (x) = λ D +∞ n=-∞ B n e -ikαnx if x ∈ D F D (x) = 0 if x / ∈ D. ( 19
)

Profile representation

We have shown that after discretization the modeled electromagnetic field is periodic with period D; this implies that the geometry that is involved in numerical modeling is also periodic of period D. This leads to build from an aperiodic surface Σ a periodic associated surface Σ which coincides with Σ in D. It is this surface Σ with elementary pattern Σ D for x ∈ D which is taken into account in the modeling. Therefore, the numerical modeling is totally identical for all surfaces Σ coinciding in D leading to the same surface Σ . This applies in particular to the following three cases:

• diffraction by a locally deformed plan surface defined by y = f (x) = Π(x/D)ha(x) (outside D the surface Σ is a plan), • diffraction by a metallic ribbon of width D where there is no diffraction outside D, • diffraction by only one pattern of a periodic surface, ie a grating, so that y = f (x + D) = f (x). As we have rigorous programs GC for the diffraction of a plane wave by gratings ie for periodic surfaces Σ . It is from GC that we will deduct the diffraction in the window D of the rough surface Σ. This model gives the diffracted energy only for λ/D. If the GC is used to calculate additional points for truncated plane wave, additional rigorous results are valid for the infinite grating and do not affect the problem for a single patternΣ D . If, for example, calculations are made to values such that ∆α = λ/2D the new results concern a 2D wide window that takes into account the multiple diffractions between two consecutive patterns Σ D of Σ .

Above the associated grating to Σ, for y > 0, the electromagnetic field is written in the form of Rayleigh development Eq. ( 5). For y = 0 coefficients A + n of Rayleigh development given by the GC coincide with B n Eq. ( 19). If we know all A + n then we know, in the approximation λ/D the diffracted field at any point of D for y > 0.

The different domains for the electromagnetic field

The modeling gives the value of the field in D. As in antenna theory [START_REF] Ishimaru | Electromagnetic Wave Propagation, Radiation, and Scattering[END_REF] according to the value of y, different domains can be defined for the electromagnetic field:

• the field inside roughness for ha(x) < y < 0 which is expressed analytically in terms of curvilinear translation coordinates (u, v) of C-Method, • the near field for 0 ≤ y 3λ which is expressed in Cartesian coordinates system as the sum of plane and evanescent waves, • the far field for 3λ y D 2 /λ (Fresnel domain) which is expressed as a finite sum of plane waves which interfere with each other,

• the field at infinity for y D 2 /λ (Fraunhofer domain ) where there is no interference at all, in any point the wave is locally a plane wave. In Fraunhofer domain, at a point P located at the distance r from the origin, in the vicinity of normal incidence the width D of the diffractive element is seen under the angle φ D/r. If φ λ/D ie r D 2 /λ then the diffracting element is seen as a point in the retaining approximation. As a result, the diffracted wave by this element is locally plane in the vicinity of P and it propagates in the direction θ, it is then written

rF ( - → r ) = A(θ)e -i - → k . - → r si r D 2 /λ. ( 20
)
This expression is valid for any point at the distance r from the origin including outside the window D, it is an extrapolation relationship for the field. As it is usual in the antenna theory, for a given r , all the points of the diffracting object being at the same distance r, amplitudes A(θ) are in phase. On the other hand due to the energy conservation field amplitude decreases as 1/r.

V. NUMERICAL MODELING

In all that follows the incident field consists of a single plane wave with incidence angle θ 0 . We choose an angular resolution ∆ sin θ = λ/D, numerical modeling is made from our GC for a grating with period D, whose elementary pattern Σ D coincides with Σ in D. In Fig. 3 the surface Σ is a plan surface with a triangular deformation. We solve the EGP by writing the boundary conditions on the periodic surface y = ha(x).

Incident field

For y = 0 the incident plane wave F - 0 (x) = e -ik sin θ0x is divided into two parts: the fraction that illuminates a portion of width L i of D centered on the abscissa x = l i and the rest of the plane wave

F - 0 (x) = Π( x -l i L i )F - 0 (x) + F - 0 (x) -Π( x -l i L i )F - 0 (x) . ( 21
)
The first term represents a pencil of light of width L i , the second term represents the complete plane wave that illuminates the entire grating except the window of width L i (Babinet's theorem [START_REF] Born | Principles of optics[END_REF]). Each of the two terms of Eq. ( 21) forms a beam. Maxwell's equations are linear, the solution to the diffraction of the plane wave is the sum of solutions obtained for each of the terms. What concerns us here is the solution for the first term, which represents response of D to excitation by a width of light beam L i centered at l i in the plane y = 0.

In the Fourier space the first term is:

F - Li (x) = Π( x -l i L i )e -ik sin θ0x (22) = ˆ+∞ α=-∞ A -(α)e -ikαx dα, (23)
with

A -(α) = L i λ sinc(π(α -α 0 )L i /λ)e +i2πli/D , (24)
where sinc(x) = sin(x)/x is the cardinal sine. With the chosen angular resolution using the discretized amplitudes of the incident wave

A - n = A -(α n ) are A - n = L i λ sinc(πnL i /D)e +i2πli/D . ( 25
)
If L i = D all A - n are equal zero except for n = 0

Diffracted field by D

For periodic problem associated with D, above the grating (y > 0 ) the diffracted field is written in the form of a Rayleigh development [START_REF] Petit | Electromagnetic Theory of Gratings[END_REF]. As we are concerned only by the field to infinity we restrict this development only to plane waves

F + ∞ (x, y) = n∈U A + n e -ikαnx e -ikβny . ( 26
)
In the plane y = 0, each of the plane waves of this development can be considered as the sum of a truncated plane wave F + (D),n of width D and the rest of the plane wave. For each diffraction order, we can write

F + (D),n (x) = Π( x D )A + n e -ikαnx , (27)
let, in the Fourier space

F + (D),n (x) = ˆ+∞ α=-∞ A + n (α)e -ikαx dα, ( 28
)
with

A + n (α) = A + n D λ sinc (π(α -α n )D/λ) . ( 29
)
A priori modeling gives the electromagnetic field only in the window D however, to a point at infinity in the Fraunhofer domain at the distance r of origin the diffracted field by D can be considered as the radiation of a diffracting point located at the origin O; we can then think that we may use the above equation to calculate the field by extrapolation to infinity in the form of a locally plane wave Eq. ( 20)). In the Fraunhofer domain the amplitude of the plane wave at a point P is the sum of all contributions of each order diffracted by the grating

A + (α) = n∈U A + n D λ sinc (π(α -α n )D/λ) , (30) let A
+ (θ) = n∈U A + n D λ sinc (π(sin θ -sin θ n )D/λ) . (31)
At the angular precision related to D, in Fraunhofer domain, for a given value of r the energy density radiated in the direction θ is then

E (θ) = A + (θ)A + (θ). (32)
We note that if θ = θ n then A + (θ n ) = A + n , diffracted energy in such direction θ n is equal to the efficiency of the associated grating.

VI. RESULTS

Results are given for an aperiodic surface, we consider only the window D. This surface is illuminated by a plane wave of incidence θ 0 , we compute the diffracted wave at infinity created by only the part of the incident plane wave that illuminates in a plane y = 0 the part of width L i of D centered on l i (Fig. 3).

Principle of calculation

In summary, from the GC the resolution of diffraction problem by a part of window D of the rough surfaceΣ is done in five steps:

1) determination of the diffraction matrix S ∞ (θ 0 ) of periodized problem with elementary pattern Σ D for the angular resolution ∆ sin θ = λ/D with incidence θ 0 , 2) computation of A - n for the incident plane wave truncated to width L i centered on a point located at the distance l i of the origin Eq. (24).

3) computation of A + n Eq. ( 7) of the diffracted wave to infinity by the grating of period D , 4) calculation of the analytical expression of the diffracted amplitude A + (θ) by interpolation in cardinal sine from A + n computed in step 3, Eq. (30), 5) calculation of the angular distribution of diffracted energy E (θ) normalized with respect to energy of the incident beam where angles are given in degrees.

Diffraction by a deformed plane surface Deformation of small height in Kirchhoff method domain To be able to simply compare our results with those that would be obtained in the Kirchhoff method the chosen deformation ABCDE consists of four line segments Fig. 4. To limit the number of parameters and simplify the interpretation of results we chose AB=BC'=C'D=DE=b = 3.2λ and a normal incidence θ 0 = 0°. The geometrical optics shows that up to a slope ϕ = 45°there is no double reflection on the entire pattern, this is the Kirchhoff domain where we have chosen the first example with h = 1.848λ and ϕ = 30°Fig. 5. this grating equal to an integer number of periods. The selected operating point is in the resonant domain : d = 2.3λ, h = 1.1λ and θ 0 = 8°. For this grating there are four real diffraction orders (-2,-1,0,+1). Efficiencies for this grating are given in Table 1

Fraunhofer diffraction by a metallic ribbon

We have chosen to determine the diffraction of a portion of width D of a perfectly conductive surface, the results relate only diffraction by this surface portion regardless of what exists outside this window. Hence the Fraunhofer diffraction pattern obtained is the same that for an infinitely thin metal ribbon of width D having the form of a single grating pattern. By Babinet's theorem is also diffraction obtained for the complementary problem with a perfectly conducting periodic screen where we practiced a slot of width of one pattern, the diffraction pattern is then the same as that of the ribbon except for the angles θ n where we find efficiencies for the infinite grating. In both cases of geometry discussed in this paper results for the ribbon can be found in Figs. 5a and9a.

CONCLUSION

Currently, there are rigorous numerical modeling to determine the diffraction of a plane wave by a grating. One may think that the results obtained for a grating formed by the infinite repetition of an elementary pattern contain the results of diffraction of the plane wave by only one of those patterns. We adopted this approach to investigate the diffraction of a plane wave by a window of width D of a rough surface considered as single pattern of a periodic surface with period D. Our numerical results show that one finds well the expected results that are given by geometrical optics, the ray tracing and optical physics.

The results of Fig. 5 concerning a weakly modulated plane surface are in perfect agreement with those provided by KA which assumes that surface is plane in the vicinity of each point. When the height of roughness increases, multiple reflections appear and when both sides of the pattern are at right angles we find the phenomena of retro reflection provided by geometrical optics, this phenomenon persists in the resonant domain as can be seen in Fig. 7.

We used the linearity of Maxwell's equations to decompose the incident plane wave and the diffracted plane waves as the sum of a truncated plane wave and its complement. This is in perfect accord with the Huygens Fresnel principle of assuming that at each point the wavefront can be considered as the limit of isotropic radiation point source. If in our modeling we decompose the incident plane wave into a sum of truncated plane waves with very small width L i and if then searches the diffracted wave as the sum of all these contributions then we found the Huygens Fresnel principle.

A final argument in favor of our approach lies in the Heisenberg uncertainty principle [START_REF] Heisenberg | The Physical Principles of the Quantum Theory[END_REF]. If we are interested in the corpuscular aspect of light it is possible to interpret the efficiency of a grating such as the probability of finding in the direction θ n an incident photon associated to a plane wave arriving under the incidence θ 0 . If one restricts the study of the diffracted wave only to incident photons which interact with a single period D of the grating they are diffracted in a direction θ with a probability density A(θ)A(θ). There is then a information on the position of the diffracted photons which come from a point of the pattern (-D/2 < x < D/2) which, by virtue of the uncertainty principle, induces an uncertainty in their direction. In the plan y = 0 the uncertainty relation is written ∆p x .∆x ≥ h where h is Planck's constant, ∆p x the uncertainty of the momentum of the photon according Ox and ∆x = D the uncertainty of the position. For a photon momentum is written p = hν/c = h/λ with ∆p x = ∆(p sin θ) = (h/λ)∆ sin θ which leads to an angular uncertainty ∆ sin θ ≥ λ/D. From the point of view of the uncertainty principle, once fixed the width D, to make additional measures for an angular resolution ∆(sin θ) less than λ/D adds no additional information. It is therefore possible to determine A(θ) only from measures A n = A(θ n ) made for angles obeying the grating formula sin θ n = sin θ 0 + nλ/D, where θ 0 is an arbitrary origin angle.

If the uncertainty principle should be applied to mea-surements, it also applies to numerical modeling. If we are looking for the Fraunhofer diffraction pattern of a window of width D illuminated by a plane wave, all information that is possible to know is contained in the diffracted amplitudes in directions θ n of associated grating. Make additional calculations provides no additional information for diffraction by the window of width D. If, for example, calculations are made to an angular resolution ∆ sin θ = λ/2D < λ/D new results obtained are only valid for L ≥ 2D, they contain the contributions of multiple reflections between two patterns which is not consistent to the originally posed problem. If the roughness are of very low amplitude effect of multiple reflections is negligible then, in this case, the new calculated points correspond to the problem initially posed however, the new calculations are completely useless.

For sake of simplicity we have chosen a perfectly conducting rough surface, currently, taking in account the conductivity is not a problem in particular in C-Method. The extension of the formalism to the finite conductivity should be done without any difficulty.

Figure 1 .

 1 Figure 1. Diffraction by a grating

Figure 2 .

 2 Figure 2. Diffraction by a rough surface

Figure 3 .

 3 Figure 3. Diffraction of a beam by a grating

Figure 4 .

 4 Figure 4. A plan with a triangular deformation

Figure 5 .

 5 Figure 5. Diffraction for small height h= 1.848λ, (a) diffraction by the complete pattern, with dotted line diffraction for h = 0, (b) diffraction by AB, (c) diffraction by BC and (d) diffraction by BD In this example the two polarizations give practically the same results. We represent the reflection from the plane surface D = AE dotted line corresponding to h = 0 and full line corresponding to the diffraction for both TE and TM for h = 1.848λ.

Figure 7 .

 7 Figure 7. Diffraction of a triangular groove for several values of ϕ, when only inclined face CD is illuminated, dotted line TE mode, solid line TM mode.

Figure 8 .

 8 Figure 8. Grating in rectified cosine

Figure 9 .

 9 Figure 9. Diffraction of a plane wave by a planar finite surface of width D = N d with N =1,3,5,9 identical patterns of width d, doted line TE mode, solid line TM mode

Figure 10 .

 10 Figure 10. Diffraction of a truncated plane wave by a surface of width D = N d with N = 3, 5, 9, 11 identical patterns of width d where only the central pattern is illuminated, doted line TE mode, solid line TM mode. Results for N = 1 are given in Fig. 9a.

  . We chose this operating point in the resonant domain to have very different results in TE and TM.

	orders	angles	efficiencies TE efficiencies TM
	-2	-46.92°0.0491	0.2205
	-1	-17.19°0.5900	0.0535
	0	8.00°0.2716	0.0703
	1	35.03°0.0894	0.6557
			Table I
	EFFICIENCIES FOR GRATING WITH PROFILE IN RECTIFIED COSINE

wave (L i = AE, l i = 0) . We note the presence of three maximum : θ = 0 for the reflections on AB and DE and two maximum near θ = ±60°for reflection on inclined plane surfaces CD and BC.

• Fig. 5b: only side AB is illuminated (L i = b, l i = -3/2b), specular reflection is almost perfect. • Fig. 5c: only the inclined side BC is illuminated (L i = b, l i = -b/2), there again it almost perfectly specular reflection for incidence of 30 °which give θ = 60°. • In Fig. 5d: both sides BC and CD are illuminated (L i = 2b, l i = 0), once again reflections are around θ = ±60°T hese results show that, as expected, KA is well verified, replacing the profile by four planar segments led to very good results.

Deformation with multiple reflections We keep the same geometry but with a hollow deformation Fig. 6 that will lead to multiple reflections. The geometrical optics shows that the double reflection appears from ϕ = 30°. For ϕ = 45°after a double reflection, light comes back into the direction of incidence. This phenomenon occurs again at ϕ = 60°after a triple reflection. The numerical experiment shows that the horizontal faces AB and DE on which there is no double reflection always reflect as a horizontal plane mirror. It is for this reason that we have represented only the reflection on the inclined face CD (L i = b, l i = -b/2) for ϕ = 30°, 45°, 60°, 75°.

• Fig. 7a: ϕ = 30 • , it is observed that the result is practically the same as in Fig. 5c, there is little coupling with the other faces. • Fig. 7b :ϕ = 45°, we observe that, with the double reflection, as expected, most of the light comes back into the direction of incidence θ = 0 which is still the case for ϕ = 60°Fig. 7c. • Fig. 7d: ϕ = 75°the effects of reflections triple and quadruple completely distort the specular reflection by CD and the effects of polarization have a important effect. In the last three cases the method Kirchhoff is unable to correctly model the diffraction.

Diffraction by one or more patterns of a grating

The rough surface Fig. 8, is a grating with period d with a(x) = -|cos(2πx/d)|. We consider a part D = N d of


Mihal Brumbulli
email: mihal.brumbulli@pragmadev.com

Emmanuel Gaudin Pragmadev

Automatic Interleaving for Testing Distributed Systems

Keywords: Distributed systems, testing, modeling, simulation, TTCN-3, SDL

The constantly ever-growing interest for large-scale distributed systems like the Internet of Things imposes many challenges for developers and researchers from many areas. The development of distributed software applications is by no means trivial, and their inherent complexity becomes apparent during testing. Indeed, testing the operation of single isolated nodes does not suffice, because it may be affected by the distribution and inter-communication between nodes. Re-writing a test case to consider distribution is neither efficient nor simple, because concurrency is never easy to implement. In this paper we present an approach that automatically interleaves execution of test cases to simulate concurrency inherent in distribution. We focus on independent test cases that might exhibit a correlation due to distributed interaction. The approach is applied in the context of standard modeling and testing languages, and enables identification of interaction points during test case execution that depend on distribution. The re-execution of the test case is then interleaved at the identified points to account for distribution.

I. INTRODUCTION

The development of distributed systems has gained much attention from industry and researchers with a variety of applications, a trend that is sure to continue in the immediate future due to the ever-growing interest for the Internet of Things [START_REF] Inc | Gartner says the Internet of Things installed base will grow to 26 billion units by 2020[END_REF]. However, the development of software applications for large-scale distributed systems is not a trivial task. This is especially true for testing, an activity which is quite challenging even for simple non-distributed systems. The operation of nodes in a distributed system is not isolated, and as such it requires for test cases to account for the distribution and interaction between nodes. The implication here is that existing test cases and their execution have to be adapted to consider distribution. This adaptation consists of (a) introducing concurrency handling into test cases, and (b) controlled concurrent execution that deals with all relevant interleavings. For the former the tests cases have to be modified, and considering that concurrency handling is never easy to implement, the effort that is required should not be overlooked. The later implies the existence of a scheduler that is able to handle all relevant interleavings.

In this paper we present an approach for automating the interleaved execution of test cases. We apply the approach in the context of standard modeling and testing languages. The system under test is described in SDL [START_REF] Itu-T | Specification and Description Language -Overview of SDL-2010[END_REF], TTCN-3 [START_REF] Etsi | Methods for Testing and Specification (MTS); The Testing and Test Control Notation version 3; Part 1: TTCN-3 Core Language[END_REF] is used for testing, and SDL-RT deployment diagrams [4] describe the distribution of system components. Test cases are executed against the system in a simulated environment extended with an interleaving algorithm.

In Section II we give an overview of related work and position our approach in respect to existing state of the art. We introduce the relevant technologies in Section III and our approach in Section IV. An example is given in Section V to illustrate the use of the presented solution. Finally, we conclude in Section VI with a discussion around the approach, its current status, and future work.

II. RELATED WORK

Testing of distributed communication systems has been approached from different angles. We identify two major groups, i.e., distributed testing and interleaved execution.

Hartman et al. in [START_REF] Hartman | A Test Execution Environment Running Abstract Tests for Distributed Software[END_REF] discuss the execution of abstract tests for distributed software. They underline the importance and benefits of interleaving test case sequences for discovering defects in the system. The approach can be applied to a single test case, it uses either concurrent or sequential synchronization of execution sequences, and there is no automation involved in establishing the synchronization points.

Schieferdecker and Vassiliou-Gioles in [START_REF] Schieferdecker | Realizing Distributed TTCN-3 Test Systems with TCI[END_REF] discuss the distribution of TTCN-3 test setups. The focus is on the underlying concepts of the language and how they support management and execution of distributed test cases. This approach implies test execution on the target. This is obviously an advantage, however, it also needs a complex synchronization mechanism for controlling execution of each of the distributed test cases.

Bloom et al. in [START_REF] Blom | TTCN-3 for Distributed Testing Embedded Software[END_REF] emphasize the fact that, prior to testing in a target environment, the software is usually tested in the host environment. They propose a simulation-based approach and focus on the semantics of time. However, it is not clear how the execution of several test instances is interleaved to account for the inherent concurrency.

Testing of concurrent software is discussed in [START_REF] Musuvathi | CHESS: Systematic Stress Testing of Concurrent Software[END_REF] and [START_REF]Iterative Context Bounding for Systematic Testing of Multithreaded Programs[END_REF]. The authors focus on multi-threaded software and propose a solution that interleaves execution in a controlled way. The approach uses neither a standard language (like TTCN-3) nor any kind of abstract notations, instead, it is based on general purpose programming languages.

We adopt the idea of controlled interleaved execution in the context of distributed systems. To do so, we consider distributed communication between instances of the system to be similar to thread interleavings, and can potentially impact the behavior induced by the execution of test cases. Furthermore, the approach is based on standard and formal languages with precise semantics. This allows simulation in a controlled environment and, what is of great interest, automatic generation and execution of all relevant interleavings. The Specification and Description Language (SDL) is a specification language defined by the International Telecommunication Union (ITU-T) in the Z.100 series [START_REF] Itu-T | Specification and Description Language -Overview of SDL-2010[END_REF]. SDL is targeted at the unambiguous specification and description of the behavior of reactive and distributed systems.

1) Structure, Behavior, and Communication: In SDL the overall design is called the system, and everything outside of it is defined as the environment. The system can be composed of agents and communication constructs. There are two kinds of agents: blocks and processes. Blocks can be composed of other agents and communication constructs. When the system is decomposed down to the simplest block, the way the block fulfills its functionality is described with processes. A process provides this functionality via extended finite state machines. It has an implicit queue for signals. A signal has a name and parameters; they go through channels that connect agents and end up in the processes implicit queues. Fig. 1 illustrates these concepts in a simple client-server application example. The client is started with mStart signal which takes as parameter the number of request the client should send to the server. The client will try to connect to the server, and in case of success it will start sending requests and waiting for replies. When done (or in case of error) it will report back (to the environment) the number of replies received from the server via the mDone signal.

2) Deployment: SDL descriptions are platform independent, i.e., they do not capture any information concerning the implementation details. For example, Fig. 1 speaks of a bClient and bServer, however it does not specify whether these agents are distributed or not. Deployment diagrams as defined in [4] can supplement the models with missing information about distribution. This approach has been used for simulating distributed applications as described in [START_REF] Brumbulli | Simulation Configuration Modeling of Distributed Communication Systems[END_REF] and [START_REF] Brumbulli | Model-Driven Development and Simulation of Distributed Communication Systems[END_REF]. We adopt the idea and simplify it by focusing only on components, nodes, and connections as shown in Fig. 2. 1 The semantics are straightforward, i.e., if components (representing SDL agents) are attached to different nodes, then they are considered distributed, otherwise local execution is implied.

B. TTCN-3

The Testing and Test Control Notation Version 3 (TTCN-3) is a standardized testing technology developed and maintained by the European Telecommunication Standards Institute (ETSI). The ETSI TTCN-3 standards [START_REF] Etsi | Methods for Testing and Specification (MTS); The Testing and Test Control Notation version 3; Part 1: TTCN-3 Core Language[END_REF] have also been adopted by the ITU-T in the Z.160 series [START_REF]Testing and Test Control Notation version 3: TTCN-3 core language[END_REF].

The abstract definition of test cases makes it possible to specify a non-proprietary test systems which are independent of both platform and operating system. The abstract definitions can be either compiled or interpreted for execution.

Fig. 3 shows a TTCN-3 module definition with a single test case that triggers the sending of 10 requests (line 15 and 19) from the client and expects 10 replies (line 16 and 21).

IV. APPROACH

A. Problem Statement

We are interested in the effects (if any) induced by distributed execution of test cases. To better understand the problem let us revisit the client-server example, and suppose that the result of the test case presented in Fig. 3 is pass, i.e., the system behaves as expected. What will happen if the number of clients is increased; will the system behave the same? There is one effective way to answer this question: test the system with multiple clients. This can be tackled by (a) rewriting the test case so that it accounts for multiple clients, or (b) execute multiple instances of the test case (one for each client) in parallel. The later needs an underlying synchronization mechanism that allows controlled execution of parallel (distributed) test cases, which is never trivial and requires additional expertise (not TTCN-3). That is why the former is more sound from a tester's perspective, because it is confined at the abstract level of testing provided by TTCN-3. However, rewriting the test case to take into account only the number of clients is not enough. Indeed, the difficult part is not in the number of clients but in their parallel (concurrent) execution due to distribution. To model concurrency all possible interleavings between test cases should be considered, which implies execution of all permutations of TTCN-3 instructions of different clients. For example, if startMessage is sent first to client 1 and then to client 2, it is important to consider also the case where it is first sent to client 2 and then to client 1.

B. Problem Analysis

Supposing concurrent execution of K test cases, with each test case consisting of n i instructions for i = 1, 2, . . . , K, the number of all interleavings is given by:

I = ( K i=1 n i )! K i=1 (n i !) (1)
If we consider concurrent execution of K instances of the same test case, then n i = N ∀i, where N is the number of instructions in the test case, and (1) can be re-written as:

I = (KN )! (N !) K (2)
This is a typical case of the state-explosion problem which makes execution of all interleavings unpractical even in a controlled and automated environment. However, not all interleavings are relevant and their number (in most cases) can be drastically reduced. Indeed, the behavior induced by a test case can be affected by distribution only if there is an interaction between nodes. This means that, if the execution of a test case does not involve any distributed communication, then distribution will not have any impact. For example, if there isn't any communication between the client and the server, then interleaving is pointless because there is no distributed behavior in the first place. That is why it makes sense to interleave execution at critical points, i.e., instructions that trigger interaction between nodes by means of distributed communication.

C. Interleaving Algorithm

We start by grouping the instructions and then interleaving the execution of the groups. The condition is that each group must include at most one instruction which triggers distributed communication. Let m j i be an instruction in the test case, where i = 1, 2, . . . , N is the index (relative order) of the instruction, and j = 0, 1 describes whether the given instruction triggers any interaction (interleaving point). A group consists of all subsequent m j i for which j ≤ 1. The following shows an example sequence and corresponding grouping:

m 0 1 , m 1 2 , m 0 3 g1 , m 1 4 g2 , m 1 5 , m 0 6 g3 , m 1 7 , m 0 8 , m 0 9 g4 , m 1 10 g5 (3)
For two test cases with the above sequence, the algorithm can generate > 700 times less interleavings. 2Interleaved execution can be automated using the following algorithm:

{K is the number of instances} {N is the number of groups} {I is the next interleaving} for i := 0 to K -1 do for j := 0 to N -1 do I[i * N + j] := i end for end for loop handleInterleave(I) i := I.length -1 while i > 0 and

I[i -1] ≥ I[i] do i := i -1 end while if i = 0 then return end if j := I.length -1 while I[j] ≤ I[i -1] do j := j -1 end while temp := I[i -1] I[i -1] := I[j] I[j] := temp j := I.length -1 while i < j do temp := I[i] I[i] := I[j] I[j] := temp i := i + 1 j := j -1 end while end loop
Each group of instructions is represented by a simple integer for ease of calculation. The handleInterleave procedure maps the integer to its corresponding sequence of TTCN-3 instructions (permutation) that is to be executed next.

D. Tool Support

PragmaDev Studio3 is a set of tools that helps specifiers, developers, and testers to manage complexity in the development of today's systems. The tools use the recognized international standards of SDL, TTCN-3, SDL-RT, and UML [START_REF] Omg | OMG Unified Modeling Language (OMG UML). Version 2.5[END_REF].

A key functionality of the tool-set is provided by the PragmaDev (Co-) Simulator as shown in Fig. 4. The co-simulator allows execution of TTCN-3 test cases against an SDL system. SDL and TTCN-3 descriptions are translated into an internal representation (byte code) to be interpreted by the executor, which in turn forwards the scheduling of events to the scheduler. We extend the functionality of the scheduler with the interleaving algorithm.

In the first phase the test case is executed against the system using the scheduler in normal mode (no interleaving involved). All TTCN-3 instructions (send statements) that trigger distributed communication between agents of the SDL system are marked during execution. Every communication between agents is checked against a deployment diagram, and if the sender and receiver of a message are attached to different nodes in the diagram, then the last TTCN-3 instruction that triggered such behavior is indeed the one to be marked.

In the second phase the test case is executed against the system in interleaving mode. The scheduler automatically creates K instances4 of the test case and enters the interleaving algorithm. On each iteration the algorithm creates an interleaved sequence of TTCN-3 instructions based on marking done in the first phase and the number of instances. The sequence is then executed like a normal TTCN-3 test case by the scheduler, and at the end the SDL system is reset to its initial state for the next iteration.

The whole process is completely automated and transparent to the tester. There is no need to rewrite the tests, but just execute them with the interleaving scheduler.

V. EXAMPLE

We have used the presented approach in testing an accesscontrol system. The system is composed of several terminals and a central unit. Each terminal has a slot for entering a card and a keypad for entering the key. This information is sent to the central unit which checks whether access should be granted to a user and notifies the terminal from where the request was issued. The user can be either an administrator or a normal one. The administrator can add or delete normal users and is identified by a special card key.

A. Structure and Behavior

Fig. 5 shows an excerpt of the SDL model of the accesscontrol system. The system is composed of two blocks, namely bLocal for terminals and bCentral for the central unit. Each block has a single process within which implements the behavior of the system. The user enters the card and the key which are sent to the central unit via the checkCardAndCode. If the card and key are those of the administrator, then the central unit enters in administration state, in which normal users can be added or deleted. On the other hand, if the credentials are not those of the administrator, the list of registered normal users is scanned for matching credentials. The employee signal is sent back in case of a match, otherwise an intruder is signaled.

B. Deployment

A simple deployment scenario for the access-control system is shown in Fig. 6. The figure speaks of two bLocal (terminals) connected to a single bCentral (central unit). We will use this scenario to automatically generate and execute all relevant interleavings.

C. Test Case

To show the applicability of automatic interleaving we start with the most basic test case for the system, i.e., try to get in and out of (without doing anything else) administrator mode as shown in Fig. 7. The user enters the card = 0 (administrator) and then the code = 007. At this point the information (card + code) is sent to the central unit. Because these are the right credentials, the user is allowed to enter in administrator mode, where he/she can add or delete users. However, none of this actions is taken, and the user just exits from the administrator mode. The request to exit is also processed by the central unit. Executing this scenario against the system using the scheduler in normal mode (no interleaving, one terminal and one central unit) will indeed result in a pass for the test case.

When execution in normal mode is finished all the required information will be available for entering the interleaving mode. Based on the description above there should have been identified two interleaving points by now: (1) after entering the last digit of the code and (2) after the request to leave the administrator mode. Indeed, these are the points during execution where a communication between the terminal bLocal and the central unit bCentral is triggered. This translates into two groups of TTCN-3 statements whose execution shall be interleaved: the first group consists of lines 34-44 and the second of lines 46-47 in Fig. 7. Based on the number of terminals bLocal in Fig. 6 and that of groups (K = 2, N = 2), the total number of interleavings to be executed is 6. It is important to note that, if the grouping algorithm is not applied, the number of interleavings will be 252. 5What we found during execution in interleaving mode was in fact quite surprising. We didn't expect to find any problems from such a simple test case, considering that it induces a very basic behavior to the system. The system was modeled so that only one terminal at a time can get administrator access, meaning that an attempt from a second terminal will fail. Indeed, this is the behavior we observed during the interleaved execution, however, what we did not expect is for the second terminal to block indefinitely waiting for a reply from the central unit (i.e., test case execution did block on line 44 in Fig. 6). We were able to immediately jump to the point in the model (using the PragmaDev Studio user interface) causing the problem. The central unit, after granting access to the first terminal, entered the administration state so that it could process any requests coming from the terminal for adding or deleting users. However, in this state it was discarding (not handling) all log-in requests, which in turn caused the second terminal to wait indefinitely for a reply. The solution was straightforward: add a new transition in the SDL state machine of the central unit to handle such requests.

We were able to identify 4 similar problems in our model of the access-control system by using the interleaving scheduler on a set of more complex test cases. The set was composed of (a) previously hand written test cases depicting typical usage scenarios and (b) a set of automatically generated test cases using the model-based approach described in [START_REF] Deltour | Model-Based Testing: An Approach with SDL/RTDS and DIVERSITY[END_REF]. It is important to note that we were able to identify these problems without writing a single line of TTCN-3 code.

VI. CONCLUSION

Testing of software applications for large-scale distributed systems is not a trivial task, because the test cases and their execution need to be adapted in order to account for the distribution and interaction between nodes.

In this paper we presented an approach for automating the interleaved execution of test cases to mitigate the complexity of concurrent behavior due to distribution. We applied the approach in the context of standard modeling and testing languages where the system under test was described in SDL, the distribution of system components in deployment diagrams, and the test cases in TTCN-3. Test cases were executed against the system in a simulated and controlled environment with an interleaving scheduler.

The presented algorithm allowed us to significantly reduce the number of interleavings without deviating from our target, that is the thorough testing of the system. Instead of using each single instruction of the test case as interleaving point, we focused only on those relevant, i.e., instructions that induce distributed behavior to the system.

Automatic interleaving was used for testing the distributed behavior of an access-control system. Even with a basic test case we were able to identify a problem in the system, a trend that continued also with more complex test cases. Furthermore, we achieved these results by reducing the number of interleavings by a factor of 42 (in the simplest case).

There are however three issues that need further discussion. First, the proposed algorithm may not always produce significantly less interleavings. This depends on the degree of distribution in the behavior, i.e., the number of interleavings grows if more inter-communication between distributed nodes takes place. This can degrade simulation performance especially with the increasing complexity of test cases. However, we believe this type of scenario to be more of an exception than the rule. Indeed, energy consumption is one of the major challenges of distributed systems composed of potentially millions of battery powered devices (e.g., Internet of Things), and reducing inter-communication between devices to an acceptable minimum is always an engineering goal. Second, although we reported a successful application of the approach with a simple example, we strongly believe that its benefits are emphasized when used for testing complex systems, which is still work in progress as the time of writing this paper. Last, the proposed approach is based on simulation and at present cannot be applied for test cases executing on target. The reason is simple: simulation allows fine-grained control over the execution, does not require any complex synchronization mechanism for interleaving, and can exploit the benefits of the proposed algorithm. The last point is very important because, without any means of reducing the number of interleavings, testing every possible combination is unpractical due to the state-explosion problem. In the end, in addition to a mechanism for controlling the interleaving, execution on target would require means to track distributed communication triggered during testing.

Fig. 1 .

 1 Fig. 1. Excerpt of the SDL model of a client-server application.

Fig. 2 .

 2 Fig. 2. Deployment model of the client-server application.

1 :Fig. 3 .

 13 Fig.3. TTCN-3 module with a test case for the client-server application.

Fig. 4 .

 4 Fig. 4. Architecture of the PragmaDev (Co-) Simulator.

Fig. 5 .

 5 Fig. 5. Excerpt of the SDL model of the access-control system.

Fig. 6 .

 6 Fig. 6. Deployment model of the access-control system.

Fig. 7 .

 7 Fig. 7. Test case for administrator access on the access-control system.

The number of values in the id attribute is the number distributed nodes, e.g., the figure implies two clients and one server.

The result can be obtained by replacing in (2): K = 2 and N = 10 prior to grouping, and K = 2 and N = 5 after grouping.

http://www.pragmadev.com

The number of instances is deduced from the deployment diagram.

In this case the number of groups is equal to that of the TTCN-3 send statements (N = 5).