Dr V Malgras

Dr J Henzie

Dr T Takei

Prof Y Yamauchi
email: y.yamauchi@uq.edu.au

Keywords: mesoporous silica, mesoporous films, perovskite nanocrystals, confinement, nanospace

Stable Blue Luminescent CsPbBr 3 Perovskite Nanocrystals Confined in

Mesoporous Thin Films

Victor Malgras [a] , Joel Henzie [b] *, Toshiaki Takei [b] , and Yusuke Yamauchi [c,d,e] * Abstract: Creating CsPbBr 3 perovskite nanocrystals with bright blue emission is challenging because their optical properties depend sensitively on structure, which can be influenced by colloidal purification procedures. Growing perovskites in mesoporous templates bypasses some of these purification issues because the size of the nanocrystal is governed by the dimensions of the pores. Mesoporous silica consisting of aligned channels with tunable diameter can be easily synthesized and used as template. In SBA-15 films the mesoporous channels are interconnected by numerous ~1.5 to 4.0 nm diameter (d) smaller pores. These openings help fill the templates with perovskite precursor solution because they prevent dead-end channels. When the perovskite solution evaporates and retreats, some of the liquid remains trapped in the interconnecting pores by discontinuous dewetting. The precursor crystallizes, generating stable d ~ 3.1 nm blue-emitting perovskite nanocrystals. The mesoporous template also serves as a protective barrier to preserve the optical properties of the CsPbBr 3 from atmospheric conditions. Besides the significant blue-shift of the luminescence peak, we observed a decrease of the longer lifetime component in the film, accompanied by a strong increase of the external quantum efficiency.

The optoelectronic properties of bulk lead halide perovskites have enabled a whole new class of inexpensive, efficient, solution processable light-emitting and light-harvesting devices. Doping and alloying approaches have long served as an effective means of tuning the optical properties of semiconductors. But nanostructuring is a parallel and complementary strategy to control the excitonic band gap, as well as relaxation, recombination and energy transfer dynamics. [1] Colloidal lead halide perovskite nanocrystals (NCs) have been generated using various routes including hot injection synthesis, [2] room temperature crystallization, [3] and various microemulsion methods. [4] Methylammonium lead bromide (MAPbBr 3 ) perovskites NCs are known for their bright green emission, possessing photoluminescence quantum yields (PLQY) close to unity. [5] Cesium lead bromide (CsPbBr 3 ) NCs are an all-inorganic perovskite counterpart that are more stable than MAPbBr 3 , and have already achieved PLQY values > 90%. [2] To our knowledge, there are few reports of blue-emitting cesium-based perovskite NCs because band gaps in the blue region require the nanoparticle to have one dimension < 4 nm. This is challenging due to material stability issues, and to the polydisperse nature of most synthetic colloidal methods. [6] Ultrathin two-dimensional (2D) CsPbBr 3 nanoplates exhibit blue emission, although their ionic nature makes them susceptible to degradation. [7] Integrating these kinds of colloids into devices is a challenging task, because their sensitivity to water will require all the synthesis, purification and fabrication to be carried out in an inert atmosphere.

To realize a simple and reliable bench-top method to generate stable blue-emitting all-inorganic perovskite nanoparticles, we used ordered mesoporous silica films as templates to direct and confine the growth of nanocrystals. Nanoporous matrices offer numerous advantages in the synthesis of nanostructured lead halide perovskite materials in particular.

[8]

Mesoporous templates provide a ligand-free approach to prevent aggregation, and depending on the dielectric properties of the template it can interconnect or isolate the NCs electronically and optically. For example, silica is ideal for studying the properties of isolated nanoparticles, because its insulating nature contributes minimal collateral electronic interactions. In addition, the natural encapsulation limits adventitious exposure to water, promoting NC stability. Compared to their powdered counterparts, mesoporous thin films yield larger periodic domains (> 1 µm) with a fixed orientation relative to the substrate. Well-ordered structures can even lead to monodisperse particles arranged according to the negative replica of the mold. We envisioned film templates as a platform towards functional composites that are "ready-for-integration", whereas mesoporous powders require further processing steps that cause additional challenges. Finally, films produce negligible scattering in the active region of the device since, unlike the typical micrometer grain size in powders, their features are far below the diffraction limit of light (<< λ).

Capillarity is the main driving force for infusing the channels of the template with perovskite precursor solution. As the diameter (d) of the channel decreases below some threshold, nanofluidic effects begin impacting viscosity and restraining mass transport, thus limiting effective filling of the template. [9] In order to circumvent this nanofluidic hindrance, we sought to exploit a well-known property the mesoporous SBA-15 family: its hierarchical porous composition. The main mesochannels are the negative replica of the cylindrical micelles arranged in a 2D hexagonal phase and separated by the aqueous domain where the hydrolitic condensation of SiO 2 occurs. Two additional categories of micropores are also present: [10] i) ultramicropores (< 1 nm): small pockets that are the imprint of the hydrophobic polyethylene oxide (PEO) subunits on the surface of the mesochannels; ii) larger pores (1.5 -4 nm): likely form when adjacent micelles come in contact or merge during the condensation of SiO 2 . The perovskite precursor solution can access the entirety of the mesoporous film because the main mesoporous channels have less resistance while the interconnecting pores prevent dead-end channels (see Scheme 1). When the precursor solution begins to evaporate, it retreats and dewets the surface of the channels. Surprisingly, this solution dewetting effect spontaneously forms 3.1 nm CsPbBr 3 perovskite NCs in the interchannel pores. Some larger NCs formed in the mesochannels, but they are a minority product. This approach enables the fabrication of blue-emitting nanocomposites consisting of stable all-inorganic lead halide perovskite NCs in a pre-synthesized mesoporous silica film.

Scheme 1. Schematic illustration of the hierarchical filling (left image) where the precursor solution infuses the network via the larger channels and smaller interconnecting pores. As the solution evaporates and dewets the channels, small amounts are trapped in the smaller pores, crystallizing to form CsPbBr 3 NCs (right image). Although some crystallization occurs in the larger main channels to form green emitting NCs, it mostly takes place in the interconnecting pores, thus the primary emission comes from the small blue-emitting NCs.

The mesoporous films are prepared by spin-coating glass or silicon substrates with a sol-gel solution of silica precursor (tetraethylorthosilicate) containing a precise concentration of triblock copolymer (EO 20 -PO 70 -EO 20 , Pluronic P123) that yields a 2D hexagonal mesophase of in-plane cylindrical micelles. The polymer and residiual solvent was then removed via thermal treatment (see supporting information for detailed procedure). An SEM image of the top surface of the pristine mesoporous silica film (Figure 1a) highlights the curved mesochannels domains that are uniformly distributed throughout the sample and appear continuous over several micrometers. After casting and annealing the CsPbBr 3 perovskite solution onto the mesoporous substrate, the surface of the sample appears to be filled with nanoparticles, evidenced by the bright contrasting regions inside the channels (Figure 1b). The average channel width and the 2D hexagonal mesoporous periodicity remains unaffected, as supported by the SAXS patterns in the inset (pore-to-pore distance of nm). The semi-elliptical shape of the SAXS pattern is the direct consequence of the distortion induced on the mesoporous structure. It is known that after calcination of the templated film, the periodicity along the vertical axis shrinks and the pores adopt an elliptical shape. [11] The STEM image in Figure 2a shows a typical region of the mesoporous film containing the perovskite nanoparticles. Large crystals arranged along the direction of the channels can be observed, suggesting an effective penetration of the precursor solution during the synthesis. The lateral length of these particles (perpendicular to channel orientation) spans from 5.0 to 8.0 nm. This broad size distribution is caused by the elliptical shape of the mesopores, observed from random orientations under TEM. The complementary high-resolution (HRTEM) images (Figure 2b-c) evince the crystallinity of the perovskite particles, with lattice fringes aligned in the direction of the encapsulating silica channels. FFT analysis of the image (Figure 2d) indicates the d-spacing is 2.94 Å, which matches the (004) interplanar distance of cubic CsPbBr 3 perovskite (Pbnm, ICSD card #97852). This observation suggests that the crystals grow in the <001> direction along the channels. The XRD pattern of the composite film reveals several weaker peaks that support this assignment (see Figure S1). There is an ongoing debate right now as to whether such perovskite NCs adopt the cubic or orthorhombic structure, [12] and the slight discrepancies in their respective lattice parameters cannot be resolved with the available instrumentation. In addition to these NCs, we noticed that the template is riddled with with much smaller particles in the 2.5 to 3.6 nm size range. This bimodal particle size distribution is highlighted by the histogram in the inset of Figure 2a. The smaller particles were not observed in the powder sample (see Figure S2), although the same precursor concentration was used. This indicates that the smaller particles may originate from how the precursor solution infiltrates and dewets the interior spaces of the template. Furthermore, the adsorption/desorption isotherms revealed that although the specific surface area of the powder sample is considerably larger, the overall contribution of the interconnecting pores in the film sample is more important (Figure S3). As shown in Figure 1, the curved mesochannel domains span the surface of the film. Thus infiltration and transport of precursor solution can only pass through the top during the spincoating process. It is possible that the centrifugal force generated by the spin-coating process helps fill the channels. But as mentioned previously, the mesochannels in SBA-15 are interconnected by micropores. [13] A network of large channels and small pores offers innumerable paths for the perovskite precursor solution to permeate the entirety of the film template. When the solution evaporates and retreats, it dewets the surface of the template. The interchannel pores have a lower surface free energy than the solution, so they likely remain filled via a discontinuous dewetting mechanism that is well explained in 2D systems and has been used to grow colloidal nanocrystals.

[14] Therefore, the interchannel pores in the films offer two advantages: (i) they offer more paths for the precursor solution to fill the film and (ii) they act as a local discontinuity to trap the solution and direct the growth of perovskites NCs that are small enough to experience quantum confinement and enable blue emission. Figure 3a shows a photograph of CsPbBr 3 perovskite NCs in mesoporous powder and mesporous film. The blue emission from the film was further confirmed by the clear displacement of the photoluminescence (PL) peak (Figure 3b) compared to the mesoporous powder and bulk perovskite samples. The bulk perovskite has a peak at 525 nm, which is shifted to 507 nm in the mesoporous powder because the quasi-spherical particles experience moderate 3D quantum confinement, expanding the band gap from 2.36 to 2.45 eV. The PL peak from the mesoporous film further shifts to 487 nm, which indicates the particles are much smaller and experience additional confinement (the reproducibility over nine different samples can be assessed in Figure S4). Although a band gap of 2.55 eV for crystals ~3.1 nm in size finds supporting evidence in past literature, [8b,15] other publications also reported different values. [2,8d,16] Compared to the bulk, the emission line width of the film sample remain narrow (96 and 108 meV, respectively), suggesting good crystal quality, while the peak for the larger particles in the powder is much broader (144 meV). Due to the location of the emission feature, we can conclude that most of the PL signal emitted from the film sample originates from the smaller crystals formed in the interchannel pores observed by TEM. The peak is asymmetrically shaped, and a double-gaussian fit (white curves) reveals the presence of a weak contribution located at 510 nm, which can likely correspond to the larger 7 nm particles. This is in good agreement with the bimodal distribution discussed earlier. In addition, the resulting excitonic band gap supports the TEM and XRD data that these are likely similar in phase as the larger perovskite NCs. It is also possible that variations in the surface stoichiometry can be partly responsible for this shift, as Br-rich compositions tend to blue-shift the PL features. The overall mechanism describing the formation of nanocrystals in confined spaces is not fully understood, but we assume that the inner walls of the interconnecting pores could impart considerable residual pressure on the surface of the NCs to quench their growth. According to previous study on methylammonium lead halide, it was evidenced that external pressure tends to induce a blue-shift in the PL emission due to amorphization. [17] We decide to rule out this effect, however, since the sharp and intense emission peak does not indicate a disordered structure.

The time-resolved PL curves in Figure 3c were fit numerically for each sample with a bi-exponential decay (see Table S1 in the supporting information). Biexponential fits are frequently used to describe the dynamics of this family of NCs. [12b,18] In our case, the mixed perovskite/air and perovskite/SiO 2 interfaces have dissimilar optical-mode densities and refractive indices and are therefore likely to provide different transition dynamics. [19] The average luminescent lifetime is assessed to be 24.0 ns, thus much shorter than that of the bulk perovskite and mesoporous powders (equal to 41.5 and 53.7 ns, respectively). Confining the growth of CsPbBr 3 in a such small volume also led to a substantial increase of the PLQY up to 31.8 % for the films versus 1.6% and 4.1% for the bulk and mesoporos powder sample, respectively. All optical parameters for the three samples are listed in Table S1. The fast component of the lifetime is slightly larger in the film than in the powder. This result seems counterintuitive at first because the radiative recombination rate typically increases inversely with the crystal size due to larger oscillator strength. The excitonic spontaneous recombination lifetime of surfactant-stabilized CsPbBr 3 perovskite NCs, however, has been reported with values ranging between 2 and 8 ns, [12b,18b] suggesting a strong dependence on factors such as the synthetic conditions, morphology and surface passivation.

On the other hand, the mesoporous film has a faster (slow component) and larger PLQY than the bulk and powder composite, which results in a strong reduction in the calculated non-radiative lifetime. Presently there is still vigourous debates about the photophysical events that are responsible for the tail of the PL time-response in this family of perovskite materials. If can be attributed to trap-assisted mechanisms involving charge capture, thermal release and delayed radiative/non-radiative decays, [20] then our results indicate the nanocrystals generated in the interconnecting pores are well-passivated single crystals that contain minimal defects, hence only a limited amount of surface recombination centers. Others have suggested that photon recycling might lead to longer apparent recombination lifetimes, accompanied by a drastic drop of the external yield compared to the internal quantum efficiency. [21] Here, our results support also this mechanism because the 3 nm crystals in the film are scattered and monodisperse, which reduces the probability for photons to be recaptured. S1.

The PL emission of the CsPbBr 3 mesoporous films should be very stable to various environmental conditions because the NCs are embedded in a composite film. To test stability, one of the film samples was kept under atmospheric conditions (relative humidity 60-70%) in the dark for five months, while another was placed under solar simulated light (λ > 300 nm) for 6 h (see Figure S5). In both samples, the PL emission spectra and PL lifetimes and PLQY remained essentially unchanged, showcasing its strong resistance to degradation due to adventitious exposure.

In summary, we demonstrated a simple route to obtain stable blue-emitting CsPbBr 3 NCs confined in the interconnecting pores present in adjoining walls of a mesoporous SiO 2 film. These pores are very small and can trap the equivalent of hundreds of yotoliters ( 10-24 liters) of precursor solution that deposits ultrasmall blue-emitting CsPbBr 3 NCs. We provided sufficient evidence, through electron microscopy and optical measurements, to demonstrate that the blue shift of the film compared to the bulk ( ~ 0.2 eV) is mainly due to strong quantum confinement, although some stress-induced crystallographic distortions are somewhat likely to play a minor role. Compared to the bulk and mesoporous powder sample, the film possesses a shorter PL lifetime and a much higher PLQY. This decrease in non-radiative recombination is combined with a strong resistance to photochemical and moisture-based degradation.

Experimental Section

The starting precursor solution consisted of a stoichiometric mixture of CsBr and PbBr 2 dissolved in DMF (0.03 M) stirred at room temperature overnight. From this solution, a bulk reference sample, a powder composite and a film composite were all prepared. The bulk sample was obtained by drop casting the precursor solution on a glass slide. The powder sample was prepared via the previously reported method: [8a] in brief, the perovskite solution was added dropwise to SBA-15 silica powder (synthesized according to [22]) before being vortex mixed vigorously (1 mL/mg). For the composite film, first a sol-gel solution consisting of a mixture of TEOS, P123, H 2 O, HCl and ethanol (molar ratio 18.5:1.5:100:0.02:110) was prepared and aged for 3 h, according to a previously published method. [23] The solution was then spin-coated on various substrates (glass, silicon) at 3,000 rpm for 15 s, aged overnight and heat-treated at 400 °C for 5 h, with a ramping of 1 °C/min. The perovskite solution is then spin-casted at 3,000 rpm for 15 s on the mesoporous film. Finally, all three samples (bulk, powder and film) are heat treated on a hotplate at 80 °C, under mild vacuum, for 30 min, followed by 150 °C, in air, for 30 min.

Table of Content

Creating CsPbBr 3 perovskite nanocrystals with bright blue emission is challenging because colloidal isolation and purification procedures cause instability. Growing perovskites in mesoporous templates bypasses purification issues because the size of the nanocrystal is governed by the diameter of the pores. Films of aligned mesoporous silica templates can be assembled on surfaces with pore diameters of ~7.5 nm. These mesopores overlap in numerous locations generating tiny, ~1.5 to 3.0-nm diameter interconnecting pores. The hierarchical nature of the porous structure enabled us to easily fill these interconnecting pores with perovskite precursor, and generate stable blue-emitting perovskite nanocrystals (~3.1 nm) that are well-distributed throughout the film. The mesoporous template also serves as a protective barrier to preserve their optical properties from atmospheric conditions. We observe a decrease of the longer lifetime component in the film, accompanied by a drastic increase of the luminescence quantum yield, resulting in a reduction of nonradiative recombination.

Figure 1 .

 1 Figure 1. SEM images of a) the pristine mesoporous silica mesoporous film with channels periodically aligned over > 1 µm, parallel to the surface; b) the CsPbBr 3 perovskite-infused mesoporous silica film with several particles (brighter contrast) enclosed inside the channels. Inset: SAXS pattern of the nanocomposite film sample shown in b), attesting that the overall hexagonal periodicity is maintained despite the perovskite inclusion.

Figure 2 .

 2 Figure 2. a) STEM image of the CsPbBr 3 -infused silica mesoporous film nanocomposite. Inset: histogram highlighting the bimodal particle size distribution, counted over 400 particles. b) Low-and c) high-magnification HRTEM image emphasizing on the CsPbBr 3 NCs encapsulated in the mesoporous silica film. d) FFT of the NCs shown in c).

Figure 3 .

 3 Figure 3. a) Photograph comparing the powder and film samples under UV light. b) PL spectra comparing the position of the main emission peak for the three samples. The solid white line represents a double-gaussian fitting of the emission peak from the film sample. c) Time-resolved PL decay (λ ex = 365 nm) of the three samples, along with their respective bi-exponential fitting. Fitting parameters are reported in TableS1.

Acknowledgements

This study was supported by an Australian Research Council (ARC) Future Fellow (FT150100479), the Suzuken Memorial Foundation, and JSPS KAKENHI (Grant Numbers 17H05393 and 17K19044). The authors would like to thank Dr. Shirai (NIMS) for helpful suggestion and discussion on materials fabrication.


Antônio Padilha

Lanari Bó

Philippe Poignet

Christian Geny

Filtering Voluntary Motion for Pathological Tremor Compensation

This paper describes an algorithm to estimate tremor and voluntary motion from measured motion data. Estimation is performed by means of an extended Kalman filter (EKF), which also estimates tremor parameters. Comparison of the proposed method with techniques described in the literature are conducted with two experimental data sets from tremor patients performing the same task, drawing a spiral. The presented algorithm may be directly applied in real-time pathological tremor compensation systems.

I. Introduction

Tremor may be defined as an involuntary, approximately rhythmic and roughly sinusoidal movement [1]. Although pathological tremor itself is not lifethreatening, it may decrease considerably the person's quality of life, since the ability to perform simple daily tasks is reduced. It is the most common movement disorder found in human pathology.

An absolutely effective treatment for pathological tremor is not yet available. Although pharmacological and surgical therapies are available, they still present limitations with reference to cost, risks, and effectiveness. In this scenario, one alternative is the use of assistive devices to reduce the effects of the abnormal motion. Among the new technologies, upper limb robotic exoskeletons [START_REF] Rocon | Design and validation of a rehabilitation robotic exoskeleton for tremor assessment and suppression[END_REF] and the use of FES [START_REF] Prochazka | Attenuation of pathological tremors by functional electrical stimulation I: Method[END_REF] have been proposed.

Within the context of assistive devices to attenuate tremor, every candidate technology must continuously estimate tremor during the device operation. This is due to the fact that pathological tremor presents high variability throughout time, particularly in terms of tremor amplitude. Our effort to address this problem has been described in [START_REF] Bó | Online pathological tremor characterization with Kalman filtering[END_REF], where an algorithm to perform online tremor characterization based on two distinct parametric models was presented.

However, the sensors used to measure the current tremor state are also sensitive to the voluntary motion performed by the patient, which is normally executed in lower frequencies, when compared to pathological tremor frequencies. Hence, in order to correctly estimate the tremor motion, the voluntary component must be extracted from the total movement. Some methods have already been proposed to address this problem. [START_REF] Gonzalez | Optimal digital filtering for tremor suppression[END_REF] presented the application of digital filters to suppress tremor from computer input, but the method is not feasible for a real-time compensation system. Regarding methods applied to real-time compensation problems, in [START_REF] Riviere | Adaptive canceling of physiological tremor for improved precision in microsurgery[END_REF] an adaptive filter was used to estimate tremor modeled as an harmonic model. The same algorithm was also applied in [START_REF] Rocon | Design and validation of a rehabilitation robotic exoskeleton for tremor assessment and suppression[END_REF]. In that paper, however, explicit estimation of voluntary motion was added.

In the current paper, a new technique to perform tremor and voluntary motion estimation is presented. The algorithm was designed to improve the filter robustness to multiple motions and different sensors. To accomplish that, models are assumed for both tremor and intentional motion and an extended Kalman filter (EKF) is used as the recursive estimation algorithm. The EKF presents the advantage of explicitly considering the measurement noise and concurrently performing estimation of the tremor and voluntary motion.

The assumption of an explicit voluntary motion model may also be important in the estimation of a particular motion or task in order to trigger a compensation pattern. It may also be an useful feature for other applications. For instance, a good estimation of intentional motion may be important in gesture recognition or different rehabilitation therapies involving tremor patients.

The paper is organized as follows. In the following section our EKF-based solution to the problem is described, including description of tremor and voluntary motion models, the estimation algorithm and parameters selection. Section III presents the experimental evaluation of the proposed method with two different data sets. The experimental data was obtained from tremor patients performing the same task, but measured with different sensors. This section also presents a comparison with the results provided by the techniques proposed in [START_REF] Riviere | Adaptive canceling of physiological tremor for improved precision in microsurgery[END_REF], [START_REF] Rocon | Design and validation of a rehabilitation robotic exoskeleton for tremor assessment and suppression[END_REF]. The last section presents the final remarks.

II. Estimating tremor and voluntary motion

It is considered in this paper that the motion performed by the patient is composed by pathological tremor and voluntary motion. Both signals are nonstationary. For the scalar case, the total motion is given by

s(k) = s t (k) + s i (k), ( 1
)
where k is a multiple of T , the sampling period, s t is the tremor component, and s i is the intentional motion. s is measured by a motion sensor, which may be an inertial sensor, an optical tracking system, digitizing tablets or similar devices. The sensor provides the measurement

s(k) = s(k) + ν s (k), (2)
where ν s is an additive white Gaussian noise, ν s ∼ N (0, σ 2 s ), that represents sensor error. Our goal in this work is to estimate online both tremor and voluntary motion components from the noisy measurements of a motion sensor. In the following subsection we describe the models adopted for each separate motion. Following, the recursive algorithm used to estimate the referred signal and the tremor parameters is presented.

A. Models

Since we are estimating both voluntary and tremor motions from the measured resultant motion, our knowledge from these different motions must be employed to separate them. We assume that voluntary motion lies mainly within low frequencies (< 2 Hz), while tremor frequencies are higher (roughly between 3 to 15 Hz).

In this scenario, one may simply choose to apply classic linear filters, but this option was not considered due to two main reasons. Firstly, the intrinsic delays added are not desirable, since the estimated tremor motion will be used by a closed-loop tremor compensation system. Furthermore, since we are interested in estimating in real-time the tremor model parameters, a model-based estimation procedure must be applied in any case.

The other alternative is to describe the motions of interest by general mathematical models, which then are used in the estimation procedure. In [START_REF] Riviere | Adaptive canceling of physiological tremor for improved precision in microsurgery[END_REF], for instance, an harmonic model is assumed for tremor, while no explicit model for intentional motion is assumed. In [START_REF] Rocon | Design and validation of a rehabilitation robotic exoskeleton for tremor assessment and suppression[END_REF], models for both motions are assumed, which is also the chosen approach in this work.

Regarding the tremor model, in [START_REF] Bó | Online pathological tremor characterization with Kalman filtering[END_REF], we have compared the use of Auto-Regressive (AR) and harmonic models to the problem of tremor online estimation only. The cited paper did not discard the use of a particular model, since both models presented good performance for tremor estimation, including the estimation of power spectrum and k-step ahead prediction. For the present case, however, since two motions are estimated simultaneously, harmonic models were chosen as the tremor model. This was mainly due to the fact that the parameters that describe the frequency in harmonic models may be explicitly bounded, while in AR models the frequency behavior is described by the combination of the different parameters. Hence, filter convergence to consistent parameters become less robust.

An harmonic model may represent any periodic signal. In the case of quasi-periodic signals, an approximation may be achieved with harmonic models. Since a nonstationary signal is considered, the following rectangular model is adopted for unidimensional motion:

s t (k) = H h=1 a h (k) sin h k t=1 ω(t) + (3) b h (k) cos h k t=1 ω(t) + ν st (k),
where ω is the fundamental frequency, a h and b h are the coefficients and H is the number of harmonics, the model order. ν st , an additive white Gaussian noise, ν st ∼ N (0, σ 2 st ), represents modeling errors. For motions in ℜ 2 or ℜ 3 , another advantage of harmonic models is that coupling between each coordinate may be represented by assuming that every tremor motion, modeled by (3), has the same fundamental frequency ω.

In this model, it is assumed that tremor may be described by harmonic frequencies, which is often the case. For motions that do not present this feature, one option is to adopt nonharmonic relations between the frequencies, as in [START_REF] Veluvolu | Bandlimited multiple fourier linear combiner for real-time tremor compensation[END_REF]. The algorithm presented in this paper may be easily adapted to this nonharmonic model.

Concerning the voluntary motion, although it is a slower movement, it does not present the regular features of the tremor motion. Hence, in this work it was modeled as low-pass filtered white noise:

s i (k) = F f =1 c f ν si (k -f ), (4)
where ν st is a white Gaussian noise, ν si ∼ N (0, σ 2 si ), and F is the model order. c f , the filter parameters, are fixed parameters tuned to represent the frequency behavior assumed for voluntary motion, i.e., a low-pass filter with cutoff frequency around 2 Hz.

B. EKF

Some algorithms were already proposed to perform the separation of tremor and voluntary motion. The Weighted Fourier Linear Combined (WFLC), described in [START_REF] Riviere | Adaptive canceling of physiological tremor for improved precision in microsurgery[END_REF], is an algorithm designed to estimate an harmonic tremor model, but may be also used to separate tremor and voluntary motion. The estimation problem is solved in two steps: the tremor frequency is estimated from the bandwidth filtered signal and the other parameters are estimated from the raw signal, s. A constant term is added in the harmonic model in order to represent voluntary motion. This solution is highly dependent on the chosen gains of the recursive estimation algorithm for different voluntary motions and sensor noises.

In order to address these problems, in [START_REF] Rocon | Design and validation of a rehabilitation robotic exoskeleton for tremor assessment and suppression[END_REF] another solution was proposed. Voluntary motion is firstly estimated from the original signal by a Benedict-Bordner filter (BBF), which provides ŝi . After, the WFLC algorithm is applied to the signal s -ŝi . The BBF is equivalent to the Kalman filter steady-state solution to a linear first-order system, but with fixed restrictions between the filter parameters. In this case, considering that voluntary motion generally does not behave as the assumed model and that the correction is given by the total measured motion, the result is similar to apply a low-pass filter to estimate ŝi , and then apply the WFLC algorithm to estimate ŝt .

In our approach, we tried to address the disadvantages of previous solutions. Hence, not only models for both motions are assumed, but also the models uncertainties are explicitly considered. Furthermore, measurement error is represented, resulting in a more robust solution. The optimal estimator for such problems is the Kalman filter (KF) for the linear case. Since the present problem is nonlinear, one alternative is to use a modification of the KF for nonlinear systems, the EKF, where Kalman equations are applied to the first-order linearization of the nonlinear system [START_REF] Simon | Optimal State Estimation: Kalman, H Infinity and Nonlinear Approaches[END_REF].

The applied filter simultaneously estimates both referred motions. The filter states are composed by those states related to tremor, x t , which are composed by the estimated motion and the tremor model recursively identified parameters,

s t (k) a 1 (k) • • • a H (k) b 1 (k) • • • b H (k) ω(k)
T and the states related to voluntary motion, x i , which are organized accordingly,

s i (k) • • • s i (k -f ) ν si (k -1) • • • ν si (k -f ) T .
To obtain the full state vector, we combine the states regarding tremor and voluntary motion:

x(k) = x t (k) x i (k) . (5)
The total number of states will be defined by the orders chosen for both motion models:

n = 1 + (2H + 1) + 2F. ( 6
)
In the KF framework, estimation is organized in prediction and correction phases. In the prediction phase, x-(k), the initial a priori estimate is computed according to

x-(k) = f (x(k -1)), ( 7
)
where f (x) is composed by Eqs. ( 3) and ( 4) and the parameters models. All parameters are modeled as random walk processes, particularly since they do present fast dynamics. For instance, the evolution of parameter ω is given by

ω(k) = ω(k -1) + ν ω (k), ( 8
)
where ν ω is a white Gaussian noise, ν ω ∼ N (0, σ 2 ω ). The error covariance matrix, P, is propagated according to

P -(k) = AP(k -1)A T + Q, ( 9
)
where A is the jacobian matrix of f (x(k-1)) with respect to x(k -1) and Q, the process covariance matrix, is a diagonal matrix whose components are the variances of each filter state.

In the correction phase, the available measurement,

y(k) = g(x(k)), (10)
where g(x(k)) is given by Eqs. ( 1) and ( 2), is used to correct the initial estimate, providing x+ (k), the final a posteriori estimate. The Kalman gain, K, scales the correction applied in the correction phase:

K(k) = P -(k)C T CP -(k)C T + r -1 , ( 11
)
where C is the jacobian matrix of g(x -(k)) with respect to x-(k). Once it is computed, the final estimate is computed according to

x+ (k) = x-(k) + K(k) y(k) -g(x -(k)) (12)
and the corresponding error covariance matrix is given by

P(k) = (I -K(k)C) P -(k), ( 13
)
where I is the identity matrix.

Regarding the parameters that may be tuned in the algorithm, another advantage of the EKF is that the filter parameters have a direct physical meaning. They are related to the associated uncertainties and to the initial conditions. The process covariance matrix, Q, represents the uncertainties related to the motion models and the models parameters. Variances related to the models are lower for imprecise models, while variances related to parameters are proportional to the time-varying nature of those parameters. As to the measurement variance r, it is considered equal to the variance σ 2 s . Concerning the initial conditions, the initial frequency estimate is critical. Fortunately, to obtain a good estimate is not an unusual procedure, being normally conducted in clinical settings. The other parameters are initialized as zero. The initial error covariance matrix, P, have similar values of corresponding states variances.

Angular rate sensor

Digitizing tablet

III. Experimental evaluation

In this section, we present an experimental evaluation of the proposed algorithm. The task of drawing a spiral was chosen. It is one of the common methods to evaluate tremor severity [1], particularly because loosing his/her writing ability is one of the most serious functional disabilities associated with pathological tremor.

Figure 1 illustrates the experimental setup, where both hand tremor motion and writing data were measured. Commercial sensors were used for both acquisitions. For the inertial measurements, we employed the IDG300, an analog device from Invensense. To record writing data, the Intuos3, a digitizing tablet with 5, 080 lpi resolution from Wacom, was used.

Two different data sets from two patients are analyzed in this paper. In order to evaluate the algorithms performance under different conditions, while in one case the signal from a angular rate sensor was used, in the other data from a digitizing tablet was employed. In figure 2 the power spectrum densities (PSD) of both data sets are shown, which illustrates the assumptions made about the motions frequency ranges. In the analysis conducted within the current section, both data sets were normalized by initial estimates of signal offset. For the presented results, though, the data is de-normalized.

The obtained results are briefly compared with the results provided by the two other cited techniques, the WFLC and the combination of BBF+WFLC. Since a perfect reference voluntary motion is difficult to measure, if not impossible, a quantitative comparison of the methods is difficult to be performed. Hence, the comparison between evaluated methods was made on a visual basis.

Concerning the parameters applied in the algorithms, the same parameters were used for both data sets. The parameters used for the EKF-based algorithm proposed in this paper are shown in the Table I. Voluntary motion has been modeled as a digital low-pass filter with cutoff frequency on 2 Hz. The coefficients are not listed on the table, since they depend on the sampling frequency.

Regarding the method proposed in [START_REF] Riviere | Adaptive canceling of physiological tremor for improved precision in microsurgery[END_REF], the same published parameters were used. Parameters used in [START_REF] Rocon | Design and validation of a rehabilitation robotic exoskeleton for tremor assessment and suppression[END_REF] were not found. Hence, an effort was made to find the parameter g (in [START_REF] Rocon | Design and validation of a rehabilitation robotic exoskeleton for tremor assessment and suppression[END_REF]) that fitted both data sets, since the algorithm performance was highly dependent on the data. The value chosen for g was 0.015 and the WFLC parameters were the same taken from [START_REF] Riviere | Adaptive canceling of physiological tremor for improved precision in microsurgery[END_REF].

A. Angular rate sensor

Using a low-cost angular rate sensor attached to a patient's hand, tremor angular velocity was measured while the referred task was performed by the patient. The patient was diagnosed with Parkinson's Disease (PD) and presented variable tremor with unstable amplitude. The sensor was used to measure flexion/extension of the wrist, which was the most severe tremor direction. Tremor frequency was around 7.5 Hz.

Figure 3 shows the acquired data and the estimated voluntary motions. The results provided by the WFLC and the EKF are similar, although the voluntary motion estimated by the WFLC is more noisy, since it is computed directly by s -ŝt . The algorithm proposed in [START_REF] Rocon | Design and validation of a rehabilitation robotic exoskeleton for tremor assessment and suppression[END_REF] was not able to filter effectively the tremor from intentional motion with the gains applied in our work. Smaller g provided better results, but introducing considerable delay if too small gains were applied.

B. Digitizing tablet

A digitizing tablet is a device commonly used in the clinical evaluation of pathological tremor. It provides two-dimensional measurements of an electronic pen. In this experiment, the patient was told to draw a spiral following a reference drawing. However, we did not assume that the reference spiral could be used as the ground truth for voluntary motion. The data presented here is from a patient diagnosed with essential tremor (ET), who had an almost disabling tremor on the dominant hand.

Tremor frequency was around 6.8 Hz.

Figure 4 shows both estimated voluntary motions on coordinates X and Y . In this case, the differences in performance were greater if compared to results presented on the previous subsection. The WFLC algorithm was not able to correctly estimate tremor, as illustrated by Fig. 5, and the error was propagated to voluntary motion estimation. Better results were obtained with smaller gains to the FLC phase of the algorithm. With respect to the method BBF+WFLC, the voluntary motion is obtained first from the original signal, s. From Fig. 4, it is shown that the result was excessively smoothed and delayed. Better performance was found with higher values for the g parameter. Regarding the results obtained from our EKF-based algorithm, it is illustrated on Fig. 4, that the correct drawing without tremor components is obtained, and on Fig. 5 that voluntary motion was correctly filtered.

C. Discussion

In order to analyze the obtained results, it is of great importance to remember that the purpose of this algorithm is to estimate both tremor and voluntary motion from a signal that contains both components. Hence, the algorithm robustness to different motions is an essential feature. In addition, its robustness to signals obtained from different sensors is desirable. In order to compare those features, in both experiments all algorithm parameters in all methods were kept unchanged. The only exception to this rule is the frequency initial estimate.

Within this evaluation framework, while the performance of the EKF-based algorithm was accurate and insensitive to different data sets, this was not the case for the WFLC and BBF+WFLC methods. The same task was executed in the two experiments, but with different sensors sampled with different frequencies. For a particular data set, the parameters of each algorithm could be tuned to provide better results than those presented on Sections III-A and III-B, but it would also decrease the performance in the other experiment, revealing the low robustness of those algorithms. The insensitiveness of the proposed algorithm was greatly due to the assumption of Eq. ( 4), as the voluntary motion model, and the explicit consideration of process and measurement uncertainties. As for the WFLC and the BBF+WFLC, it was noted that they are not only sensitive to different voluntary motion dynamics and amplitude ratio between tremor and intentional motion, but also to different sampling frequencies. In our comparison, however, care was taken to adjust the algorithm for the right sampling frequency. For the WFLC, for instance, the Butterworth filter applied during the processing was designed accordingly.

The performance of our algorithm with different parameters was not deeply evaluated. The goal was to obtain a sufficiently good performance and compare it with the other methods. Whenever possible, parameters were of the same order from those used in the other methods. H, for instance, was set to 1, since it was the value used in [START_REF] Riviere | Adaptive canceling of physiological tremor for improved precision in microsurgery[END_REF]. Also, the algorithms performance in tracking tremor frequency was not compared, since this problem has been mainly addressed in [START_REF] Bó | Online pathological tremor characterization with Kalman filtering[END_REF]. One possible disadvantage of the presented algorithm is that it is more complex than the other methods. However, considering currently available computational units and the scalar corrections are applied, it is a negligible issue.

Finally, it is important to point out the potential consequences of errors in the estimated tremor and voluntary motions. The applications of the provided estimated motions (outputs of the algorithm presented in this paper) may be used in two classes of problems. Firstly, the estimated intentional motion may be used in applications such as filtering tremor from writing, gesture recognition. In those cases, the overestimation of intentional motion or the existence of high frequency components may prevent the practical use of the measured data. In a second scenario, the estimated tremor motion may be used as an input to an active tremor compensation system. In this case, a typical application uses an estimate of tremor frequency and amplitude to regulate compensatory action. If tremor is modeled as a first-order harmonic model, tremor amplitude may be given by

P t = b 1 -ia 1 2 4 .
Within this context, errors in tremor amplitude estimation from tremor data such as those presented in Fig. 5 may also prevent the practical use of such systems.

IV. Conclusions

For pathological tremor compensation applications, such as a FES-based system, estimation of both tremor and voluntary motions from the resultant composed measured motion is of great importance. In this paper, an algorithm to address this problem is presented. Models are assumed for both motions and the estimation is performed with an EKF.

The performance of the proposed method was evaluated with two experimental data sets from tremor patients performing the same task, drawing a spiral. The motion was measured with two sensors (angular rate and digitizing tablet) and the results were compared with two other methods described in the literature [START_REF] Riviere | Adaptive canceling of physiological tremor for improved precision in microsurgery[END_REF], [START_REF] Rocon | Design and validation of a rehabilitation robotic exoskeleton for tremor assessment and suppression[END_REF]. The proposed method presented a better performance, particularly in terms of robustness to different data sets, being a good candidate for application in real pathological tremor compensation systems.

Fig. 1 .

 1 Fig. 1. Graphical representation of the experimental procedure. While the patient performed the task, both the hand motion and the resulted writing were recorder.

Fig. 2 .

 2 Fig. 2. PSD of both data sets used in this work. For the tablet data, position and velocity in X are shown.

Fig. 3 .Fig. 4 .

 34 Fig.3. Angular velocity data (black) and the voluntary motion estimated by our method (green) and the two other evaluated algorithms.

Fig. 5 .

 5 Fig. 5. Tremor in coordinate X estimated by all three evaluated methods.

TABLE I

 I Parameters used for our EKF-based algorithm.

	H	σ 2 s t	σ 2 ω	σ 2 a h ,b h	F	σ 2 s i
	1	5.10 -8	5.10 -4	5.10 -3	2	5.10 -2

Antônio P. L. Bó (e-mail: bo@lirmm.fr) and Philippe Poignet are with the LIRMM UMR 5506 CNRS UM2, France. Antônio P. L. Bó is also with Neuromedics, France. Christian Geny is with the CHU Montpellier, France.


Marc Chaumont
email: marc.chaumont@lirmm.fr

Psychovisual Rotation-Based DPTC Watermarking Scheme

The Dirty Paper Trellis Codes (DPTC) watermarking is a very efficient high capacity scheme. Nevertheless, it suffers of two drawbacks: its security weakness and its CPU computation complexity. In a previous work, we treat those drawbacks by proposing a more secure embedding space and a faster embedding technique. In the current paper, we study the robustness of the scheme with the use of a psychovisual mask and a correcting code. First, we explain how to integrate the psychovisual shaping inside of the DPTC watermarking scheme. Second, we treat of the correcting code integration. Two different masks are used, and robustness comparisons are drawn between the previous DPTC and the two protected psychovisual DPTC. One can observe that for a given SSIM degradation, the robustness results are better for the two protected psychovisual DPTC approaches.

INTRODUCTION

The generation of watermarking schemes named informed schemes or side information schemes appeared around 1998 when Costa's work has been rediscovered [START_REF] Costa | Writing on dirty paper[END_REF]. The two principal techniques' categories for multi-bits watermarking are the lattice codes also named quantized based watermarking schemes (DC-QIM [START_REF] Chen | Quantization index modulation: A class of provably good methods for digital watermarking and information embedding[END_REF] and SCS [START_REF] Eggers | Scalar Costa Scheme for Information Embedding[END_REF]) and the Dirty Paper Trellis Codes (DPTC) [START_REF] Miller | Applying Informed Coding and Informed Embedding to Design a Robust, High Capacity Watermark[END_REF].

The original DPTC algorithm is known for its good robustness and its high embedding capacity but owns two strong weaknesses: the informed embedding step uses a Monte Carlo approach which is very CPU time consuming and the scheme own security weakness facing collusion attack [START_REF] Bas | Practical Security Analysis of Dirty Paper Trellis Watermarking[END_REF]. In this paper we present a psychovisually secure, robust and rapid approach.

In a previous work [START_REF] Chaumont | A Novel Embedding Technic for Dirty Paper Trellis Watermarking[END_REF], we give a solution, named rotationbased DPTC (RB-DPTC), in order to embed rapidly and securely. In Section 2, we remind its principle. In this paper, we propose to improve RB-DPTC by using a psychovisual mask and a correcting code. We thus expose in Section 3 the psychovisual shaping and briefly present the used correcting code. Those improvements lead to two new watermarking schemes: the protected with rudimentary's mask RB-DPTC (PR-RB-DPTC) and the protected with Xie and Shen's mask RB-DPTC (PXS-RB-DPTC). Finally, in Section 4, we compare the three approaches with a fixed robustness criteria and with a fixed psychovisual metric.

THE ORIGINAL ROTATION-BASED DIRTY PAPER TRELLIS CODES (RD-DPTC)

In this section, we present the embedding space and the embedding approach.

Embedding space

The embedding space is obtained by first, a wavelet transform of the image, and second, a projection of the host signal x of dimension N wlt (x is the concatenation of sub-bands coefficients except LL sub-band's coefficients) onto N sec carriers of same dimension.

Carriers are noted

u i with i ∈ [1, N sec ].
Note that a projection is just a scalar product. The obtained vector v x may then be used for the informed-coding and informed-embedding (see Section 2.2).

The carriers u i are built from normalized bipolar pseudorandom sequences. For computational complexity reasons, carriers are nor orthonormalized neither drawn from a Gaussian distribution. This is not a weakness since in high dimension, carriers are orthogonal and Gaussian property is not essential. Moreover, in order to reduce projections complexity to a linear complexity (and even with better robustness) we are using a Space Division Multiplexing approach [START_REF] Chaumont | A Fast Embedding Technique For Dirty Paper Trellis Watermarking[END_REF]. It consist in first, shuffling the wavelet coefficients of x and second, dividing x in quasi-equal region' size and attribute a carrier by region.

Informed-coding and informed-embedding

Once that the projection of the host vector x onto carriers u i , i ∈ [1, N sec ] are proceeded we obtain the host vector v x . We then run the informed-coding which is the same as the original one [START_REF] Miller | Applying Informed Coding and Informed Embedding to Design a Robust, High Capacity Watermark[END_REF]. The informed-coding takes as input the host vector v x and the message m, to be embedded, and return a codeword c * . This vector c * (of size N sec ) is the closest one to v x among vectors coming from the codebook C , and representing the m message. For more details see [START_REF] Miller | Applying Informed Coding and Informed Embedding to Design a Robust, High Capacity Watermark[END_REF] or [START_REF] Chaumont | A Novel Embedding Technic for Dirty Paper Trellis Watermarking[END_REF].

After proceeding to the informed-coding, the codeword c * is extracted. The objective of the informed-embedding is to push the host vector v x into the Voronoï region of c * in order to obtain the watermarked vector v y . Our approach consists in dichotomously reducing the angle between the host vector v x and the codeword c * until obtaining the smallest absolute angle regarding all the other angles; this angle is noted θ f . Then, one penetrate inside the Voronoï region with a given angle θ R . Our informed embedding is thus a rotation of v x of an oriented angle equal to max(θ f + θ R , ( v x , c * )). The marked vector v y is then obtained.

One then compute the watermark vector v w = v y -v x , retroproject it onto carriers in order to obtain the watermark signal w and then compute the watermarked signal y = x + w. The inverse wavelet transform of y gives the watermarked image. At the extraction we project wavelet coefficients onto secret carriers and then retrieve the closest codeword (and thus the message) from the codebook C .

PSYCHOVISUAL SPACE AND CORRECTING CODE

General scheme with psychovisual shaping

Psychovisual masks may be used in order to shape the watermark signal. Roughly speaking, the watermark signal power will be reduced in the flat regions and increased in the textured or contours regions. Figure 1 shows the general embedding scheme with the use of a psychovisual mask. Referring to that figure, there is three majors steps: 1. the construction of a psychovisual space x psy . In this space, coefficients are psychovisually equivalents. One thus, can fairly embed. This psychovisual space is such that:

∀i ∈ [1, N wlt ], x psy [i] = x[i]/α[i],
with α the psychovisual mask;

: ∀i ∈ [1, N wlt ], w[i] = w psy [i] × α[i]
. This should reduce degradation, due to watermarking, in the visible regions. For example, flat regions should own a small α value in order to reduce the power of w in those regions. 3. the shaped watermark embedding. This embedding is such that:

∀i ∈ [1, N wlt ], y[i] = w[i] + x[i].
During the decoding, the decoder should re-compute the psychovisual mask. Note that this is a sensitive part since the recomputed mask should be the same than the one used during embedding. The decoder extract the wavelet (watermarked and attacked) signal, divide it component-wise by α, project coefficients onto secret carriers and then retrieve the closest codeword (and thus the message) from the codebook C .

Rudimentary psychovisual's mask

Before using a psychovisual mask issued from the literature, we experiment the watermarking scheme behavior, if using a rudimentary mask. Results show that with a simple approach, the psychovisual impact is really good. The mask is obtained by, first, filtering the image with a sort of high-pass filter, second, apply a wavelet transform to the filtered image, and select the wavelet coefficients, and third, linearly scale their absolute value between two predefined scalar α min and α max1 . The four steps are detailed above:

1. ∀i ∈ [1, N], I f iltered [i] = 0.1 + 10 * |I[i] -I gauss

Xie and Shen's Pixel-Wise Mask

In [START_REF] Xie | Toward Improved Wavelet-Based Watermarking Using The Pixel-Wise Masking Model[END_REF], Xie and Shen propose an improvement of the well known pixel-wise masking model of Barni et al. [START_REF] Barni | Improved Wavelet-Based Watermarking Through Pixel-Wise Masking[END_REF]. The pixel-wise mask (PW-M) predicts the Human Visual Systems response to disturbance at a given coefficient position in the wavelet domain.

Let's note I θ l a sub-band at resolution level l ∈ [1, L], with orientation θ ∈ {a, h, v, d}. The PW-M value for a position (i, j) for a sub-band at level l is defined as the product of three terms:

α θ l (i, j) = Θ(l, θ ).Λ(l, i, j).Ξ(l, i, j) 0.2 ,
where the three terms are used to evaluate the sensitivity to noise changes depending on the band, the local brightness, and the local texture activities, respectively. The noise sensitivity function is computed as the following:

Θ(l, θ ) = √ 2 , i f θ = d 1, otherwise .      1.00, i f l = 1 0.32, i f l = 2 0.16, i f l = 3 ...     
.

The local brightness is computed by taking the reconstructed approximation I a l (i, j) sub-band at the level, l, of interest2 , as the following:

Λ(l, i, j) = 1 + L (l, i, j), with L (l, i, j) = 1 -L(l, i, j), i f L(l, i, j) < 0.5. L(l, i, j), otherwise.
with L(l, i, j) = I a l (i, j) 255 .

The texture activity in the neighborhood of a coefficient, is computed as the following:

Ξ(l, i, j) = [I h l (i, j)] 2 + [I v l (i, j)] 2 + [I d l (i, j)] 2 3 .VAR{I a l (i, j)}.
where the first contribution of the product is the local mean square value of the DWT coefficients at the same resolution level, l, and the second contribution, VAR{I a l (i, j)}, is the local variance of the corresponding reconstructed approximation sub-band in a small neighborhood (3 × 3 window) corresponding to the location (i, j).

In order to compute α vector, we scale the values of each α θ l (i, j) in the range [α min , α max ] as explain in section 3.

Error Correcting Code

In order to add more robustness to the transmitted message m, one encode it with a very simple convolutional code 2-memory 1/2-rate.

The state machine of this code is given in Figure 3.a. When an input bit enter into the coder, it causes a transition to a new state and outputs two bits. The two outputs bits are the resulting coded bits. The state machine can also be represented as it evolves in time with a trellis diagram. Figure 3.b shows the trellis associated to the state machine of Figure 3.a. Usually, a trellis is constructed by placing all the states in column and each possible transition are drawn with an arc between states at t time and states at t + 1 time. By convention, bold arcs represent the 1 input and non-bold arcs the 0 input.

The convolutional coder takes a sequence of bits in input and generates an output sequence thanks to the state machine. The output sequence form a codeword representing the message m and is embedded into the image thanks to our rotation-based DPTC. The decoder, when receiving the possibly degraded codeword, finds the closest codeword and returns it. In order to find this closest codeword, the Viterbi algorithm is used [START_REF] Viterbi | CDMA: Principles of Spread Spectrum Communication[END_REF].

RESULTS

Tests have been proceeded on the first 100 images of the BOWS2 data-base3 with images resized to 256 × 256 4 . Those images are 8-bits grey-level images and are personal photos.

The trellis structure owns 128 states with 128 arcs by states. Outputs arcs labels are drawn from a Gaussian distribution and there are N arc = 12 coefficients by output arc for the RB-DPTC approach and N arc = 10 coefficients by output arc for PR-RB-DPTC and PXS-RB-DPTC. The used payload is payload = 1 bit of the (non-coded) message for 64 pixels which is the same as the original DPTC algorithm [START_REF] Miller | Applying Informed Coding and Informed Embedding to Design a Robust, High Capacity Watermark[END_REF]. The number of embedded bits is thus 1024 bits for RB-DPTC and 2048 for the two improved RB-DPTC (since correcting code rate is rateECC = 1/2). Wavelet transform is a 9/7 Daubechies with l = 3 decomposition levels. Except the LL subband, all the other sub-bands are used to form the host signal x. With 256 × 256 images, the wavelet space size is thus N wlt = 64 512 coefficients. For the RB-DPTC, there are 1024 bits to embed and the number of outputs coefficients for a trellis' arc N arc = 12 coefficients. The private space size is thus N sec = 1024 × 12 = 12 288. For the two improved RB-DPTC there are 2048 bits to embed and the number of outputs coefficients for a trellis' arc is N arc = 10 coefficients. The private space size is thus N sec = 2048 × 10 = 20 480.

Lets remark that for the RB-DPTC, one bit of the message influences N arc × N wlt /N sec = 63 wavelet coefficients and with the two improved RB-DPTC, one bit of the message influences N arc ×N wlt /N sec /rateECC = 63 wavelet coefficients. Thus, the two major differences are: 1. the spreading ratio N wlt /N sec . The spreading is higher for RB-DPTC (ratio = 1 coefficient of the secure space is spread on 5.25 coefficients of the wavelet space) than for the improved RB-DPTC (ratio = 1 coefficient of of the secure space is spread on 3.15 coefficients of the wavelet space); 2. the number of outputs coefficients for a trellis' arc. The RD-DPTC (N arc = 12) should be more robust than the two improved RD-DP (N arc = 10). A contrario, the "adds" of a correcting codes for the two improved RD-DPTC improves the robustness to error and the "add" of a psychovisual mask should allows to increase watermarking strength. Four kinds of robustness attacks have been proceeded: Gaussian noise attack, Gaussian filtering attack, valumetric attack and jpeg attack. The Bit Error Rate (BER) is the number of erroneous extracted bits divide by the total number of embedded bits. The BER is computed for each attack. Two protocols are used in order to evaluate the robustness. For the first protocol, the inside angle penetration (which stand for the robustness parameter) is fix in order to obtain a quasi-equivalent PSNR degradation and for the second protocol, we choose to fix the SSIM degradation to 98% [START_REF] Wang | Image Quality Assessment: From Error Visibility to Structural Similarity[END_REF] 5 .

Three algorithms are competing with the two different evaluation protocols:

• the RB-DPTC detailed in [START_REF] Chaumont | A Novel Embedding Technic for Dirty Paper Trellis Watermarking[END_REF]. For the first protocol, the inside angle penetration is fixed to θ R = 0.11 radian and we are obtaining a mean embedding PSNR of 42.15 dB. • the PR-RB-DPTC using the rudimentary mask (see Section 3.2) and the convolutional code 2-memory and 1/2 rate (see Section 3.4). For the first protocol, the inside angle penetration is fixed to θ R = 0.09 radian and we are obtaining a mean embedding PSNR of 41.75 dB. • the PXS-RB-DPTC using the Xie and Shen's mask (see Section 3.3) and the convolutional code 2-memory and 1/2 rate (see Section 3.4). For the first protocol, the inside angle penetration is fixed to θ R = 0.001 radian and the mean embedding PSNR is 42.38 dB.

For the first evaluation protocol, we can observes on Figures 4.a, 4.b, 4.c and 4.d that below 1% BER, it is better to use the PR-RB-DPTC approach. Indeed, for very low attack power, comparing to the two other approaches, there are no bit errors. A contrario, at higher attack power, the BER is stronger for PR-RB-DPTC and PXS-RB-DPTC than for the original RB-DPTC. This could be explained by the fact that we are using correcting code in the two new approaches. When errors are too numerous, the performances of the correcting code collapse. We should nevertheless note that their is a very strong variance of PSNR over the 100 images: RB-DPTC PSNR ∈[35.4, 61.7]; PR-RB-DPTC PSNR ∈ [29.9, 56.9]; PXS-RB-DPTC PSNR ∈ [33.2, 60.8]. Those high fluctuations are too strong and may lead to wrong interpretation. Since psychovisual mask have been used, in order to better compare the three watermarking schemes, we decide to fix the degradation criteria to a SSIM value of 98%.

For the second protocol, for each scheme, we dichotomically approach the SSIM 98% value by modifying the penetration angle 6 . Figure 5.a, 5.b, 5.c, 5.d, illustrate the BER for the different attacks with a SSIM = 98% embedding. Note that a SSIM of 98% gives images of very good quality. We can observe that for low power attacks, the two improved approach (PR-RB-DPTC and PXS-RB-DPTC) give better results than the original RB-DPTC (there is an improvement of the BER of a factor up to 10% for gaussian, filtering and scaling attacks). For those low power attacks, the approach using the rudimentary mask (PR-RB-DPTC) is a little bit better than the one using the Xie and Shen's mask (PXS-RB-DPTC). For higher power attack, except for the Gaussian attack, the results are comparable for the three approaches. The weak robustness to jpeg attack was already observed in previous work of [START_REF] Chaumont | A Novel Embedding Technic for Dirty Paper Trellis Watermarking[END_REF] and [START_REF] Chaumont | A Fast Embedding Technique For Dirty Paper Trellis Watermarking[END_REF] and is inherent to the approach since the host vector x is in part, made up of high frequency wavelet coefficients. We can remark that the same sensitivity to jpeg attack is observed with quantized-based approaches when all the DCT coefficients are used [START_REF] Li | Using Perceptual Models to Improve Fidelity and Provide Resistance to Valumetric Scaling for Quantization Index Modulation Watermarking[END_REF].

CONCLUSION

In this paper, we propose improvements to a previous work on DPTC. Our previous approach named rotation-based DPTC (noted RB-DPTC), treated two major problems of security and complexity for DPTC. In this paper, we essentially improve psychovisual degradation due to watermark embedding and add a correcting code in order to reduce the BER due to the use of a psychovisual mask. More precisely, we explain how to integrate a psychovisual mask inside of the original RB-DPTC, we briefly present the used correcting code and explain how to modify the trellis parameters. Results show that for a fix SSIM value, there is a strong improvement (of a factor 10) at low power attack. Future work should investigate: the link between SSIM metric and penetration angle, the construction of robust psychovisual masks, the use of RDM concept, the robustness to jpeg attack and sensitivity to regression attack.

Figure 1 :

 1 Figure 1: Embedding scheme with a psychovisual mask

  [i]| with I gauss the result of a Gaussian filtering of image I. 2. Apply a wavelet transform to I f iltered and order coefficients in a vector β . 3. Compute β min and β max such that β min = min i∈[1,N wlt ] |β [i]| and β max = max i∈[1,N wlt ] |β [i]| and then compute α by scaling absolute values of β between [α min , α max ] (knowing β min and β max ). The rudimentary mask for the Barbara 8-bit image (Figure 2.a), crop to a 512×512 image, is given in Figure 2.c. with pixels values multiply by 255. The Figure 2.b gives the 3 levels wavelet 9/7 Daubechies decomposition centered on the 127 value.

  2. Figures 2.d illustrate the Xie and Shen mask apply on Barbara image.

Figure 2 :Figure 3 :

 23 Figure 2: Illustration of psychovisual's masks

Figure 4 :

 4 Figure 4: BER for attacks on RB-DPTC, PR-RB-DPTC and PXS-RB-DPTC. Embedding is proceeded with a fix penetration angle.

Figure 5 :

 5 Figure 5: BER for attacks on RB-DPTC, PR-RB-DPTC and PXS-RB-DPTC. Embedding is proceeded with a with a fix SSIM of 98%.

α min = 1 and α max = 8 for PR-RB-DPTC and α min = 1 and α max = 4 PXS-RB-DPTC.

In our implementation, the I a l images are build from a Gaussian pyramid.

BOWS2 data-base is located at http://bows2.gipsa-lab.inpg.fr/.

The images sub-sampling has been proceeded with xnview program and using Lanczos interpolation.

SSIM is a classical measure well correlated to the Human Visual System. The SSIM values are real positive numbers lower or equal to 1. Stronger is the degradation and lower is the SSIM measure. A SSIM value of 1 means that the image is not degrade. Lets remark that we are using the C++ implementation by Mehdi Rabah at http://mehdi.rabah.free.fr/SSIM/.

For few images, SSIM may not be tune to the 98% value depending on the original angle between v x and c * , and the angle of the Voronoï region.

Acknowledgment

This work is in part supported by the VOODDO (2008-2011) project of the french ANR (Agence Nationale pour la Recherche).