Julien Garcia

Thierry Laroche

William Daniau

Emile Carry

Gilles Martin

S Ballandras
email: sylvain.ballandras@femto-st.fr

J.-M Friedt
email: jmfriedt@femto-st.fr

Jean-Michel Friedt

temperature effects on periodic transducers using a FEA/BEM approach

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.

INTRODUCTION

The control of thermal stability of acoustic wave transducers still is one of the major problems in developing advanced solutions for passive radio-frequency (RF) signal processing. Controlling the thermal frequency drift of such devices enables one to reduce the transition bandwidth of such filters, then increasing their selectivity and allowing for high density spectrum channeling. Different approaches have been tested to efficiently address this problem, yielding complex transducer architectures combining highly coupled single crystal substrates with various overlays. Also new kinds of wave excitation and guiding have been imagined and tested in that purpose, such as interface waves. The same kind of challenge occurs for film bulk acoustic resonators (FBARs), the very serious alternative to SAW when addressing telecommunication standard frequency rising for instance, for which a combination of positive and negative thermal drift material allows for reducing the overall temperature sensitivity of the final devices.

We have been developing a mixed finite element/boundary element modeling code capable to simulate almost any kind of 1D-periodic as well as 2D-periodic transducer structures, accounting for acoustic radiation above and below the waveguide and of course different nature of material losses (viscoelastic, dielectric). In this work, this code has been enhanced to account for the thermal evolution of material characteristics as well as wave-guide thermal expansion to simulate the temperature dependence of guided (or partially guided or even radiated) modes in layered transducer structures. The basis of these developments consists of a material coefficient perturbation according to the celebrated Campbell & Jones approach [START_REF] Campbell | A method for estimating crystals cuts and propagation direction for excitation of piezoelectric surface waves[END_REF], with an appropriate strategy to expanse the meshed region along temperature. In these developments, we neglect thermo-differential stresses. This represents a severe limitation of this approach, but reveals reasonable for achieving reliable predictions of the actual SAW Temperature Coefficients of Frequency (TCF) as demonstrated by Pastureaud & al [START_REF] Pastureaud | Prediction of the thermal sensitivity of surface acoustic waves excited under a periodic grating of electrodes[END_REF]. The principle of the calculation then consists in computing the harmonic admittance of a given transducer at different temperatures and to automatically extract parameters such as resonance and antiresonance frequencies (and also the stop-band width when accessible). The thermal sensitivity of these parameters then is obtained by fitting the results considering quadratic or cubic temperature dependence according to the regarded configuration.

In this paper, we expose the basic theory of our temperaturedependent model, and we show various results validating the proposed approach for bulk acoustic wave (BAW) and SAW devices. We then propose to simulate more complex transducer structures such as thick guiding layers, or even interface wave transducer, to illustrate the capability of the model to address complex guiding structures without restriction.

II. THEORETICAL BASIS

A. FEA-BEM Model

As already reported in many references (the most detailed description of the model will be found in ref. [START_REF] Ballandras | Finite-element analysis of periodic piezoelectric transducers[END_REF]), FEA can be performed for periodic devices with rather simple modifications of the basic algebraic formula associating the displacement and electrical fields to the boundary solicitations. Figure 1 shows the typical configuration for which the model was developed. We assume an inhomogeneous region of the wave-guide treated by finite element techniques above semi-infinite or finite thickness homogeneous regions (the substrate) treated by boundary integral methods. Note that the substrate can be composed of a single material or a stack of layers provided flat interfaces between the involved materials. Considering the harmonic excitation illustrated by [START_REF] Campbell | A method for estimating crystals cuts and propagation direction for excitation of piezoelectric surface waves[END_REF], we relate all the degrees of freedom (dof) on boundary ΓA to those on boundary ΓB, yielding the following expression (1) in which u represents the mechanical displacements and φ the electrical potential. The excitation parameter γ allows to account for the harmonic excitation figure applied to the array (see ref. [START_REF] Ballandras | A mixed finite element/boundary element approach to simulate complex guided elastic wave periodic transducers[END_REF] for a comprehensive definition of this parameter). This relation then is used to reduce the number of independent dof of the FEA model. This is performed without changing the total number of dof of the problem and by simply using a variable change operator C depending on the excitation parameter γ. This provides the following form of the FEA algebraic system to solve

{ u Γ B φ Γ B } = { u Γ A φ Γ A } e -j2πγ
t [ C u  γ  0 0 C φ  γ  ][ K uu -ω 2 M uu K uφ K φu K φφ ][ C u  γ  0 0 C φ  γ  ] { v φ } = t [ C u  γ  0 0 C φ  γ ] { F Q } (2)
where K and M are respectively the FEA stiffness and mass matrices, ν and φ the independent dof of the problem and F and Q the right hand side boundary forces and electrical charges. The upper script * denotes a complex conjugation. Since K can be complex, the left hand side matrix in ( 2) is complex general (hermitian if K∈ℜ), but sparse. These properties are considered when solving the problem. The algebraic system (2) is solved, allowing for the computation of the total electrical charge under active electrodes simply by nodal value summation. The harmonic (frequency and γ dependent) admittance then is equal to the electrode current for a unit voltage excitation (φ0=1V.).

We now regard the case of acoustic radiation on one border of the meshed domain. In that purpose, the general variationnal equation is considered, limited to the purely elastic problem without any loss of generality

∭   ρω 2 u i δu i - ∂ δu i ∂ x j C ijkl ∂ u l ∂ x k dV = ∬ Γ δu i T ij n j dS (3)
in which δui is the variationnal unknown and nj the normal to the boundary Γ (bordering the domain Ώ) on which the radiation boundary condition may partially occur. The mass density is represented by ρ and Cijkl holds for the elastic constants. Equation ( 3) is written in 3D but of course its restriction to 2D problems does not induce any fundamental difficulty. The right hand side of (3) then is considered separately. In this matter, one can relate the dynamic stress Tij to the displacement uk in the spectral domain (denoted by ∼ ) by using a Green's-function-based relation which generalizes the usual surface stress relation widely used in SAW modeling as follows

 T ij =  G ijk s1 ,ω  u k (4)
Note that in this equation, the Green's function asymptotically behaves like s1 when s1 tends to infinity and this property must be carefully accounted for. Equation ( 4) allows one for considering any flat boundary for the application of the radiation conditions, even if tilted in the (x1, x2) plane. Using the now well-established periodic Green's function formalism, the right hand side of ( 4) is expressed as with p the period of the problem and l the number of the current space harmonic. The classical FEA interpolation procedure then is applied to [START_REF] Lichtensteiger | Ferroelectric size effects[END_REF], yielding the final expression of the boundary radiation operator [START_REF] Ballandras | Finite-element analysis of periodic piezoelectric transducers[END_REF].

∬ Γ δu i T ij n j dS = ∬ Γ δu i  x  ∑ l=-∞ ∞ G ijk  γ+l,ω n j p e -j π 2 p  γ+l  x-x'  u k  x' dx'dx (5)
The contribution of the total radiating boundary to the global algebraic system to be solved then consists in a frequency and excitation parameter dependent matrix X(ω,γ) related to both dof and variationnal unknowns and consequently computed in the left hand side of (2). Even for materials without any losses, the final algebraic system is general complex without any specific mathematical characteristics. Also the sparse nature of the algebraic system is degraded due to the connection of all the dof one to the other via the radiation coupling.

B. Thermal perturbation

Temperature effects can be simulated rather efficiently along the so-called Campbell & Jones approach [START_REF] Campbell | A method for estimating crystals cuts and propagation direction for excitation of piezoelectric surface waves[END_REF]. It consists in computing the effective constants of materials once temperature has changed, along a Taylor development. Physical properties of the structure (elasticity, piezoelectricity, dielectric permittivity…) actually change with temperature, and the structure expand in a complicated manner as different materials may be involved. Such phenomena perturb the acousto-electric behavior of all the elements of the structure (the substrate and the layers which compose the wave-guide) and consequently yields a thermal drift of the synchronism frequency of the device. This frequency variation can be expressed as a relative function consisting in two major contributions in the vicinity of the reference temperature T0=25°C. The relative phase velocity variation which reads

dV V T = V  T -V T 0  V  T 0  (6)
The relative propagation distance variation is written as follows de e

T 

As the frequency is given as a first approximation by the relation f=V/2e for BAW (with e the actual thickness of the bulk resonator) and f=V/λ for SAW (with λ=2p, i.e. twice the mechanical period), the relative temperature sensitivity of the frequency reads

df f  T  = dV V  T  -de e  T  (8)
where λ replaces e fort SAW. For each considered temperature poin, the electromechanical problem is solved simply by changing the effective physical constants (elastic, piezoelectric and dielectric) and the mass density composing the structure. In this matter, one has to compute the thermal sensitivity coefficients of all the physical constants involved in the model in all the space directions as follows

c ijkl T =c ijkl T 0  ∑ n=0 N 1 c ijkl T 0  ∂  n  c ijkl ∂ T T -T 0  n = ∑ n=0 N κ ijkl  n ΔT n (9) L x i T = ∑ n=0 N ∂ n  L x i ∂ T T -T 0  n = ∑ n=0 N η x i  n ΔT n ( 10
)
where the sensitivity coefficients κ (n) ijkl and the thermal expansion coefficient η (n) are experimentally determined and accessible in the literature up to the third order at least for Quartz. For most materials however, only first or at best second order coefficients have been measured. Concerning the mass density, it is related to the thermal expansion of the elementary crystal lattice Vme (T) as follows:

ρ T = ρ T 0  V ol T  with V me T =L x 1 L x 2 L x3 (11)
However, in the case of FEA-BEM computation, the FEA part must be expanded to actually account for the period or the layer thickness changes corresponding to [START_REF] Slobodnik | Microwave Acoustics Handbook , vol 1A : Surface wave velocities[END_REF]. In the case of wave-guide composed of different materials (as for SAW with Aluminum electrodes atop thick single crystal substrates), a rigorous analysis of the problem needs a finite element analysis of the structure expansion accounting for differential thermal stresses. However, as in ref. [START_REF] Pastureaud | Prediction of the thermal sensitivity of surface acoustic waves excited under a periodic grating of electrodes[END_REF], we assume here that the substrate imposes its expansion to any film deposited atop its surface, but these materials do expand in the direction normal to the propagation surface along their own physical properties. This assumption is confirmed by experimental observations [START_REF] Lichtensteiger | Ferroelectric size effects[END_REF] on thin epitaxial PZT layers and was found efficient in [START_REF] Pastureaud | Prediction of the thermal sensitivity of surface acoustic waves excited under a periodic grating of electrodes[END_REF] to accurately simulate actual devices.

Following this process, the resonance frequency of any mode of the infinite grating or of the BAW resonator is computed for each temperature point to derive an analytical frequencytemperature law. We generally assume that this law can be established according to the numerous results obtained for BAW and SAW devices as a third order polynomial function which reads

Δf f = θ α  T -T 0   θ β  T -T 0 2  θ γ  T -T 0  3 ( 12
)
where θ α, θ β and θ γ respectively represent the first, second and third order TCF computed at room temperature T0 (25°C). This computation process only has a sense if it allows for accurately identifying the resonance frequency at each considered temperature. In that purpose, we use a parametric isomorphic function describing the admittance close to the resonance. The parameters of the model then are computed using a 3-point fit method described in [START_REF] Fusero | SSBWto-PSAW conversion in SAW devices using heavy mechanical loading[END_REF]. A simple approach for fitting the resonance frequency consists in fixing γ to a value corresponding to the actual excitation conditions applied to the grating (typically ½ for SAW) and then fit the admittance changes versus frequency only. The resonance then no more depend on the excitation parameter, yielding the following isomorphic fitting function

Y ω =Y r ω Y p ω ω p -ω ( 13
)
where ωp represents the resonance frequency, and Yp(ω) and Yr(ω) hold for the pole magnitude and the other contributions respectively. The 3-point fitting method quickly converges and reveals stable for well-posed problems. It can also be applied to identify anti-resonance parameters, simply changing the admittance to impedance in (13). This approach holds faor SAW and BAW as well. Finally it must be noted that all these developments are valid also for wave guide simulated by periodic FEA only, in and 3D as well. The method finally can be extended for any FEA structure submitted to quasi-static thermal perturbation provided the assumption of negligible thermo-differential effects remains valid.

III. VALIDATION RESULTS

A. BAW on Quartz

The very first case to treat in order to check the validity of such a calculation corresponds to BAW in quartz, for which second order compensation effects on AT cut is well known and predicted by simple analytical models. In that case, we consider a very simple mesh corresponding to a 2D plane surface resonator. The thickness was arbitrarily fixed to 50 µm to provide a resonance frequency near 35 MHz considering AT quartz. The width of the resonator also was arbitrarily fixed, but this figure had no impact on the results as the electrodes extend along the whole surfaces, Setting γ to an integer value (arbitrarily 1) yields the simulation of an infinitely long transducer, as it is generally the case for simple analytical models. For achieving the computation, Slobodnik et al. elastic, piezoelectric dielectric constants and density have been used [START_REF] Slobodnik | Microwave Acoustics Handbook , vol 1A : Surface wave velocities[END_REF] but the considered effective thermal coefficients of those physical constants were those of Bechmann et al. [START_REF] Bechmann | Higher Order Temperature Coefficients of the Elastic Stiffnesses and Compliances of Alpha-Quartz[END_REF]. Figure 2 shows the implemented mesh and fig. 3 a&b show the evolution of the resonance frequency for pure shear BAW on AT and BT quartz cuts respectively. The comparison with analytical model results reveal a good agreement between the two approaches, yielding a first validation element. Only FEA was used ion this first validation step.

B. Rayleigh wave on Quartz

The second unavoidable case to be addressed by the model corresponds to Rayleigh waves on (ST,X) Quartz and the sensitivity of the temperature compensation to electrode loading, as it also consists in a so-called "classics" which was treated along various approaches. In that case, we used the same data set as in section III.A, adding the contribution to Aluminum for which elastic and temperature coefficient was found in [START_REF] Landolt-Börnstein | Numerical data and functional relationships in science and technology, Group III, Crystal and solid state physics[END_REF]. In that case, we combined FEA and BEM, both being thermally perturbed as described above. Figure 4 shows the implemented mesh and fig. 5 shows the evolution of the Rayleigh wave frequency for 2 different metal thickness (100 and 200 nm), the wavelength being arbitrarily set to 10µm. It clearly shows that the turnover temperature shifts down with the electrode thickness, which actually corresponds to experimentally observed phenomena and agree well with theoretical results reported in [START_REF] Pastureaud | Prediction of the thermal sensitivity of surface acoustic waves excited under a periodic grating of electrodes[END_REF].

IV.THERMAL SENSITIVITY OF COMPLEX WAVE-GUIDE

Passivation layers are used to protect the surface of SAW devices from any dusts or surface pollution that could counteract the IDT operation. To illustrate the capability of our model to correctly predict the response of complex surface structures and their particular thermal sensitivity, we have theoretically and experimentally check the influence of an amorphous silica (SiO2) coating on a two port surface transverse wave (STW) resonator, yielding a so-called Love wave resonator. These devices consists in a four-port STW synchronous resonator operating close to 500 MHz (grating period p=5 µm, relative electrode height h/λ=1.5%, metallization ratio a/p=0.5) built on AT cut of quartz (propagation along Z') covered by a 0.5 µm thick amorphous silica layer. It presents 2 resonances corresponding respectively to the beginning and the end of the frequency stop-band. The stop-band width only can be accurately predict if considering a conformal deposition, i.e. the corrugation due to the metal strip is reproduced at the silica surface, as shown in fig. 6. A comparison between the harmonic admittance and the transfer function of the abovementioned resonator is reported in fig. 7, showing a good agreement between theory and experiment. We then have applied our TCF calculation process, allowing for identifying the behavior of each independently. The results are reported in fig. 8, showing that the two resonance does not exhibit exactly the same temperature behavior, although their polarization is about the same (phase quadrature). This actually corresponds to what is experimentally observed, even if the actual turnover temperatures occur at smaller temperature (cf Fig. 9). This of course can be controlled by changing the effective elastic temperature constants, as they are not known as accurately as those of quartz.

Figure 6. Mesh of the Love-wave grating, the metal ratio is set to 0.5, the mechanical period equals 5µm, the electrode thickness is set to 150 nm and the silica overlay is 0,5µm thick.

Figure 7. Superposition of harmonic admittances computed using the present approach and the above with the experimentally measured admittance of a passivated 4-port STW resonator

V. CONCLUSION

We have proposed an approach based on finite element analysis to simulate the thermal effects on the resonance frequency of complex wave guide. The basic computation has been validated for standard plane BAW resonators on Quartz.

The possibility to integrate radiation effects using boundary elements based on Green's functions also has been validated for Rayleigh wave on quartz too, as the most difficult problem generally consists in accurately predicting the temperature compensation effects. This work is applied also to address Love wave device thermal behavior for which no simple model is capable to account for the influence of the IDT covered by th silica overlay. This work will be extended in the near future to address the problem of FBAR and SMR also requiring an accurate computation of the unperturbed model to correctly simulate the resonance conditions.

Figure 1 .

 1 Figure 1. Typical 1-D periodic 2D wave-guide structure addressed by the model

Figure 2 .

 2 Figure 2. Mesh of the BAW resonator

Figure 4 .

 4 Figure 4. Mesh of the SAW grating, the metal ratio is set to 0.5, the mechanical period equals 5µm and the electrode thickness varies from 100 to 200 nm.

Figure 5 .

 5 Figure 5. Thermal sensitivity of SAW under the above grating (fig.4) for 2 different electrode thickness (100 and 200 nm).

Figure 8 .

 8 Figure 8. Thermal sensitivity of Love-wave under the grating of fig.6 for the two different modes emphasized in fig.7.

Figure 9 .

 9 Figure 9. Experimental temperature drift of the beginning and end-of-thestopband frequencies of passivated Love wave resonators.


I Boullet

J Marozeau

S Meunier
email: meunier]@lma.cnrs-mrs.fr

Temporal summation of loudness as a function of frequency and temporal pattern

This research is part of a study designed to propose a new model of loudness applicable to impulsive sounds. The purpose of the research was to study the effects of duration, frequency and temporal pattern of the signal on loudness. The aims were to determine: 1) critical duration as a function of frequency, 2) how does the loudness decrease below this critical duration, 3)and whether the results depend on temporal pattern (rectangular signal, signal with exponential decay). To measure loudness we selected the method of adjustment for its good precision/duration compromise. The stimuli used in this study were five pure tones, with durations of 5 ms to 1 s. We will present the values of critical durations that we obtained for rectangular and impulses signals. The time course of loudness change will be describe both for rectangular and exponential-decay patterns. The data of this study, added to results from the literature, are used to make a synthesis on loudness of impulsive sounds. Moreover, new results are given on more natural sounds that have mostly an exponential decay.

INTRODUCTION

The sensation of loudness of stationary pure tones can be easily predicted as a power function of intensity [START_REF] Stevens | On the psychophysical law[END_REF]. This relation is independent of the duration when the stimulus is longer than one second. Below this temporal threshold, the loudness decreases in a monotonic way with duration. Munson [Mun47] proposes as a model of this effect, a leaking integrator which stores at least momentary the energy. The relation between the intensity of an infinitely long and stationary sound (I ∞ ) and the intensity of a tone burst with a duration T which is perceived with the same loudness can be describe as follows:

I ∞ = I T (1 -e -t τ ), ( 1
)
where τ is a time constant. Many studies (see [START_REF] Scharf | Handbook of perception and Human Performance[END_REF] for a complete review) were dedicated to estimate this time constant, although results often differ (from 23 ms [H.59] to 120 ms [M73]), but the usual value seems to be around 100 ms. This time constant can also be used to estimate the threshold above which the duration does not affect anymore the loudness. Many different values have been proposed, but they also differ among studies (from 15 ms [START_REF] Dt | Loudness as a function of signal duration[END_REF] to 500 ms [START_REF] Berglund | Loudness as a function of the duration of auditory stimulation[END_REF]). The effect of level and frequency on this time constant has also been widely studied. Results again are not unanimous. If a good compromise seems to emerge concerning 1/5 the effect of the level with a time constant of 200 ms at hearing threshold and 100 ms above, no unanimous result has been found for the influence of frequency. The great variance of the results can be explained by the difference of nature of stimuli (their temporal shape, frequencies, level, spectral composition) and by the experimental protocol used to estimate loudness [FM96].

The present study is a part of a research project designed to elaborate a better predictor of the loudness of impulsive sounds. The temporal shapes of the stimuli used here are chosen to mimic as close as possible natural impulsive sounds, with a short attack time, no sustain and an exponential decay. To the authors' knowledge no or only few studies [ISA citation PLEASE] have used this kind of sound. Additionally, rectangular temporal patterns have also been used. The aims of this study is dedicated (1) to estimate the threshold above which the duration does not influence anymore the sensation of loudness of impulsive sounds, (2) to study its decrease below this threshold and (3) to evaluate whether the results depend on temporal pattern (rectangular signal, signal with exponential decay).

EXPERIMENT I

Stimuli

The stimuli were composed of pure tones of 125 Hz, 350 Hz, 1 kHz, 3 kHz and 8 kHz. Two types of temporal shapes have been used: a rectangular, and an impulse. Each type of shapes varies in duration. The rectangular shape is composed of an onset and offset made of a half Hanning window ramp of 5 ms each and a sustain part lasting 0, 12, 40, 90, 210 or 490 ms. The impulsive shape is composed of a linear slope that goes form 0 to a maximum value in 5 ms, followed by an exponential decay characterized by the duration needed for the amplitude to decrease from 90 % to 10 % of the max. The duration used were 2, 5, 10, 22, 50, 100, 220, 500 or 1000 ms. The number of stimuli were thus 75 (6 * 5 + 9 * 5). The maximum level of each stimulus was set at 72.5 dB SPL.

Methods

Loudness evaluations were collected using an adjustment method [START_REF] Gescheider | Psychophysics. Method and Theory[END_REF], in which the listener has to change the level of a comparison sound to set it to the same loudness of stimuli separately. Each stimulus was judged twice: once the comparison sound was presented before the stimulus, and once after. The comparison sound was a narrow band noise around 1kHz with a bandwidth of 120 Hz with rectangular temporal envelope of 500 ms. A silence of 400 ms was inserted between the two sounds. All stimuli were presented randomly. Ten subjects participated in this experiment, including 6 women and 4 men with normal hearing. This experiment takes place in an anechoic chamber. Stimuli were played via a Genelec 1031A loudspeaker.

Results

For each stimulus, averages were computed from the two judgments of the ten subjects. Figure 1 shows the measured loudness levels as a function of duration for rectangular (left panel) and impulsive (right panel) sounds. In order to keep a common definition of duration, the time when the amplitude is above -6dB of its maximum was used (as proposed by [FM96]).

These graphs show that the loudness increases as expected exponentially with the duration, but much more interesting is the great similarity between bursts of sounds and impulses. The temporal shape of the stimuli seems not to interact with the effect of duration on loudness. Unfortunately, because the durations used are not exactly the same, there is no way to evaluate statistically the similarity between the two sets of data. However it is possible to run an analysis of variance (ANOVA) for each type of shape separately. Two repeated-measured ANOVA were performed with the frequencies (5) and the durations (6 for burst, 9 for the impulsive) as main factors. For burst sounds, as expected the two main factors were significant at 5 % [F (4, 40) = 10.86; p < 0.001 for the frequency and F (5, 50) = 120.9; p < 0.001 for the duration]. But much more interesting no significant effects were found for the interaction between the two factors [F (20, 200) = 1.14; p = 0.3]. This result suggests that no effects of frequency on influence of duration on loudness were found. All the lines in the left panel of figure 1 can be so considered as parallels. For impulsive sounds, the two main factors were also significant [F (4, 40) = 7.47; p < 0.001 for the frequency and F (8, 80) = 36.43; p < 0.001 for the duration]. But the interaction was also significant [F (320, 320) < 3.04, p < 0.001]. This result can be ascribable mainly to the short 125-Hz stimuli. Figure 1 shows explicitly that the loudness of these stimuli were judged much louder than excepted for short durations. This can be explained by an error in the conception of these sounds, since a short temporal impulsive envelope cannot properly modulate a low frequency signal. When this frequency is discarded from the analysis, the interaction becomes also not significant [F (24, 240) < 1].

EXPERIMENT II

Because no explicit temporal threshold was found, a new experiment was run with stimuli with a wider duration range. Also because no effect of the frequency has been found, only one frequency at 1-kHz was used. The two shapes were kept. The sustain part of the rectangular sounds lasted either 10, 22, 50, 100, 220, 360, 500, 1000 or 1500 ms. Durations needed for the impulses amplitude to decrease from 90 % to 10 % of the max were either 2, 5, 10, 22, 50, 100, 220, 360, 500, 1000, 2000, 4000, or 8000 ms. Figure 2 shows the loudness evaluation as function of the log of the duration (defined again as length were the stimulus were above -6dB of its maximum values) for the rectangular (left panel) and the impulsive (right panel) sounds.

In order to describe the effect of duration on loudness, the equation 1 which describe the relation intensity-duration was transformed to a loudness-duration relation as follows: where L comp is the level of the comparison sound, and L T ∞ the level of the tested sound at a duration above the temporal threshold. The time constant for rectangular sounds was 100 ms, which corresponds perfectly with the literature. Time constant of 77 ms for impulses is more delicate to comment because the figure shows explicitly that the model does not fit well the data. We can observe that the loudness still increases in a monotone way with duration but the leaking integrator does not model well this relation. This can be explained partly by the apparition of another mechanism used by subjects to judge the global loudness of a time varying sound [START_REF] Meunier | Loudness of sounds with temporal variable intensity[END_REF].

L comp = L T ∞ + 10 * log10(1 -e -t τ ) (2) 3/5

CONCLUSIONS

These experiments were dedicated to study the effect of duration on the loudness of impulses sound. The first experiment reveals that for short sound (under 500 ms) no effect of the frequency was found. The loudness of impulses sounds seems to be affected by the duration approximately in the same way that the rectangular sounds. When considering a wider scale of duration, differences between bursts of sounds and impulses appear. The first one can be described efficiently with a leaking integrator with a time constant of 100 ms, and a threshold above which the duration does not affect anymore the loudness can be determined approximately around 500 ms. For impulses sounds, no threshold was found. Loudness still increases with duration above 1 second and more data have to be collected to model this phenomenon.

Figure 1 :

 1 Figure 1: Results of the experience I, for the rectangular on the left panel and the impulsive sounds on the right

Figure 2 :

 2 Figure 2: Results of the experience II, for the rectangular on the left panel and the impulsive sounds on the right

ACKNOWLEDGEMENTS

We would like to thank Georges Canevet for his useful help and advices.


Quoc Thong Nguyen
email: quoc-thong.nguyen@univ-ubs.fr

Yves Delignon

Anh Thu Phan

Probabilistic Modelling of Printed Dots at the Microscopic Scale

Keywords: Probabilistic model, Bernoulli process, Metropolis Hastings within Gibbs, Microscopic printing, Markov chain

Microscopic analysis of paper printing shows regularly spaced dots whose random shape depends on the printing technology, the configuration of the printer as well as the paper properties. The modelling and identification of paper and ink interactions are required for qualifying the printing quality, for controlling the printing process and for application in authentication as well. This paper proposes an approach to identify the authentic printer source using micro-tags consisting of microscopic printed dots embedded in the documents. These random shape features are modelled and extracted as a signature for a particular printer. In the paper, we propose a probabilistic model consisting of vector parameters using a spatial interaction binary model with inhomogeneous Markov chain. These parameters determine the location and describe the diverse micro random structures of microscopic printed dots. A Markov chain Monte Carlo (MCMC) algorithm is thus developed to approximate the Minimum Mean Squared Error estimator. The performance is assessed through numerical simulations. The real printed dots from the common printing technologies (conventional offset, waterless offset, inkjet, laser) are used to assess the effectiveness of the model.

Introduction

Preventing and discouraging unauthorised printed materials is the expected achievement of the authentication methods. Papers dealing with this subject can be categorised in two strategies. The first one consists in embedding an extrinsic signature such as secure tags in the prints [START_REF] Picard | Towards fraud-proof id documents using multiple data hiding technologies and biometrics[END_REF][START_REF] Villán | Multilevel 2-d bar codes: toward high-capacity storage modules for multimedia security and management[END_REF]. While the second one, which is the motivation of this paper, consists in characterising the intrinsic features of printing and scanning process [START_REF] Chiang | Printer and scanner forensics[END_REF]. As examples of intrinsic features, the signature consisting on the banding artefact coming from fluctuation of the optical photoconductor angular velocity in case of laser printers [START_REF] Ali | Intrinsic and extrinsic signatures for information hiding and secure printing with electrophotographic devices[END_REF], dimple effect which is specific to inkjet printers or also texture features [START_REF] Ryu | Electrophotographic printer identification by halftone texture analysis[END_REF], or the analysis of quality signature of the unique print to differentiate one printer technology/supplier from another [START_REF] Oliver | Use of signature analysis to discriminate digital printing technologies[END_REF]. Authentication can be performed based on paper statistics, for example, the paper statistics is extracted by capturing microstructure images from different carton packages under different light conditions and with different cameras. These microstructure images are considered as digital fingerprinting on which the authentication frameworks are derived [START_REF] Voloshynovskiy | Towards reproducible results in authentication based on physical non-cloneable functions: The forensic authentication microstructure optical set (famos)[END_REF]. In [START_REF] Sharma | Paperspeckle: microscopic fingerprinting of paper[END_REF], the authors proposed the intrinsic feature of paper based on a texture speckle pattern, a random bright/dark region formation at the microscopic level when light falls on to the paper. The extracted signature of printer is a unique pattern of misplaced toner powder on each paper [START_REF] Zhu | Print signatures for document authentication[END_REF]. For electrophotographic printers, the intrinsic signals extracted from banding artefact are used to identify the device [START_REF] Mikkilineni | Signature-embedding in printed documents for security and forensic applications[END_REF]. Printer source identification is an important forensic tool.

The techniques using characters for identification of source scanner and printer have been proposed [START_REF] Khanna | Survey of Scanner and Printer Forensics at Purdue University BT -Computational Forensics: Second International Workshop[END_REF]. Source identification using printed character in different languages has been recently studied [START_REF] Tsai | Digital forensics of printed source identification for chinese characters[END_REF][START_REF] Tsai | Japanese character based printed source identification[END_REF][START_REF] Tsai | Digital forensics for printed character source identification[END_REF][START_REF] Tsai | Source identification for printed arabic characters[END_REF][START_REF] Tsai | Digital forensics of microscopic images for printed source identification[END_REF].

Besides the role in print authentication, depending on the purposes of the problems, there are a number of applications concerning the printing models. A specific model is considered in the scale of the experiments/applications. For instance, in the human vision scale, halftoning technique is a method that takes advantage of the optical illusion of human vision to display continuous tone grayscale image with only black or/and white dots [START_REF] Ulichney | Review of halftoning techniques[END_REF][START_REF] Velho | Digital halftoning with space filling curves[END_REF][START_REF] Pappas | Model-based halftoning[END_REF]. The scale of the model is the resolution of the printer. An image of print not only contains the properties of a printing process but also includes image degradation; Baird proposed a physics-based model in which the parameters defining imaging defects are described [START_REF] Baird | Document image defect models[END_REF]. Under the microscopic scale printing, a probabilistic model for the average coverage of the toner particles ink of electrophotographic printing process is introduced by Norris et al. [START_REF] Norris | Printer modeling for document imaging[END_REF]. However, the model has a limited validation, and the estimation procedure has not been considered. For authentication purpose, we propose a probabilistic modelling of microscopic printed dots.

The shapes of microscopic printed dots can be considered as the intrinsic feature of the printing process. At the microscopic scale, a dot is a random pattern whose shape depends on the technology, the setting of the printer, the ink quality and/or the paper properties. From a statistical point of view, the digital acquisition of these random dots can be modelled as a spatial interaction binary model based on an exponential power kernel that depends on location parameter and shape parameters. The accurate and consistent model can help to create an authentic signature to improve the authentication procedure. The model proposed in our previous works [START_REF] Nguyen | Printer technology authentication from micrometric scan of a single printed dot[END_REF][START_REF] Nguyen | Printer identification from micro-metric scale printing[END_REF] consists in a spatial distribution model that simulates the randomness of the shape of printed dots. The random shape is manipulated by the parameters in this probabilistic model. In [START_REF] Nguyen | Modélisation de points imprimés à l'échelle micro-métrique[END_REF], the numerical estimation procedure using MCMC algorithm is developed. However, in order to cover the diverse random shape of printed dots, the surrounding pixels or the conditional probabilities has to be taken into account.

In this paper, we propose a spatial interaction binary model that describes more precisely the random shape of the microscopic dots by taking into account the neighbourhood of the pixels. In particular, the Markov property between pixels in the image is considered. This property helps to describe more accurately the intrinsic features of the print, which can thus improve the discrimination power. An MCMC algorithm is developed for estimating the model parameters. Finally, we analyse the performance of the proposed unsupervised identification algorithm. The rest of this paper is organised as follows. Section 2 is devoted to visually describe the randomness of shape of microscopic dots. The formulation of spatial stochastic model is discussed in Section 3. The Metropolis-Hastings within Gibbs algorithm for the estimation is developed in Section 4. In Section 5, the Posterior Cramér -Rao bound used to evaluate the performance of the algorithm is derived. Section 6 is devoted to the performance of the new algorithm which are analysed through simulations, and the benefit of such a model and estimation algorithm are illustrated in the case of real microscopic scale printed dots. The basic authentication of printer identification using the model is also described in Section 7.

Concluding remarks are given in Section 8.

Printed dots at the microscopic scale

As aforementioned, we consider the shape of microscopic printed dots as an intrinsic feature to perform authentication on printer identification. This section is devoted to describing the random shape of microscopic printed dots. We then demonstrate intuitively that the shapes are different between different printers, which can benefit printers discrimination. Let us start by the observation of printed dots at the microscopic scale. A single printed dot is created by printing an image containing one black pixel. In our study, the microscopic printed dots are captured using an optical Zeiss Microscope with an AxioCam camera. The camera has a maximum solution of 2464 x 2056 (5Mp) with a sensor size of 8.5 mm x 7.1 mm and a pixel size of 3.45 µm x 3.45 µm. Since the focus of the study is on the physical shape of coverage area of the ink, the collected images are only in black and white.

The images on the left column in Fig. 1 illustrate the dots printed by offset, laser, inkjet printers respectively. Random shapes of these dots are observed under microscope scale. The physical size of the dots is very small. For example, size of a dot from the offset printer (1200 dpi) is theoretically 21.16µm. The printing resolution is an important element in this study. If the resolution is too high, the outcome will be similar to the dust noise from the printer. For instance in our study, the dots printed from a laser printer under resolution 1200 dpi are not distinguishable from the noise.

The resolution can neither be too low, since the shape of the dots are no more random except for the edge. Therefore, the resolution has to be chosen manually for a specific printer.

Visually, from the left column of Fig. 1, we can see that the shapes of printed dots from different printers are dissimilar. two different laser printers at the same resolution 600 dpi. On the right column in Fig. 1, each image illustrates the frequency of appearance of black pixels at each position (profile) from 100 captured samples. As observed, the probabilities of emergence of the black pixels are lower when the pixels are further from the centres. The shape of the profile contains the characteristics of the printed dots, which means the shape of the dots contains the identity of the source printer in some sense. The model in [START_REF] Nguyen | Printer technology authentication from micrometric scan of a single printed dot[END_REF][START_REF] Nguyen | Printer identification from micro-metric scale printing[END_REF] introduced the parameters controlling the width of the droplets, the density of the black pixels and the average number of black pixels. However, the interaction with the surrounding pixels was not taken into account. In particular, when a pixel is black, its neighbour is more likely to be blackened. In addition, when the pixels 100 are further from the centre of the dots, this influence of the black pixel(s) on their vicinity is also smaller. Including these properties could help the model to have a more precise simulation, thus describing in more details the identity of source printer. These properties are detailed in the next Section.

Spatial interaction binary model with Markov chain 105

Let S be a finite set corresponding to N pixels of a binary image U that displays K dots. We consider the random process U = (U s ) s∈S , where each U s is a Bernoulli random variable which takes its value from {0, 1}, respectively black state and white state.

Let V s be a neighbourhood of the pixel s, the geometric shape of V s is independent of s ∈ S. U is modelled by a Markov random field

P (U s = u s | (U t = u t ) t =s ) = P (U s = u s | (U t = u t ) t∈Vs ) , (1)
i.e. the probability that the pixel s takes the value u s conditional on all other pixels of the image equals to the probability of u s conditional on the values of the pixels in the neighbourhood V s . This probability also has to depend on the distance of the pixel s to the centre µ, which means P (U s = 0 | (U t = u t ) t =s ) decreases when the pixel s goes further from the centre, i.e. the Euclidean distance ||s -µ|| increases.

110

A two-dimensional image can be transformed into a one-dimensional chain by

Hilbert curve [START_REF] Hilbert | Ueber die stetige abbildung einer line auf ein flächenstück[END_REF] (see Fig. 3). Because the targeted interactions (statistical dependencies) between pixels within the microscopic dot are well preserved in the spatial structure by Hilbert curve, the Markov chain with Hilbert path is used in many image processing applications [START_REF] S.-I. Kamata | Interactive analysis of multi-spectral images using a hilbert curve[END_REF][START_REF] Fjortoft | Unsupervised classification of radar images using hidden markov chains and hidden markov random fields[END_REF][START_REF] Dridi | Em-based joint symbol and blur estimation for 2d barcode[END_REF]. Moreover, Pseudo-Hilbert Scan algorithm can be applied to arbitrarily-sized image [START_REF] Zhang | A pseudo-hilbert scan algorithm for arbitrarily-sized rectangle region[END_REF]. In convenient index, let U = (U 1 , U 2 , ..., U N ) be a vector of random variables ordered according to the Hilbert curve transformation of the image. In case the neighbourhood V s is restricted to one pixel, the image U becomes a Markov chain The Markov chain (2) encompasses the parameters of the initial probability

P (U = u) = P (U 1 = u 1 ) N i=2 P (U i = u i | U i-1 = u i-1 ). (2)
P (U 1 = 0), (3)
and transition probabilities

P (U i = 0 | U i-1 = u), (4)
for 1 < i ≤ N . The parametric probabilities for (3), (4) are

P (U 1 = 0) = 1 - K k=1 (1 -p k (s 1 )) (5)
P (U i = 0 | U i-1 = u) = 1 - K k=1 (1 -p k (s i ))λ (1-u) . (6)
The function p k (.), named kernel, measures the ability to blacken the surrounding pixels of the k-th dot. The kernel p k (s) has to decrease with respect to the distance ||s -µ k ||, so that, as requisite property, the initial and transition probabilities also decrease with respect to the distance. The support of p k (.) is a subset of (0,1].

The positive parameter λ characterises the interaction between pixels when the neighbour pixel is black, U i-1 = 0. With λ = 1, the probabilities (6) only depend on the relative position to the centres, P (U i | U i-1 ) = P (U i ), which is the approach proposed in [START_REF] Nguyen | Printer technology authentication from micrometric scan of a single printed dot[END_REF][START_REF] Nguyen | Printer identification from micro-metric scale printing[END_REF]. In case λ < 1, the smaller λ leads to a higher chance of being black of current pixel U i , since P (U i = 0 | U i-1 = 0) is greater with smaller λ (see Eq. ( 6)). When the neighbour pixel is white, U i-1 = 1, the probabilities depend on the distances of the current pixel to the dots centres. With these properties, the dots are more compact, and the density of black pixels at the areas closer to the centres is higher. The new proposed parameter λ helps to describe a wider variety of shape of dots.

The transition probabilities depend on p k (s) which relates to the position of the pixels, this implies that U is a non-stationary Markov chain. The expected number of black pixels of the image, denoted n b , is calculated as

n b = N i=1 E (δ(U i )) = N i=1 P (U i = 0), (7)
with δ(x) = 1 if x = 0, δ(x) = 0 elsewhere. The probability P (U 1 = 0) is computed by [START_REF] Ryu | Electrophotographic printer identification by halftone texture analysis[END_REF], while the other terms are calculated by [START_REF] Oliver | Use of signature analysis to discriminate digital printing technologies[END_REF] with the recurrence relation and the law of total probability

P (U i = 0) = P (U i = 0 | U i-1 = 0)P (U i-1 = 0) +P (U i = 0 | U i-1 = 1)P (U i-1 = 1), (8)
with 2 ≤ i ≤ N . Let us note that with smaller λ, P

(U i = 0 | U i-1 = 0) is greater,
that leads to the value of P (U i = 0) also becoming greater. This implies that smaller λ statistically increases the number of black pixels. In the following, the kernel p k is chosen as a Gaussian power kernel

p k (s) = η exp    - s -µ k 2 σ 2 β    . (9)
Fig. 4 illustrates the simulation of one dot with σ 2 , β equal to 70 and 5 respectively 130 and various values of λ and η. It can be seen that the dots with small λ are more continuous, in other words, the black pixels make their neighbourhood more likely to be black. The parameter η controls the amplitude of the kernel and the density of black pixels increases with η. When η is large, the density of black pixels near the centre is almost independent of λ.

θ(u) = θp(θ | u)dθ, (10)
where p(θ | u) is the posterior distribution. To approximate the integral [START_REF] Mikkilineni | Signature-embedding in printed documents for security and forensic applications[END_REF], an MCMC algorithm is applied [START_REF] Robert | Monte Carlo Statistical Methods[END_REF]. Since the vector of parameters is of high-dimension, it is more efficient to use Gibbs sampling procedure which consists in sampling successively each parameter from their conditional posterior distribution. Unfortunately, since no closed form expression is available for these distributions, we propose to use a Metropolis-Hastings within Gibbs algorithm [START_REF] Robert | Monte Carlo Statistical Methods[END_REF]. The chosen strategy is to run a random walk as proposal distribution:

θ * = θ (t) + t , ( 11
)
where t is a zero mean normal random vector with a diagonal covariance matrix.

We generate θ 2 = log(σ 2 ), θ 3 = log(β) instead of σ 2 , β to ensure their positivity; a transformation is also applied for η and λ to generate θ 4 , θ 5 respectively, i.e. Algorithm 1 Sampling algorithm with random walk

1: Initialise θ (0) = θ (0) 1 , θ (0)
2 , θ Generate θ b ∼ N (θ

(0) 3 , θ (0) 4 , θ
(i-1) b , σ 2 b ) 5:
Compute the acceptance ratio α(θ

(i-1) b , θ b ) = min    p θ b | θ (i-1) -b , u p θ (i-1) b | θ (i-1) -b , u , 1    6:
Decide

θ (i) b = θ b with probability α θ (i-1) b , θ b θ (i-1) b otherwise 7:
end for 8: end for Where θ

(i-1) -b = θ (i) 1 , . . . , θ (i)
b-1 , θ

(i-1) b+1 . . . , θ (i-1) 5
.

The initialisation of θ proceeds as following. The centres are initialised by the careful seeding [START_REF] Arthur | k-means++: The advantages of careful seeding[END_REF]. The other parameters are respectively initialised as σ 2 = 0.6 × N πK , β = 2. Since η and λ are in (0, 1], intuitively, we may initialise both of them equal to 0.5. From Bayes' law, we have:

p (θ b | θ -b , u) = p (u | θ b , θ -b ) p (θ b , θ -b ) p (θ -b , u) , (12)
where p (u | θ b , θ -b ) is the likelihood and p (θ b , θ -b ) is the prior. Since all parameters are independent, the acceptance ratio in Alg. 1 is rewritten as

α(θ (i-1) b , θ b ) = min    p (u | θ b , θ -b ) p b (θ b ) p u | θ (i-1) b , θ -b p b θ (i-1) b , 1    , (13)
with p b the prior of parameter θ b . The prior distributions are chosen with fixed pa-rameters by experience. Each centre {µ k } has multivariate normal distribution, σ 2 is log-normal distributed, β follows gamma law, while both η and λ are beta-distributed.

Let us note that when the priors are closed to uniform, the acceptance ratio mainly depends on the likelihood ratio. The choice of the variances of the priors depends on the confidence of the initialisation and the experience.

Posterior Cramér -Rao bound

Let u be the observation vector and θ be the r-dimensional vector of random parameters (with r = 2k + 4). The joint probability density of the pair (U, θ) is p(U, θ).

Let θ(u) denote an estimate of θ. The Posterior Cramér-Rao Bound (PCRB) states that [START_REF] Tichavskỳ | Posterior cramér-rao bounds for discrete-time nonlinear filtering[END_REF][START_REF] Van Trees | Detection, estimation, and modulation theory[END_REF]:

E[( θ(U) -θ)( θ(U) -θ) T ] ≥ I -1 , (14)
where I is r × r Fisher information matrix (FIM)

I = E ∇ θ log p(U, θ)∇ T θ log p(U, θ) = -E ∆ θ θ log p(U, θ) , (15)
where ∆ θ θ := ∇ θ ∇ T θ is the second derivative operator and ∇ θ is the gradient operator with respect to θ.

Using the fact that p(U, θ) = p(U | θ)p(θ), the expression of the FIM in (15) can be expressed as:

I = -E ∆ θ θ log p(U | θ) -E ∆ θ θ log p(θ) = I d + I p , (16)
where I p stands for the a priori information matrix and

I d is the "standard" FIM (I d (θ)
is used in the derivation of the CRB) averaged over the a priori distribution of θ:

I d = Θ I d (θ)p(θ)dθ. (17)
In many special cases, when only the i-th component of θ, say θ i , is of interest, the PCRB is expressed as

E[( θi (U) -θ i ) 2 ] ≥ [I -1 ] ii (18)

Experimental results and performances

In this section, we aim at assessing both the performance of the estimators and the adequacy of the model. The Metropolis-Hasting within Gibbs algorithm is analysed, the performance of the algorithm is evaluated from the images generated by the parametric model. The experiment on the real captured dots is also conducted.

Synthetic Data

The simulated images are generated with a size of 64 × 64. As mentioned previously, the priors are fixed, the expected values of the prior are set by the initial values, and the variance are chosen large in order to consider a highly uncertain prior knowledge. The estimation is implemented with 20000 MCMC iterations and a burnin period of 5000, which means that 15000 last samples are used to approximate the parameters over 20000 generated samples.

Analysis of the estimation with respect to β

In this analysis, only β is varying. The simulated images of single dot are generated. The true centres µ of the dot is [START_REF] Tichavskỳ | Posterior cramér-rao bounds for discrete-time nonlinear filtering[END_REF][START_REF] Tichavskỳ | Posterior cramér-rao bounds for discrete-time nonlinear filtering[END_REF]. With σ 2 = 50, η = 0.8, λ = 0.6, along with three values of β, β ∈ {0.8, 2, 5}, Tab. 1 shows the mean squared errors of the estimators. Larger value of β leads to smaller errors for σ2 and η. In contrast, the errors of β and λ are larger with large β. With a large value of β, the simulated dots are different mainly near the edges. It means that the information of λ is less significant, which makes the estimation of λ less consistent. While with smaller β, the contribution of λ to the connected black pixels is more visible (see Fig. 5). On the other hand, the simulated dots look similar for large values of β (see Fig. 6), which also makes the estimation of β less consistent. The mean squared errors of the estimators with various values of η, σ 2 = 50, β = 2, λ = 0.6.

Analysis of the estimation with respect to λ

The simulations of one dot were conducted with β = 2, σ 2 = 50, η = 0.8, along with three values of λ, λ ∈ {0.2, 0.6, 0.9}. Tab. 3 shows the MSE of the estimators.

The estimators of σ 2 , β, η have decreasing MSE when λ increases. Since the number of black pixels decreases when the value of λ increases, the information for the estimation of σ 2 , β, η exposes clearer on the realisations. On the other hand, the estimation of λ is better for small values, because the connection between black pixels strongly emerges.

σ 2 β η λ 0.2 66.06 0.247 0.026 0.0017 0.8 45.12 0.23 0.011 0.0068 0.9 36.57 0.199 0.0113 0.0181 Table 3: The mean squared errors of the estimators with various values of λ, β = 2, σ 2 = 50, η = 0.8.

Posterior Cramér-Rao bound for single dot case

The PCRB is used to assess the precision of the estimators. Since the posterior information matrix I in ( 16) cannot be calculated analytically, it is approximated by Monte Carlo simulation. The random values of vector of parameters θ generated N times by the priors. For each generated vector of parameter, M images are simulated to approximate the standard Fisher information matrix I d (θ) in [START_REF] Velho | Digital halftoning with space filling curves[END_REF]. In Tab. 4, we demonstrate the approximated values of the Posterior Cramér-Rao Bound in case of single dot with N = 100 generated random values of vector θ. The values of priors are as follows: centre (µ x , µ y ) has multivariate normal distribution with mean [START_REF] Tichavskỳ | Posterior cramér-rao bounds for discrete-time nonlinear filtering[END_REF][START_REF] Tichavskỳ | Posterior cramér-rao bounds for discrete-time nonlinear filtering[END_REF] and covariance matrix diag(5, 5); σ 2 is log-normal distributed with mean 50 and variance 100; β follows gamma law with mean 2 and variance 1.5; η and λ are beta-distributed B(3, 1.5). For each value of θ, M = 40 images are generated to approximate matrix I d (θ). Specifically, the size of each image is 64 × 64. With each generated image, the parameters are estimated using Algo. 1. In total, 40 × 100 images are used to approximate the Mean Squared Error. As we observe in this example, σ 2 is estimated least accurate.

The performance with many dots

The impact of the number of dots, K, is analysed. Instead of one dot on each image as the analyses above, each image contains multiple simulated dots. The experiment is set up with σ 2 = 50, β = 2, η = 0.8, λ = 0.6, and all generated dots are completely separated from each other. As demonstrated in Tab. 5, the estimators are clearly more consistent for the images with more dots. The accuracy and consistency of β, η, λ is crucial in authentication. Therefore, analysing many dots is a good mean to have a more consistent estimation of the parameters. Since more dots give more information about the shape, which are determined by the parameters of the model, the estimators are more consistent. The standard deviation of the estimators with different number of cluster-dots K, σ 2 = 50, β = 2, η = 0.8, λ = 0.6.

Estimation results from printed dots

The dots from four printers of four popular technologies are taken: conventional offset, waterless offset, laser and inkjet printings. The resolutions of each printing process for each sample are respectively 1200 dpi with offset printers, 600 dpi with laser printer, and 720 dpi with inkjet one. For each printer, 10 images of size 64x64 of one dot are collected for the experiment. Fig. 7 represents the images of the real dots and the realisation generated from the estimated parameters obtained from the developed MCMC method. Metropolis-Hasting within Gibbs algorithm is conducted with 30000 MCMC iterations with 10000 as burn-in period. The realisation from the Spatial binary model [START_REF] Nguyen | Printer identification from micro-metric scale printing[END_REF] is illustrated to compare the improvement as well. Tab. 6 gives the numerical values of the estimators by MCMC estimation method as well as the standard deviation values. Fig. 7 visually compares the printed dots, the realisations from the proposed Markov spatial model and from the spatial binary model that does not take into account the neighbour dependence developed in [START_REF] Nguyen | Printer technology authentication from micrometric scan of a single printed dot[END_REF][START_REF] Nguyen | Printer identification from micro-metric scale printing[END_REF]. The results are very promising when the interaction is taken into account. The realisation from Markov spatial model is more realistic to the real dot than the one from spatial binary model. With the offset and inkjet images, the fragment near the border is resembled better, which makes them more similar to the real dots. The shape of the dot from inkjet printer is more complicated to simulate due to the concave parts.

With the print from laser printer, we can see the substantial improvement using the Markov spatial model w.r.t. the spatial binary model. The simulated dot is much more realistic and continuous.

In Fig. 8 and Tab. 7, we give other estimation of four printed dots with the resolution 600 dpi from HP-600 M620 laser printer. We can again have a visual comparison to show the improved simulation using the Markov spatial model. The capability of the proposed model to capture more details can clearly increase the power of discrimination between printers. Moreover, since the characteristics of printed dots are represented by the parameters, the proposed model may significantly reduce training data than using machine learning generative models to characterise the printed dots.

Application in printer source identification

The capability of simulating the random shape of printed dots of the model and the performance of the estimation are demonstrated in previous Section. In this part, we use the typical random shape which is characterised by the parameters as a "signature" to identify the printer source. The vector of three parameters (β, η, λ) can be used as a characteristic of a printer with a specific setting. These three parameters are dimensionless, which means that the size of image of dots does not affect the value of these parameters. With the model, we are able to set up an authentic tag from an exact printer by setting specific values for the parameters of the printer (ink toner, paper...). Under a specific setting, a particular printer statistically creates the same shape of microscopic dots. In this experiment, we applied the maximum likelihood classification method to discriminate the printed dots from different sources [START_REF] Foody | Derivation and applications of probabilistic measures of class membership from the maximum-likelihood classification[END_REF][START_REF] Asmala | Analysis of maximum likelihood classification on multispectral data[END_REF]. Fig. 9 summarises the procedure of the proposed printer identification method.

The estimations of the shape parameters are displayed in Fig. 10 that shows respectively β vs. λ, β vs. η and λ vs. η. Each printer is 3D Gaussian cluster in the space of parameters. Their location and shape are respectively given by their means (Tab. 6)

and their covariance matrices (Tab. 8). For each printer (conventional offset, waterless offset, laser and inkjet), we take 10 other images for the discrimination experiment. In total, 20 images for each printer with 10 images are used to approximate the parameters of 3D Gaussian cluster. In this specific experiment, all the dots are printed on uncoated paper, and other conditions are identical during the experiments. In this case, the normal distributions of the shape parameters enable to use the maximum likelihood classification method. All the samples are correctly classified (100%) with respect to their printer sources. The confusion matrix is presented in Tab. 9 Other experiment is also conducted on two laser printers, a HP laser printer and a Dell laser printer. The resolution of 600 dpi is set up for both printers, and the samples are printed on uncoated paper. The relation between estimated parameters of each printer is given in Fig. 11 which clearly separates two printers. Their location and shape of the 3D Gaussian clusters of each printer are respectively given by the means in Tab. 10 and the covariance matrices in Tab. 11. Similar to the previous experiment, 20 images are captured from each printer with 10 images are used to approximate the parameters of 3D Gaussian cluster. The maximum likelihood classification method also gives a completely accurate classified result (Tab. 12). These results point out the relevance of the model for printed dots at the microscopic scale.

As demonstrated in Fig. 8, Tab. 7 and the analysis in the simulation experiment section, the estimation of parameters is more consistent with many dots. Therefore, the analysis of many dots will help to improve the accuracy of discrimination an exact printer source of a document from the other printers. In practice, a micro-tag consisting in printed dots can be embedded in an authentic document.

Concluding remarks

In this paper, we have proposed a novel statistical model that considers the spatial interaction of printed dots at microscopic scale. This model allows to describe more accurately the formation of the shape of the binary dot. The spreading of the ink on the surface was described as the emergence of the black pixels and the interaction between pixels that depends on the distance between the centres of the dots. All these features are modelled by non-stationary Markov chain.

An estimation method for the parameters based on the Metropolis-Hasting within Gibbs algorithm is proposed. The performance of this approach is assessed, not only

Fig. 2 Figure 1 :

 21 Figure 1: Binary image of 4 single dots printed by an offset printer (a), the empirical probabilities of appearance of black pixels in the given position (profile) of 100 samples of one printed dot with the resolution 1200 dpi (b). The dots from a laser printer with the resolution 600 dpi (c), and the profile (d). The dots from an inkjet printer at 720 dpi (c), and the profile (f).

Figure 2 :

 2 Figure 2: The dots printed by two laser printers, Dell laser printer (a) and HP laser printer (b). The conditions are identical, resolution is 600 dpi.

Figure 3 :

 3 Figure 3: Scanning of 8 × 8 image

A

  Bayesian framework is employed to estimate the model parameters. Let us recall our problem, the parameters of interest in (5), (6) are θ U = {µ k } , σ 2 , β, η, λ which belong to a 2K + 4 dimensional space. The Minimum Mean Squared Error (MMSE) estimator θ of θ U given the observed data u is

Figure 4 :

 4 Figure 4: The simulation of the model (5), (6) with Gaussian power kernel with σ 2 = 70, β = 5 and different values of λ, η.

η = 1 1+e θ 4 , λ = 1 1+e θ 5 .

 45 For convenient notation, we rewrite the vector of interested 140 parameters as θ = [θ 1 , θ 2 , θ 3 , θ 4 , θ 5 ]. The Metropolis-Hastings within Gibbs using random walk is summarised in Algo. 1.

2 : 3 :

 23 for i = 1 to N iter do for b = 1 to 5 do 4:

(a) λ = 1 ( 2 Figure 5 :

 125 Figure 5: The role of λ with small β = 0.8.

(a) β = 7

 7 (b) β = 7.5 (c) β = 8

Figure 6 :

 6 Figure 6: Realization with three different values of β.

Figure 7 :

 7 Figure 7: The estimation from four printing processes and the realisations from the models.

Figure 8 :Figure 11 :

 811 Figure 8: The real dots vs. the regenerated image, the approximated centres are demonstrated.

Table 1 :

 1 The mean squared errors of the estimators with various values of β, σ 2 = 50, η = 0.8, λ = 0.6.

	σ 2	β	η	λ
	0.5 100.9 0.2041 0.0065 0.0058
	0.8 45.12	0.23	0.011 0.0068
	0.97 49.85 0.1845 0.032	0.01

Table 2 :

 2

Table 4 :

 4 The approximated PCRB and the Mean Squared Error (MSE) of each parameter in case of single dot.

Table 5 :

 5

		1 dot 2 dots 4 dots 7 dots
	σ 2 6.03	3.08	2.42	1.92
	β	0.25	0.19	0.18	0.15
	η	0.08	0.04	0.04	0.03
	λ	0.08	0.06	0.05	0.03

Table 6 :

 6 The estimators of the real dots by MCMC.

	Conventional Offset Waterless offset
	Avg.	Std.	Avg.	Std.
	β 3.71	1.4	6.55	1.74
	η 0.96	0.05	0.99	0.0004
	λ 0.67	0.15	0.52	0.13
		Laser	Inkjet
		Avg. Std. Avg. Std.
	β 1.21 0.92 2.28 0.51
	η 0.59 0.32 0.99 0.01
	λ 0.30 0.06 0.20 0.04

Table 7 :

 7 The estimators of σ 2 , β, η, and λ from four printed dots.

Table 8 :

 8 The covariance matrices calculated from the estimators of four printers.

	Actual/Predicted Con. offset Wat. offset Laser Inkjet
	Con. offset	20	0	0	0
	Wat. offset	0	20	0	0
	Laser	0	0	20	0
	Inkjet	0	0	0	20

Table 9 :

 9 Confusion matrix for ML classification of four different types of printer: conventional offset, waterless offset, laser and inkjet.

Table 10 :

 10 The mean values of the estimators of HP laser printer and Dell laser printer.

	0.0147	0.001 -0.0003	0.243 0.054 -0.003
	0.001	0.009	0.0009	0.054 0.028 0.001
	-0.0003 0.0009 0.0007	-0.003 0.001 0.0007
		(a) HP		(b) Dell

Table 11 :

 11 The covariance matrices calculated from the estimators of HP and Dell laser printers.

	Actual/Predicted HP laser Dell laser
	HP laser	20	0
	Dell laser	0	20

Table 12 :

 12 Confusion matrix for ML classification of two laser printers.

Acknowledgments Q-T. Nguyen wants to dedicate this work to Prof. Y. Delignon, 2nd author of the paper who passed away before this work can be published.


Samia Moumade
email: samia.moumade@iae-aix.com

Aurélie Hemonnet-Goujot

Pierre Valette-Florence

THE ERA OF "WOKE-WASHING": HOW CONSUMERS PERCEIVE BRAND ACTIVISM?

HAL is

PERÇOIVENT L'ACTIVISME DE MARQUE ?

Résumé : Alors que la recherche consacrée à l'activisme de marque (AM) se développe, nous ne savons pas comment les consommateurs perçoivent l'AM et son influence sur la relation avec les marques. Cette recherche exploratoire vise à comprendre la perception et les réactions des consommateurs face aux campagnes d'AM à travers une double perspectivecognitive et affective. Premièrement, il contribue à mieux délimiter la définition d'AM de marque à partir de construits étroitement liées. Deuxièmement, basée sur une technique projective AOL, d'un côté affectif, cette étude révèle le rôle potentiel de la distance psychologique dans la perception par le consommateur des efforts de la marque soutenant une cause sociale et souligne la légitimité de la marque acquise en s'engageant dans l'AM. Du côté cognitif, cette étude montre les incitations managériales derrière les marques à s'engager dans l'AM, tout en offrant aux consommateurs un moyen d'exprimer leur identité et leurs valeurs.

Mots clés : Activisme de marque; Relation consommateur-marque; Qualité de la relation avec la marque

THE ERA OF "WOKE-WASHING": HOW CONSUMERS PERCEIVE BRAND ACTIVISM?

Abstract: While research dedicated to brand activism is growing, we still do not know how consumers perceive brand activism (BA) and how it influences their relationship with brands. This exploratory research aims at gaining a better understanding of the perception and reactions of consumers towards BA campaigns through a dual perspective -cognitive and affective. First, it contributes to better delineate the construct of BA from closely related constructs. Second, based on an AOL projective technique, from an affective perspective, this study reveals the potential role of psychological distance in consumer's perception of the brand's efforts in supporting a social cause. It underlines the acquired brand legitimacy through engaging in brand activism. From a cognitive perspective, this study shows the managerial incentives behind brands in engaging in brand activism, while also providing a way for consumers to express their identity and values. Keywords: Brand Activism; Consumer-Brand Relationship; Brand Relationship Quality