Yann Clénet

Tristan Buey

Eric Gendron

Zoltan Hubert

Fabrice Vidal

Mathieu Cohen

Frédéric Chapron

Arnaud Sevin

Pierre Fédou

Gaële Barbary

Bruno Borgo

Jean-Michel Huet

Alexandre Blin

Olivier Dupuis

Julien Gaudemard

Saber Ben Nejma

Damien Gratadour

Vincent Deo

Florian Ferreira

Simone Thijs

Vincent Lapeyrère

Jordan Raffard

Fanny Chemla

Bertrand Le Ruyet

Arielle Bertrou-Cantou

Milan Rozel

Youssef Younès

Gérard Rousset

Gérard Zins

Emiliano Diolaiti

Paolo Ciliegi

Vincent Garrel

Sebastian Rabien

Josef Schubert

Michael Hartl

Veronika Hörmann

Richard Davies

MICADO-MAORY SCAO Preliminary design, development plan & calibration strategies

Keywords: ELT, MICADO, MAORY, SCAO, preliminary design, development plan, calibrations

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Stand-alone Relay Optics with MICADO Calibration Unit

INTRODUCTION: MICADO OVERVIEW

MICADO is the ELT first light-imager. [START_REF] Davies | The MICADO first light imager for ELT: overview and operation[END_REF][START_REF] Davies | MICADO sci-ops: shaping its adaptive optics systems[END_REF] Working in the near-IR (0.8-2.4 µm) at the ELT diffraction limit, it will offer four observing modes:

• Standard imaging: with 1.5 & 4 mas pixel scales, the corresponding FoV will be 19 & 51 arcsecond [START_REF] Davies | MICADO sci-ops: shaping its adaptive optics systems[END_REF] . More than 30 broad-band & narrow-band filters will be available

• Astrometric imaging: it drives MICADO design, with a gravity invariant implementation, a fixed mirror optical design, state-of-the-art ADC and dedicated astrometric calibration and data pipeline

• High contrast imaging: it will use the central detector and will be enabled via a classical configuration of focal plane coronagraphs and Lyot stops, as well as pupil plane vAPP coronagraphs and sparse aperture masking. Pupil tracking will be available for angular differential imaging 3

• High resolution spectroscopy: it will provide coverage of a wide wavelength range simultaneously (1.15-1.345 µm, 1.48-2.45 µm or 0.845-1.48 µm) at a resolution of ∼20000 on faint compact or unresolved sources. Three slits will be available: 3 ×16 mas, 15 ×20 mas (for sky subtraction along the slit), 3 ×48 mas (wide)

MICADO will benefit from the Multi-Conjugate Adaptive Optics (MCAO) correction provided by MAORY, 4 for which MICADO's design is optimized. MICADO will also benefit from a SCAO correction developed under MICADO's responsibility and jointly by MICADO and MAORY. This SCAO module will be integrated in MAORY. In case of a delay in the development of MAORY, MICADO will be commissioned at the telescope first-light in a stand-alone configuration with a specific optical relay and with the SCAO module alone (Fig. 1). In any case, in a phased approach of the AO integration at the ELT, SCAO will be the first AO mode to be tested with MICADO at the telescope. MICADO's planning is the following:

• 2018, Nov.: (successful) preliminary design review • 2020, Oct.: final design review • 2024, June: preliminary acceptance Europe

• 2025: preliminary acceptance Chile SCAO first, followed by preliminary acceptance Chile MCAO

MICADO-MAORY SCAO PRELIMINARY DESIGN

SCAO Green Doughnut and dichroic plate

The Green Doughnut (GD) is a volume hosting in its upper part the SCAO module and in its lower part the MAORY MCAO NGS WFS (Fig. 1 and 2 left). The SCAO GD structure is made of 1) the optical bench supporting the SCAO WFS and the SCAO calibration unit, 2) its feet connecting the bench to the MICADO cryostat interface and 3) below the bench, the SCAO dichroic plate assembly (Fig. 2

and 3 left).

The SCAO GD structure has stringent specifications:

• SCAO envelop limited to 2.8m×0.4m

• Mass budget limited to 620 kg • First eigen frequency greater than 50 Hz

• Low CTE to ensure optical alignment over the ESO-specified environmental conditions These specifications are fulfilled thanks to a carbon fiber bench with a baffling and pods contributing to the structural stability.

The SCAO dichroic plate is removable for MCAO operations. Its tilt (18 • ) allows to escape the MCAO NGS patrol FoV with an acceptable vignetting while keeping the induced astigmatism limited. This ∅340mm×30mm dichroic plate has a slight wedge on its rear surface to compensate for the lateral shift of the pupil it introduces at the entrance of the cryostat

SCAO WFS

The SCAO WFS is based on a (double) pyramid WFS. It is sensing the wave-front between 0.589 and 0.96 µm with a ∅2 FoV. Critical functionalities of the WFS are: pupil control & guiding, field patrolling, atmospheric dispersion compensation, NCPA management, wave-front sensing itself. We describe hereafter how these functionalities are addressed in the SCAO design (Fig. 3

right).

Pupil guiding is handled thanks to:

• A pupil imaging lens (PIL) system, located between the WFS camera and the pyramid, in charge of the managing the small pupil offsets

• A pupil steering mirror, located at the entrance of the WFS, in charge of managing large (and rare) pupil offsets Right: Zoom on the WFS components (the pupil imaging camera is not represented).

• A K-mirror stabilizing the pupil in rotation, whose stringent specifications (∼0.1 λ/D shift over 6 • rotation) have motivated specific prototyping 5

A field selector at the entrance of the WFS allows:

• To select, finely position and track the AO target in the SCAO 6 ×20 patrol FoV

• To correct for the differential image position between IR and visible

• To manage dithering/offsets and differential tracking Near an intermediate pupil plane an assembly gathers several components:

• A phase plate assembly, dedicated to the compensation of NCPA between the WFS and the cryostat, supporting three phase plates (one for the dichroic astigmatism, one for this astigmatism and the fixed cryostat low order aberrations, and one neutral)

• An ADC, made of two rotating prisms

• A neutral density positioner

• The modulation mirror, specified at ±0.1 mrad @ 500 Hz and 1 mrad @ 50 Hz

The sensing is performed thanks to:

• An achromatic 20mm×20mm×39mm double pyramid in a F/30 beam

• An ESO in-development ALICE camera based on the E2V CCD220 detector with 240×240 pixels, allowing to image the four pyramid pupils of 96 pixel-diameter each

The WFS also integrates:

• A pupil viewing camera, sampling finely the pupil for AIT & maintenance

• Calibration sources, including for reference slopes, at the entrance of the WFS

SCAO calibration unit

The SCAO calibration unit (SCU) is the SCAO test, calibration and maintenance facility. During its lifetime it will know successively an AIT and an on-sky configuration.

In AIT configuration (Fig. 4 left), the SCU will be used during the AIT phases in Meudon and in the MPE hall in Garching. It will allow to perform the integration and performance tests of the SCAO module. For this purpose, it will be equipped with:

• A newly developed and recently delivered ALPAO 64×64 actuator deformable mirror (DM), currently under tests in Meudon (Fig. 4 right). It will be used in the final AIT phase in Meudon for the validation of the SCAO performance 6, 7

• A slow tip-tilt mirror (TTM) allowing with the 64×64 DM to mimic the ELT M4+M5 configuration

• Turbulent phase screens, pupil masks and calibration/alignment sources SCU will be turned into its final on-sky configuration when MICADO moves to ESO X4 hall in Garching, removing the turbulent phase screens, the pupil masks and the TTM, and replacing the 64×64 DM by a low order DM (so-called XDM) dedicated to NCPA calibration (Fig. 3 left).

SCAO real-time computer

SCAO real-time computer (RTC) design is based on developments made in Observatoire de Paris in the framework of the Green Flash project that are now going on with the COSMIC initiative (gathering Observatoire de Paris, Microgate and the Australian ANU and SUT) that aims at building an RTC platform, modular to adapt to various AO systems 8 SCAO RTC will follow the functional architecture defined by ESO, dividing the RTC into a Hard Real-Time Core (HRTC), a Soft Real-Time Cluster (SRTC) and a Communication Infrastructure.

A classical control law is considered using MVM with up to 25.5k measurements (expanded space control or so-called "full pixel" control 9 ) and 5.4k commands at 500 Hz. Considering the allowed 285 µs latency for computation, it leads to a required 1700 GB/s memory bandwidth in single precision for the HRTC. The maximum data transfer rate for telemetry (considering calibrated pixels, residual modes, DM commands, etc) is 2 Gb/s (in double precision) while the required compute throughput for the supervisor, related to the command matrix optimisation task, performed about every minute, is 314 GFLOP.

The HRTC design is based on hardware (HW) acceleration using NVIDIA GV100 GPU, equipped with HBM2 high memory bandwidth. It will make use of a Microgate FPGA-based HW interface and the GPU direct feature allowing to directly transfer the WFS camera data to the GPU. The HRTC software (SW) stack is divided in 3 levels: an interface module (called Octopus), the core compute engine (called Marlin) and the user interface (called Kraken) 10 (Fig. 5

right).

The SRTC HW will be made of a GPU server while the various functional blocks of the SRTC SW will be assembled around the telemetry middleware to send / receive data from the HRTC module. On top of this, the control middleware is used by the supervisor SW to send commands to the HRTC software to perform observing sequences and optimize operations.

SCAO performance estimation

SCAO performance has been estimated with our GPU-based COMPASS platform that allows fine sampling (e.g. 512×512 pixel support for our SCAO final pyramid image) and fast simulation even at the ELT scale. It is installed on a computer equipped with 2 CPU Intel(R) Xeon(R) CPU E5-2630 v4 @ 2.20GHz (10 cores each) with 64GB RAM and 8 NVIDIA Titan X (Pascal) GPU cards.

For MICADO PDR, we performed AO simulations 11 to • Demonstrate the better performance of the pyramid WFS vs. Shack-Hartmann WFS

• Demonstrate the gain in performance using (Korkiakoski 12 ) pyramid optical gains (though computed with an unrealistic method, with 2 interaction matrices obtained by simulation: one diffraction limited, one under close loop conditions)

• Demonstrate the interest of the "full pixel" centroiding method that virtually cancels any static misalignment of the pyramid WFS pupils and relaxes dramatically the specifications on the pyramid angles

• Optimize the pyramid WFS parameters in terms of -Pupil sampling: performance plateau between ∼88 et ∼100 pixels for the diameter of the four pyramid pupils. Our current baseline is 96 pixel diameter pupils.

-Number of filtered modes vs. guide star magnitude: ∼10 for a 12 magnitude reference source, ∼200 for a 15 magnitude reference source and up to 3000 for a 17 magnitude reference source, under median atmospheric conditions -Modulation amplitude vs. guide star magnitude: current baseline is 5λ/D

• Estimate island effect residuals (differential pistons over the ELT pupil petals), when simply forcing phase continuity over the pupil by slaving the actuators pair-wise on each side of the spiders: 125 nm rms under median atmospheric conditions • Estimate windshake residuals from simulated ELT guide probe measurements: 97 nm rms under median atmospheric conditions

• Validate our NCPA management strategy, using a phase plate in the WFS design together with reference slopes+pyramid optical gain

We then derived the SCAO performance for various turbulence conditions and point-like reference source magnitudes as well as a global SCAO error budget (Fig. 6). The complete SCAO error budget is obtained from these end-to-end AO simulations and additional instrumental error terms (130 nm rms left for these instrumental contributors to fulfil the SCAO performance specifications of SR=60%, i.e. 250 nm rms).

Since PDR, major developments in terms of AO simulations are concerning:

• The pyramid optical gains estimation, with the development of our new "CLOSE" method that allows a robust estimation, in real-time, of these gains from previous WFS measurements [START_REF] Deo | CLOSE: a self-regulating, best-performance tracker for modal integrator based AO loops[END_REF] (see e.g. https: //tinyurl.com/singleStarCLOSE).

• The performance estimation for double/extended reference sources, [START_REF] Vidal | Analysis the MICADO-MAORY SCAO performance[END_REF] demonstrating basically that the SCAO performance depends on the double star components on sky orientation and separation (Fig. 7 right and https://tinyurl.com/doubleStarCLOSE) when the difference of magnitude between the double star components is larger than 1.75 the performance is equivalent to the single star case the performance is preserved with reference sources with a diameter up to 0.8 (Titan case, https: //tinyurl.com/TitanCLOSE)

• The development of a new island effect management strategy based on a MMSE reconstructor filtering piston per segment together with an estimation of the missing differential piston component based on Kolmogorov statistics 14 (Fig. 7 left)

MICADO-MAORY SCAO DEVELOPMENT PLAN

The CANARDO on-sky demonstration and the SCAO prototyping activities

During the manufacturing phase, in 2021, our development plan includes an on-sky demonstration of our WFS core (double pyramid + PIL + modulation) + RTC prototype using the ALPAO 64×64 actuator DM and an OCAM2 camera (ALICE not yet delivered). It will be performed on the CANARY bench at WHT. The goals of this demonstration are to check the proper functioning of the WFS hardware and to test the algorithms to be implemented in the SRTC, in particular their robustness to variable turbulent conditions (optimized modal gains, NCPA & island effect management, etc).

For that purpose, the manufacturing of the WFS core will be launched in advance to the project manufacturing phase. Tests of the WFS core + RTC prototype will done beforehand in our ATLAS AIT lab.

In parallel, pyramid WFS and K-mirror prototyping is performed in our Sesame lab together with the development of our field selector proto-flight model.

AIT milestones

In Meudon, the SCU (in its AIT configuration) and WFS will be integrated on separate optical benches but tested together. Performed tests will be: alignment and mounting procedures, internal interfaces verification, functional tests, NCPA procedure validation, software (ICS, SRTC) tests and continuous improvement, performance measurements. After the GD structure test with a specific test-bed, the SCU (AIT configuration) and WFS will be integrated on the GD to perform tests similar to the aforementioned ones (in particular AO performance).

In Garching, first tests with MICADO will be in the MPE hall (always with the SCU in its AIT configuration): external interfaces verification (mechanical, cabling, software), common optical alignment, NCPA measurements to build the final phase plates, performance measurements with MICADO. Tests with MICADO in its final configuration at the ESO X4 hall will be done with the SCU in its on-sky configuration. SCAO will be then shipped to Chile.

In Bologna, AIT with MAORY will be done with duplicated or dummy parts to check mechanical interfaces and interlock safety systems.

MICADO MAORY SCAO CALIBRATION STRATEGIES

Calibrations during AIT in Europe

These calibrations will be two-fold:

• Fixed and recurrent general calibrations at sub-system or system level, in particular: matching MICADO vs.

SCAO coordinate systems, phase plate manufacturing, SCU photometry, MICADO plate scale, MICADO slits and coronagraph masks positions in SCAO field selector space, XDM & HODM command matrices, reference slopes for full turn of K-mirror and ADC, PIL interaction matrix, WFS detector calibration • AO loop related calibration: WFS calibration against a modal basis larger than M4 with a large interaction matrix allowing to cope with any DM (provided we know its mode influence functions)

Calibrations during commissioning

These calibrations will either:

• Not use ELT observing time: calibrations done during AIT in Europe will be redone during commissioning for verification (see above)

• Use ELT observing time, requiring a fully functional ELT:

-Verifications of various settings/parameters supposedly known from design (signs/orientations of offsets, of pupil, of K-mirror rotation, etc)

-Verifications related to M4/M5: M4 modes/actuators on-sky orientation, scale of M5 motion, ESOprovided M4 modal influence functions (spatial aspects), RTC latency to M4

-Non-AO calibrations: field selector scale + registration with respect to MICADO detector, secondary guiding adjustment, transmission, etc -AO-loop calibrations: closing the loop with an increasing number of modes, checking for the pupil guiding, M4 registration tracking, pyramid optical gain tracking -TLR verifications: SCAO performance (SR vs. magnitude & seeing)

Calibrations during operations

SCAO is designed with the goal of minimizing the calibration phases during the instrument lifetime at the telescope and then lowering the impact on the observatory operation. In particular, the SCAO calibration unit delivered in the SCAO module is meant to allow us to perform calibrations internally without interfering with any activity at the ELT.

During daytime:

• NCPA can be re-calibrated, on demand according to the on-sky image quality (estimated calibration time: 3h/yr)

• Verification of SCAO good functioning by running each day a ∼10 min template During night time:

• M4 influence functions: we rely entirely on ESO providing these as-built M4 influence functions and no calibration of the influence functions is foreseen

• M4 interaction matrix: this calibration is performed during commissioning at night (no calibration source is available in the telescope) and is not planned to be re-calibrated (relying on ESO confidence in M4 + telescope stability)

• Optimization tasks done in background to observations: modal gains computation for command matrix optimisation and NCPA compensation, M1 & M4 registration, WFS detector background, windshake/vibration management, M1 missing segments

Figure 1 .

 1 Figure 1. MICADO at the telescope in its standalone configuration, i.e. with the SCAO module alone and a dedicated optical relay.

Figure 2 .

 2 Figure 2. Left: The GD with the baffled SCAO parts on the top and the MAORY MCAO NGS WFS at the bottom. Right: Bottom view of the SCAO bench, with the deployable dichoic plate assembly.

Figure 3 .

 3 Figure 3. Left: Opto-mechanical design of the WFS (blue rays) and SCU in its final on-sky configuration (magenta rays). Right: Zoom on the WFS components (the pupil imaging camera is not represented).

Figure 4 .

 4 Figure 4. Left: The SCU in AIT configuration, with a breadboard supporting the AIT components. Right: ALPAO 64×64 DM @ Meudon (co-funded by Région Ile-de-France/DIM ACAV and CNRS/INSU).

Figure 5 .

 5 Figure 5. Left: overview of the different controls and data transfers at play for SCAO, involving both the instrument control software (and the corresponding ESO instrument control system framework modules) and the RTC. Right: HRTC core architecture, featuring the different HRTC SW levels. 10

Figure 6 .

 6 Figure 6. SCAO performance including turbulence, NCPA compensation, windshake and island effect. Left: Performance vs. the guide star magnitude for the different ESO atmospheric conditions. Right: Corresponding error budget for R=12.

Figure 7 .

 7 Figure 7. Left: Performance of various island effect management strategies: actuator slaving (blue), MMSE-based 14 (red), ideal differential piston sensor (green). Black line is the performance in case of no island effect correction. Right: K-band performance for a R=11 double star reference source vs. double star components separation and angle.[START_REF] Vidal | Analysis the MICADO-MAORY SCAO performance[END_REF]

Figure 8 .

 8 Figure 8. Left: MICADO-MAORY SCAO development plan till 2021, emphasizing the prototyping activities and the path to the CANARDO experiment. Right: MICADO-MAORY SCAO development plan after 2021, emphasizing the successive AIT phases.

ACKNOWLEDGMENTS

The purchase of the ALPAO 64×64 actuator deformable mirror has been funded by a grant from DIM ACAV/Région Ile-de-France and thanks to the financial support of INSU. Its tests and characterization have been supported also by INSU and Observatoire de Paris through the funding of J. Raffard's position.


Guillaume Cathelain

Bertrand Rivet

Sophie Achard

Jean Bergounioux

François Jouen

Dynamic Time Warping for Heartbeat Detection in Ballistocardiography

Monitoring vital signs of neonates can be harmful and lead to developmental troubles. Ballistocardiography, a contactless heart rate monitoring method, has the potential to reduce this monitoring pain. But signal processing is uneasy due to noise, inherent physiological variability and artifacts (e.g. respiratory amplitude modulation and body position shifts). We propose a new detection of heartbeats using a dynamic time warping (DTW) template matching method. A heartbeat template is automatically trained and the most similar potential heartbeats in the signal are classified as heartbeats. Data normalizing and DTW steps overcome amplitude and time variability issues. The DTW template matching algorithm has been tested on a 20 to 50 minutes-long BCG records of ten healthy adults in noisy conditions. It is suitable for medical, real-time and low-cost applications.

Introduction

The physiological state (cardiac and respiratory activities) of children in hospital environment is usually evaluated by electrocardiography (ECG) and photoplethysmography (PPG). The former uses electrodes on the thorax and on the limbs extremities; the latter uses pulse oximetry probes that tweak either a finger or a toe. Monitoring children is difficult because of their high mobility that leads to repetitive artifacts and false alarms. Moreover, electrodes peel off the epidermis and pulse oximetry probes can regularly detach from the finger or toe. Several factors are involved: the child twitches and pulls the probes, the electrode adhesive is worn or the nursing staff could not properly set up the electrodes. Children pain, which may be monitored by heart rate variability (HRV) [START_REF] Appelhans | Heart rate variability and pain: Associations of two interrelated homeostatic processes[END_REF] or infrared thermography [START_REF] Brummelte | Procedural pain and brain development in premature newborns[END_REF], affects the neuro-motor and cognitive development especially for preterm neonates [START_REF] Alalwani | Neonatal Pain Recognition using LBP descriptor and Wavelet Thresholding Technique[END_REF]. Finally, this equipment is expensive to install and maintain; given the elements described above, it must be changed several times during hospitalization. Consequently, a nonintrusive apparatus is necessary to detect children heartbeats in a harmless way.

Ballistocardiography is a non-intrusive monitoring method for cardiac activity. It was invented at the end of the 19th century but supplanted by the first electrocardiographs, improving precision and robustness at that time. Today, with new sensor technologies and digital signal processing, this technology is gaining renewed interest. Ballistocardiography's principle relies on measuring ballistic forces [START_REF] Kim | Ballistocardiogram: Mechanism and Potential for Unobtrusive Cardiovascular Health Monitoring[END_REF]: during ventricular systole, blood is ejected from the left ventricle through the aortic arch, generating a pulsed cardiac ballistic force. This mechanical phenomenon allows contactless measurements of a pressure variation (e.g strain gauge) or a deformation (e.g accelerometers) of the patient bed or mattress. A ballistocardiogram (BCG) is a record of this mechanical phenomenon characterized by H to N peaks; however a BCG signal is generally affected by noise such as respiratory, movement artifacts and hardware limitation. Few studies [START_REF] Lee | Ballistocardiogram of baby during sleep[END_REF] focus on pediatric ballistocardiography, where BCG signals are noisier. Compared to adults, the BCG signal amplitude of a 3kg infant is about 30 times lower due to low weight and low cardiac contractile force [START_REF] Erkinjuntti | Use of the SCSB method for monitoring of respiration, body movements and ballistocardiogram in infants[END_REF].

Specific digital signal processing algorithms have been developed for detecting heartbeats, beat-to-beat heart rate and heart rate variability (HRV) in BCG signals using time domain or time-frequency domain methods [START_REF] Inan | Ballistocardiography and seismocardiography: A review of recent advances[END_REF]. In template matching methods, heartbeats' shape is manually [START_REF] Jansen | Monitoring of the Ballistocardiogram and with the Static and Charge Sensitive and Bed[END_REF][START_REF] Shin | Automatic Ballistocardiogram (BCG) Beat Detection Using a Template and Matching Approach[END_REF] or automatically [START_REF] Paalasmaa | Adaptive heartbeat modeling for beat-to-beat heart rate measurement in ballistocardiograms[END_REF] modeled and crosscorrelated to the signal. However the cross-correlation function may be difficult to threshold because of amplitude and time variability of the IJK complex.

In this context, we present a new method to detect heartbeats based on a template matching algorithm improved by DTW and heartbeat normalization. The resulting detection is robust to noise and heartbeats variability. The remaining of this paper is organized as follows: Section 2 presents the materials and Section 3 details the proposed method before results in Section 4. Finally, Section 5 concludes this study.

2.

Materials

A MM0 3000 medical bed has been equipped with a Murata SCA11H BCG sensor. The sensor is placed on top of a MMO Aerospacer overlay in a head-to-foot direction as shown in Figure 1. This wireless sensor is based on a SCA61T accelerometer with a 150 µm/s 2 /√Hz low noise density. The analogue output is AC coupled, anti-aliasing filtered and amplified before digitization at 1 kHz.

The digital signal is filtered with third order Butterworth filters, specifically a 25 Hz low-pass filter and a 2 Hz high-pass filter. They are applied forward and backward in order to prevent phase distortion. Lastly, the signal is decimated to 200 Hz sampling frequency.

Methods

The ultimate goal of the methods is to accurately detect heartbeats in BCG, for heart rate and HRV analyses.

Template initialization

In template matching methods, a template helps recognize similar patterns in the signal; it is often defined manually. In our case, a heartbeat template is automatically initialized for each participant in the following manner. An envelope detection algorithm based on a Hilbert transform is applied on the ten first motionfree seconds of the BCG signal. Local minima, spaced by a minimum delay Δt depending on a hr estimate of the heart rate, are detected. The HRestimate equals 60 times the fundamental frequency of the BCG envelope's Fast Fourier Transform.

Segments starting from one of these local minima to another one are very likely to include one heartbeat each and are thus called inter-beat intervals IBI. Global minima of IBIs are probably J peaks, as illustrated in Figure 2. In Figure 3, IBIs are J peaks synchronized, z-score normalized and averaged altogether. Normalizing IBIs is necessary to preserve the shape of heartbeats and to overcome amplitude variability issues. The initial heartbeat template is parametrized by its length L and the J peak index iJ.

Heartbeat candidates selection

In order to recognize patterns similar to the initialized template, the BCG signal is sliced into L-long IBIs called heartbeat candidates. Local minima of the BCG signals are detected as potential J peaks of heartbeat candidates. Heartbeat candidates are IBIs shortened to L with J peak at index iJ. There are potentially about 25 heartbeat candidates per second because low-pass cut off frequency is 25 Hz; features like minimum distance between peaks or minimum height of peaks can help reduce the number of candidates.

Dynamic Time Warping

Heartbeat candidates must be compared with the template in such a way that time variability issues are overcome. Dynamic Time Warping (DTW) is a method for measuring similarity between two time series which are globally or locally delayed. By aligning the sequences in time, DTW is independent of the non-linear time variability of the sequences. It is widely used for pattern recognition in time series. It has been used in BCG for evaluating the potential changes in the signal induced by different respiratory patterns and posture [START_REF] Martín-Yebra | Evaluation of respiratory-and postural-induced changes on the ballistocardiogram signal by time warping averaging[END_REF].

Optimization, in terms of normalization, endpoints weighting and warping constraint has been performed [START_REF] Mueen | Extracting Optimal Performance from Dynamic Time Warping[END_REF]. Firstly, the initial template and heartbeat candidates are z-score normalized and multiplied by a tapered cosine window in order to smooth endpoints and force a zero value. Secondly, the width of the warping path is constrained to 5 samples by a Sakoe-Chiba band, which corresponds to 25 ms local delays. In fact, as the maximum frequency of the signal is limited to 25 Hz, adjacent local minima are minimally spaced by 40 ms. If time warping were longer than this delay, heartbeat candidates might be confounded.

Heartbeat detection

Heartbeat candidates are sorted in a buffer by their DTW distance with the template. Starting from the smallest DTW distance, heartbeats candidates are discarded from the buffer as soon as they are closer in time than Δt to a heartbeat candidate with smaller DTW distance.

Beat-to-beat delays, heart rate series and HRV are computed from J peak position of the remaining heartbeat candidates.

Template training

The initial template is an average of IBIs synchronized with their probable J peaks for the ten first seconds. Its shape may differ from actual heartbeats, consequently the template must be refined. In the template training step, IBIs are updated with the new detected J peaks and the trained template is the average of IBIs synchronized with these J peaks. The previous steps of optimization, DTW execution and heartbeat detection are then iterated to make sure that no false positive or false negative occurred when the template differed from actual heartbeats. Iterations are executed until the template converges (two iterations are generally enough for a ten-second-long BCG signal). An example is given in Figure 4.

Signal-to-noise ratio estimation

Noisy conditions must be quantified to check the algorithm performance. SNR estimation was estimated for each motion-free BCG signal as the average sample correlation coefficient between pairs of sequential heartbeats [START_REF] Mccall | Standing ballistocardiography measurements in microgravity[END_REF].

Results

The DTW template matching algorithm presented in this paper has been tested on the BCG database described in 2.2. A medical expert manually checked every heartbeat using a dedicated data visualization software. A detected heartbeat is classified as a true positive only if its J peak is at the same location of the J peak of a real heartbeat. Over all records (M = 14139 real heartbeats), false positives (FP), true positives (TP) and false negatives (FN) are reported in Table 1. Specificity was not computed because true negatives are ill-defined. Overall sensitivity and positive predictivity are 95.6 % and 96.8 % respectively, with records' SNR ranging from 2.0 to 14.6. If we separate records whose SNR are lower than 3 (four records) from those higher than 3 (six records), sensitivity and positive predictivity will be 85.3% and 88.5% in the noisiest case; 98.7 % and 99.2 % in the second case.

Discussion

We presented a template matching method for automatically detecting and determining the shape of heartbeats in BCG signals, which solves amplitude and time variability issues by dynamic time warping and normalization. This algorithm was successfully tested on a database with varying conditions of noise. Synchronized gold standard physiological measurements such as ECG were missing, so detected heartbeats were checked manually. This process is tedious especially in noisy conditions where it is difficult to manually label heartbeats. Moreover, the heartbeat detection step of the algorithm forces a minimum delay between two consecutive heartbeats based on an estimation of the average heart rate. This strong assumption might be locally refined using Short-Term Fourier Transform depending on the length of the signal. In addition as heartbeats are selected by their ascending DTW distance with the template, one badly labeled heartbeat often leads to multiple false positives in its neighborhood; this issue could be solved by adding features for the heartbeat detection step.

The adult BCG database, realized in noisy ecological conditions, served as a proof of concept for this algorithm, which will be further tested on infants and neonates BCG databases during experimentations in Raymond Poincaré Hospital -APHP in Garches, France. The presented algorithm enables the possibility to detect heartbeats in BCG signals recorded by usual industrial MEMS sensors that are noisier but far less expensive.

Figure 1 .

 1 Figure 1. Position of the BCG sensor. Ten healthy adult volunteers were monitored for 20 to 50 minutes-long naps during lunch break in Ecole Pratique des Hautes Etudes. This process has been validated by an ethical commission in conformity with the European data legislation. Volunteers were asked to down on the bed and rest in supine position.

Figure 2 .

 2 Figure 2. Flowchart of the methods. Prior to the detection algorithm, motion-free BCG signals were manually isolated. In this paper we focus on these signals. Figure 2 illustrates the methods. Final steps of template training and Signal-to-Noise (SNR) ratio estimation are detailed in Sections 3.5 and 3.6.

Figure 3 .

 3 Figure 3. BCG envelope detection and segmentation.

Figure 4 .

 4 Figure 4. Initialization of the template.

Figure 5 .

 5 Figure 5. Refinement of the template.

Table 1 .

 1 Truth table for detection of heartbeats.

	M = 14139	Heartbeat	No heartbeat
	Detected	TP = 13523	FP = 433
	Undetected	FN = 616

Acknowledgments

The database recording would not have been possible without the support of Ecole Pratique des Hautes Etudes and its volunteers. Medical beds and overlays were lent by MMO, France. The whole study has been financed by Institut Carnot Cognition and PSL University.


Fateme Ghayem

Bertrand Rivet

Christian Jutten

Rodrigo Cabral Farias

GRADIENT-BASED ALGORITHM WITH SPATIAL REGULARIZATION FOR OPTIMAL SENSOR PLACEMENT

Keywords: Sensor placement, Signal extraction, Signal to noise ratio, Alternating optimization, Penalty method

In this paper, we are interested in optimal sensor placement for signal extraction. Recently, a new criterion based on output signal to noise ratio has been proposed for sensor placement. However, to solve the optimization problem, a greedy approach is used over a grid, which is not optimal. To improve this method, we present an optimization approach to locate all the sensors at once. We further add a constraint to the problem that controls the average distances between the sensors. To solve our problem, we use an alternating optimization penalty method. As the associated cost function is non-convex, the proposed algorithm should be carefully initialized. We propose to initialize it with the result of the greedy method. Experimental results show the superiority of the proposed method over the greedy approach.

INTRODUCTION

Whenever a physical phenomenon is under study, it is required to collect some data using specified number of sensors. Several considerations, including reduction of expenses, having ergonomic design, and reducing computational costs, make it necessary to limit the number of sensors. As such, it is essential to place as minimum sensors as possible in those positions that maximize information gathered by the sensors. This is called optimal sensor placement, which has been studied in various domains, with different solutions depending on the application, see, e.g., [START_REF] Berry | Sensor placement in municipal water networks[END_REF][START_REF] Meo | On the optimal sensor placement techniques for a bridge structure[END_REF][START_REF] Hermens | Development of recommendations for SEMG sensors and sensor placement procedures[END_REF][START_REF] Chen | Automatic sensor placement for model-based robot vision[END_REF][START_REF] Daniel | Sensor placement for on-orbit modal identification and correlation of large space structures[END_REF][START_REF] Martinez | Optimal sensor placement and motion coordination for target tracking[END_REF][START_REF] Krause | Efficient sensor placement optimization for securing large water distribution networks[END_REF].

One specific application of the optimal sensor placement problem is to extract a source signal given some measurements obtained from a limited number of sensors. In the classical kriging approaches [START_REF] Shewry | Maximum entropy sampling[END_REF][START_REF] Cressie | The origins of kriging[END_REF][START_REF] Krause | Near-optimal sensor placements in gaussian processes: Theory, efficient algorithms and empirical studies[END_REF], this problem is studied in two steps. In the first step, the optimal sensor placement problem is solved such that the sensor recordings provide a good estimation of the spatial map. Afterwards, in the second step, an extraction approach such as a linear extraction model is used to extract the source signal from the sensor measurements.

Recently, by assuming a stochastic model, a robust sensor placement criterion based on the maximization of the aver-age signal to noise ratio (SNR) has been proposed [START_REF] Ghayem | Optimal sensor placement for signal extraction[END_REF]. This method works by maximizing the average SNR of the desired source, and has shown a significant improvement compared to the classical kriging approaches in terms of the output SNR. To solve the maximization problem, a greedy approach is used in which the sensors are added one by one from a fixed grid of candidate sensor positions. Two limitations can be observed: first the sensor location are restricted to be on a predefined grid and secondly the greedy approach is sub-optimal. Consequently, to be precise, the grid should be fine enough leading to a high computation cost. To solve these issues, we propose a first order optimization-based approach that in contrast to the one-by-one strategy adopted by the greedy method on a grid, optimizes all the sensor positions at once and does not discretize the search space. Finally, since placing 2 sensors very close to each other may not be feasible, e.g., due to the physical size of the sensors, a regularizing term is added to avoid to choose too close sensor positions. As the cost function is non-convex, to avoid sub optimal local minima, the proposed algorithm is initialized by the solution of the greedy approach.

The rest of this paper is organized as follows: in Section 2 a review on the problem of optimal sensor placement for source extraction, as well as the criterion suggested in [START_REF] Ghayem | Optimal sensor placement for signal extraction[END_REF] are presented. The details of the proposed method are discussed in Section 3. Section 4 presents the performance of the proposed method and compares it with greedy approach, before the conclusion in Section 5.

BACKGROUND

In this section we present an overview on signal extraction for optimal sensor placement suggested in [START_REF] Ghayem | Optimal sensor placement for signal extraction[END_REF]. Let y(x, t) denote the observation recorded at location x and at time t which is modelled as

y(x, t) = a(x)s(t) + n(x, t), (1)
where a(x) is the spatial gain between the signal of interest s(t) and the sensor at location x ∈ R D , and n(x, t) is a spatially correlated additive noise. The noise is assumed to be uncorrelated with the source. Assuming that M sensors are located at positions X M = {x i } i∈{1,•••,M } , linear source extraction amounts to design a vector f ∈ R M to estimate the source as ŝ(t) = f T y(X M , t), where y(X M , t) = [y(x 1 , t), . . . , y(x M , t)] T . To find the best f , a classical criterion is the output signal-to-noise ratio (SNR) defined by [START_REF] Ghayem | Optimal sensor placement for signal extraction[END_REF] SN R(f

) = E[(f T a M s(t)) 2 ] / E[(f T n M (t)) 2 ], (2)
where a M = [a(x 1 ), . . . , a(x M )] T denotes the vector of spatial gains and n M (t) = [n(x 1 , t), . . . , n(x M , t)] T the vector of noise. The spatial gain a(x) of the source of interest is modelled as a stochastic Gaussian process to take into account the uncertain knowledge on it:

â(x) ∼ GP(m a (x), k a (x, x )), (3)
where m a (x) is the mean function and k a (x, x ) is the covariance function. Furthermore, we consider a zero mean Gaussian process with covariance matrix k n (x, x ) to model the noise. The criterion to optimize is then obtained [START_REF] Ghayem | Optimal sensor placement for signal extraction[END_REF]:

J(X M ) = (m a M ) T (R n M ) -1 m a M + Tr((R n M ) -1 R a M ), (4)
where R a M ∈ R M ×M and R n M ∈ R M ×M are the covariance matrices of the spatial gain and the noise respectively, whose (i, j)th elements are k a (x i , x j ) and k n (x i , x j ). Also,

m M = {m a i } i∈{1,•••,M } is the set of means at locations {x i } i∈{1,•••,M } ,
and Tr(•) is the trace operator. The optimal sensor locations are finally obtained as

XM = arg max X M J(X M ). (5)
Directly maximizing (4) in a grid requires a combinatorial search, which leads to a high computational cost. Therefore, a greedy approach has been introduced in [START_REF] Ghayem | Optimal sensor placement for signal extraction[END_REF] that selects the M sensors by sequentially selecting N < M sensors at a time. Assuming that K sensors have already been placed, to choose the locations of the next N sensors, the following criterion is optimized:

J(X N |X K ) = E[â T K+N (R n K+N ) -1 âK+N | X K ], (6)
where

K + N means {X N X K } and thus âK+N ∈ R K+N can be divided as âK+N = [â T K , âT N ] T .
Once the sensor locations XM are obtained, the source of interest is extracted by using the following separation vector [START_REF] Ghayem | Optimal sensor placement for signal extraction[END_REF]

fM = (R n M ) -1 m a M . (7)

PROPOSED METHOD

In this section we present our proposed framework to solve the optimization problem for sensor placement. Unlike the greedy approach, our proposed method directly provides the positions of all the required number of sensors. By considering a one dimensional situation, we want to minimize f (x M ) = -J(x M ). In order to control the average distances between each pair of the sensors, we constrain the sum of the squared distances to be greater than a threshold. Furthermore, due to the spatial constraints of the boundaries, we consider a normalized case where 0 ≤ x i ≤ 1. Therefore, we study the following minimization problem:

min x M f (x M ) s.t. Dx 2 2 ≥ 0 ≤ x i ≤ 1 i ∈ {1, 2, . . . , M } , (8)
where

D ∈ R M (M -1) 2
×M is a matrix that enumerates all the possible combinations of positions in pairs of size two. For instance, if the number of sensors is M = 3, then

Dx =   1 -1 0 1 0 -1 0 1 -1     x 1 x 2 x 3   =   x 1 -x 2 x 1 -x 3 x 2 -x 3   .
To solve (8), we define an auxiliary variable z M = Dx M , and reformulate (8) as the following problem:

min x M ,z M f (x M ) s.t.      z M ∈ A , z M = Dx M , 0 ≤ x i ≤ 1 i ∈ {1, . . . , M }, (9)
where

A = z M ∈ R M z M 2 2 ≥
. To solve (9), we use the penalty method [START_REF] Bertsekas | Nonlinear Programming[END_REF], by adding the constraint z M = Dx M as a penalty to the target function with the penalty parameter α:

min x M ,z M ∈A f (x M ) + 1 2α z M -Dx M 2 2 s.t. 0 ≤ x i ≤ 1 i ∈ {1, . . . , M }. ( 10
)
To solve [START_REF] Krause | Near-optimal sensor placements in gaussian processes: Theory, efficient algorithms and empirical studies[END_REF], we do the minimization over x M and z M alternately. At iteration l, first, the cost is optimized over z M , fixing x M to its current estimate x

M . That is:

z (l) M = argmin z M ∈A 1 2α z M -Dx (l) M 2 2 . (11)
The solution to the above minimization is a projection onto the set A as follows:

z (l) M =    Dx (l) M , if Dx (l) M 2 2 ≥ Dx (l) M Dx (l) M 2 2
, otherwise.

For the second step, the variable z M is fixed as in [START_REF] Bertsekas | Nonlinear Programming[END_REF], and we do the minimization over x M as follows:

x

(l+1) M = argmin x M f (x M ) + 1 2α z (l) M -Dx M 2 2 s.t. 0 ≤ x i ≤ 1, i ∈ {1, . . . , M }. ( 13
)
Since the constraint is a quite simple convex set, to solve (13), a projected gradient descent is used: after a gradient descent (GD) update, the result is projected onto [0, 1]. That is, by defining the cost function as follows:

g(x M ) = f (x M ) + 1 2α z (l) M -Dx M 2 2 , (14)
the gradient step to optimize (13) is

x (l+1) M ∂xi = [∂m a j /∂x i ] j , and ∂C ∂xi = [∂C ij /∂x i ] (i,j) in which C represents any covariance matrix, and C ij corresponds to its (i, j)th entry. This way, we have computed the gradient of f (x (l) M ) over the i th sensor position, providing thus the expression of the gradient vector

∇ x M f (x (l) M ).
To determine µ in (15), we use a backtracking line search strategy [START_REF] Bertsekas | Nonlinear Programming[END_REF]. After updating x M using (15), any element of x

(l+1) M is projected into [0, 1].
Finally, to solve [START_REF] Krause | Near-optimal sensor placements in gaussian processes: Theory, efficient algorithms and empirical studies[END_REF] we start with an initial point and alternate between the projection step [START_REF] Bertsekas | Nonlinear Programming[END_REF] and the GD step (15). As done in standard penalty methods [START_REF] Bertsekas | Nonlinear Programming[END_REF], the problem (10) should be solved for a decreasing sequence of α, e.g., as {α 0 , α 1 , . . .} where α k+1 = ηα j , with 0 < η < 1. For each fixed value of α, we perform a few iterations between (12), (15) projecting x. Moreover, iterations corresponding to α j+1 are initialized by the final estimate found for α j . The final algorithm to solve (8) is summarized in Algorithm 1.

Since the problem ( 9) is non-convex, its initialization is important to find an appropriate minimizer. We propose to initialize the algorithm with the solution obtained by the greedy approach [START_REF] Ghayem | Optimal sensor placement for signal extraction[END_REF]. In this way, the algorithm is more likely to end up with a good local minimum.

Algorithm 1 Alternating minimization (AM) for solving [START_REF] Cressie | The origins of kriging[END_REF] 1: Inputs:

x (0) , z (0) , µ 0 , α 0 , Q 2: Initialization: Set µ = µ 0 , l = 0 3: for j = 1, 2, • • • , Q do 4:
while stopping criterion not met do 5:

t (l) = ∇g(x (l) ) 2 2 6: while g(x (l) -µ∇g(x (l) )) > g(x (l) ) -µ 2 t (l) do 7: µ ← β • µ 8:
end while 9:

x (l+1) = x (l) -µ∇g(x (l) )

10:

Project

x (l+1) into [0, 1] 11:
Perform projection ( 12) to obtain z (l+1)

12:

l ← l + 1

13:

end while 14:

α j+1 = η • α j 15: end for 16: Output: x (l)

NUMERICAL EXPERIMENTS

In this section, the numerical setup is first presented. Then the influence of the initialization, the effect of the regularization based on the sensors distances and of the smoothness of the spatial gain are presented in Subsections 4.2, 4.3 and 4.4, respectively.

Numerical setup

Synthetic data are generated in a 1D space, where the range of the sensor locations is normalized between 0 and 1. We consider a prior on the spatial gain and noise to be generated from GP(m(x), C(x, x )), with a square exponential covariance function C(x, x ) = σ 2 exp(-(x -x ) 2 /(2ρ 2 )). The mean of the noise is set to be 0. The mean of the gain is given by m a (x) = i=1 γ i sin di (w i πx), where, γ i , d i and w i are the i th elements of the vectors G = [0.1, 0.3, 0.5, 0.7, 0.9], D = [1, 1, 3, 1, 2], and W = [25, 26, • • • , 29], respectively. The smoothness parameters ρ n , and ρ a , and the variances σ n and σ a as well as the size of the spatial grid for greedy initialization take different values for each experiment. Also, we set α 0 = 1, Q = 50, η = 0.5, µ 0 = 1, and β = 0.5.

Influence of the initialization

We set the size of the spatial grid to be 100. Two different values of the uncertainty on the spatial gain are considered, σ a = 1 and σ a = 3. The noise variance σ n is accordingly set such that the SNR becomes 0.8 dB. The smoothness of the uncertainty on the spatial gain (ρ a ) is set to ρ a = 0.001, which corresponds to an uncertainty with almost no spatial correlation. The spatial smoothness of the noise ρ n is set to ρ n = 0.01ρ a . The lower bound on Dx 2 2 is also considered to be = M (M -1) 2 × 10 -3 . The true value of the output SNR (2) computed with estimated extraction vector fM (7) versus the number of desired sensors M are depicted in Fig. 1. If no prior information is used for the initialization, one can initialize the sensor locations regularly-spaced between 0 and 1. On the contrary, one can use the previously proposed greedy approach, where each sensor are added one by one [START_REF] Ghayem | Optimal sensor placement for signal extraction[END_REF], to initialize the sensor locations. Firstly, one can see that the greedy initialization leads to a better extraction of the source s(t) than using regularlyspaced initial locations for the sensors before applying our proposed method to adjust the sensor locations. Indeed, the difference of output SNRs varies between 10dB for a single and 5dB for 25 sensors. Moreover, the proposed method to adjust the sensor locations leads to improve the SNR of about 3dB to 5dB compared to the greedy approach. Indeed, this result is expected since the proposed method tackles the optimization of the sensor locations all at the same time instead of one after the other as in the greedy method. It is also worth noting that the output SNR is worse by applying the proposed method with a regularly-spaced initialization than by just choosing the sensor locations by the greedy method proposed in [START_REF] Ghayem | Optimal sensor placement for signal extraction[END_REF] with no additional adjustment.

Regularizing sensors distances

Figure 2 shows the effect of regularizing sensor distances and how it can help to control the average distances between pairs of sensors. In this part, all the parameters are set as in the previous section with σ a = 1, except that here we consider a tighter grid of size 320. Also, the number of desired sensors are set to be M = 15. For the proposed method, two different values for the lower bound are considered: = 0.5 and = 1.

The second and the third sub-figures demonstrate the effect of in tuning the average distances between the sensors. The final SNR values for each approach, from top to bottom, are 29.22 dB (initial SNR), 32.22 dB ( = .5), and 31.53 dB ( = 1), respectively. However, increasing leads to only a slight decrease of the output SNR while increasing the average distance between the sensors.

Effect of the smoothness parameter ρ a

In this part we study the performance for different smoothness levels of the uncertainty on the spatial gain (ρ a ). We consider an almost difficult situation with the uncertainty parameter of the spatial gain to be σ a = 5. We also considered SNR to be Fig. 3. Effect of the smoothness parameter of spatial gain ρ a .

2 dB to set σ n . The rest of the parameters are set similarly as in subsection 4.2. In Fig. 3, the output SNR versus the degree of the spatial gain smoothness is depicted. We can see that as the signal becomes more non-smooth, the performance of the greedy approach deteriorates much faster than the proposed method. This is due to the presence of highly informative sensor positions in between grid points, which cannot be chosen by the greedy approach.

CONCLUSIONS AND PERSPECTIVES

In this paper, we addressed optimal sensor placement for signal extraction by maximizing output signal to noise ratio. In contrast to the greedy approach proposed in [START_REF] Ghayem | Optimal sensor placement for signal extraction[END_REF], the proposed method adjust all the sensors locations at once instead of choosing them one at a time. To this end, a gradient-based method is proposed to search for the sensor locations over the whole space. A constraint, controlling the average distances between sensors, is also considered to avoid to choose too close sensors (e.g., depending of their size). Due to the non-convexity of the cost function, the proposed algorithm is initialized with the solution of the greedy approach. Experimental results demonstrate that the proposed method provides about 3 dB improvements over the greedy approach. Also, thanks to the new constraint, the proposed method is shown to be able to control the average distances between the sensors. In future works, an explicit constraint on each distance between pair of sensors will be studied as well as other global optimization algorithms to avoid convergence to a local optimum.

Fig. 1 .

 1 Fig. 1. Influence of the initialization. Output SNR vs. the number of sensors. left: σ a = 1, right: σ a = 3.