O Mareschal

S Loiseau

F Verjus

L Valbin

G Lissorgues

R Bouregba

G Poullain

S Saez

C Dolabdjian

MODELING AND FABRICATION OF PIEZOELECTRIC ALUMINUM NITRIDE RESONATOR AND ITS APPLICATION IN OSCILLATORS

Keywords: Piezoelectric resonator, oscillator, aluminum nitride, elongation mode

 operating at MHz frequencies. This resonator is composed of a piezoelectric Aluminum Nitride (AlN) layer sandwiched between two aluminum (Al) electrodes. The resonators are manufactured on a silicon substrate. In opposition to common resonator, TFEAR works in extensional (elongation) mode excited via AlN d 31 piezoelectric coefficient. A Finite Element Analysis is performed to simulate the static and modal behaviors. The simulation is used to determine the geometry of each material so as to reach the desired frequency range and to compare theoretical and experimental results at different resonance frequencies. A modified Van-Dyke Butterworth (MBVD) equivalent circuit model is developed, including substrate losses. Samples have been measured and quality factor as high as 2500 in air have been measured. A first oscillator based on TFEAR resonator was designed and fabricated. Results are discussed.

INTRODUCTION

More and more industrial electronic applications, like Universal Serial Bus (USB), Radio Frequency Identification (RFID), or Personal Digital Assistant (PDA), needs smaller, lower power consumption and cheaper systems. Usually, quartz is used for clock reference oscillators. Conventional quartz resonators are millimeter sized, and are not silicon compatible. A solution is to develop Micro-Electro-Mechanical Systems (MEMS) resonators which are micrometer sized and mostly integrated circuit compatible. Bulk Acoustic Wave (BAW) resonators are now used in some applications in the GHz range frequencies. The aim of TFEAR is to be an alternative for quartz crystal resonators in the MHz range frequencies. MHz frequencies electrostatic resonators have been presented in earlier papers [START_REF] Nguyen | MEMS technology for timing and frequency control[END_REF]. Disadvantage of those solutions is the high voltage value needed to polarize the resonators. Piezoelectric resonators do not need this voltage and consequence is lower power consumption.

In this paper, we first address the TFEAR principle, particularly its excitation mode. Physical model will be presented in the first section. The second section deals with the fabrication process. Process flow will be described with a particular attention for the AlN deposition. Then, the MBVD equivalent circuit, used to mimic the TFEAR behavior, will be presented. Finally, results of characterization for oscillator application will be shown and compared to theory and simulation.

PRINCIPLE Thin Film Elongation Acoustic Resonator (TFEAR)

TFEAR consists of a multilayer beam constituted of insulator-metal-piezoelectric-metal materials. Transverse piezoelectric coefficient d 31 of the AlN piezoelectric layer is used [START_REF] Tsubouchi | AlN material constant evaluation and saw properties on AlN/AlO/ and AlN/Si[END_REF]. Consequently, a driving voltage applied on the piezoelectric thin film allows beam elongation. The name "TFEAR" comes from the propagation of bulk acoustic waves along the length of the beam. TFEAR consists of AlN film (1 µm) sandwiched between two Al electrodes (0.2 µm each). AlN material has those advantages to be a non-contaminant piezoelectric material (in opposition to the lead present in the PZT) and presents high acoustic wave velocity [START_REF] Muralt | Is there a better material for thin film BAW applications than AlN?[END_REF]. Two designs were investigated during this study. The first consists of a T-shaped beam anchored at its base end and the second consists of an X-shaped beam anchored at each base ends, as it is shown in Figure 1. For each design, the mask allows to make many samples which present various lengths and widths between 75 µm and 350 µm, and 20 µm or 50 µm respectively. The range of the frequency is tuned by the length of the beam forming the resonator.

Design

Process steps

The resonator is manufactured on (100) oriented high resistivity silicon wafers. TFEAR is fabricated using conventional integrated circuit manufacturing tools. Figure 2 presents a schematic process flow. First, the silicon substrate is thermally oxidized (a). Pads (AlSi 1% Cu 0.04% ) are deposited by DC-magnetron sputtering tool (b). Bottom electrode (AlSi 1% Cu 0.04% ) is also deposited by DCmagnetron sputtering tool (c). Thereafter, AlN layer is deposited by DC or pulsed DC reactive sputtering in nitrogen and argon atmosphere (d). Top electrode is deposited like bottom electrode (e). These four layers are patterned by dry etching for the AlSi 1% Cu 0.04% and by wet etching for the AlN layer. Finally, TFEAR is released by designing of the cavity (f) using anisotropic and isotropic dry etching for the silicon oxide layer and the silicon substrate, respectively. The Figure 3 shows an example of a SEM view of a TFEAR fabricated on the above conditions.

AlN deposition

Regarding the piezoelectric AlN layer, under-layers have to be optimized to have the best crystallized material. Indeed, the most crystallized material leads to the highest piezoelectric coefficient d 31 characterizing the beam.

Thermal oxide shows better results compared to oxide grown by Low Pressure Chemical Vapor Deposition (LPCVD) and oxide grown by Plasma Enhanced CVD (PECVD). Then thermal oxide was chosen for the insulator layer. Designs of experiments were used to optimize the deposition parameters of AlSi 1% Cu 0.04% and AlN layers in order to get the highest crystallinity of AlN with DC or pulsed-DC magnetron sputtering. The aim was achievement of the best crystallinity of the AlN characterized by X-Ray Diffraction (XRD) and thus, the best piezoelectric property and the highest Q factor with the complete resonator. The best results were achieved with the pulsed-DC magnetron sputtering tool (OERLIKON clusterline 200).

PHYSICAL MODEL

Simulations, made with ANSYS software, have been performed to analyze different designs and to predict expected resonance frequencies. The simulation takes into account the electrodes and the part of the substrate (Si) where beam is anchored.

According to geometrical parameters, the reachable resonance frequency f s , of such sample is in the 10 MHz-50 MHz range according to the relation:

1 2 s E f L ρ = ( 1
)
Where L is the length of the beam, E is equivalent Young modulus of the stacked layers, and ρ is the beam equivalent density. These two properties are balanced with the thickness of each layer and can be expressed by: ( )

top

ELECTRICAL MODEL

To describe resonator electrical behaviors, the model based on Van-Dyke Butterworth equivalent circuit used [START_REF] Lakin | Equivalent Circuit Modeling of Stacked Crystal Filters[END_REF] is presented in Figure 4. R m , L m and C m represent the equivalent resistance, inductance and capacitance of the motional branch and C 0 the static capacitance of TFEAR electrodes. Additional parallel resistor R p represents the dielectric losses of the AlN layer. C sub and R sub represent substrate and connection pads elements.

C sub can be evaluated by estimating the capacitance of pads connections and substrate and by testing the substrate behavior with DC polarization. The quality factor Q of the equivalent circuit is given by:

1 m m m L Q R C = . ( 4
)
Figure 4: MBVD equivalent circuit of TFEAR.

CHARACTERIZATION AND RESULTS

Characterization of resonators

On-wafer level measurements are performed using semi automated prober Electroglass 2001X with thermal chuck. Impedance of TFEAR is measured using Agilent 4194A with impedance probe kit. TFEAR is driven by a 25mV signal to work in linear state. The resonance frequencies were previously estimated by simulation which permits to sweep over a frequency range around those values for each sample. The amplitude and phase of the impedance is recorded. Because one of the pads and the substrate are connected to the ground, the electrical model presented in the Figure 4 can be simplified, considering R' p =R sub ||R p and C' 0 =C sub ||C 0 represent the total resistance and capacity of the resonator measured on the substrate, respectively.

Matlab post-processing allows to extract equivalent circuit parameters and fit is performed. The Figure 5 presents measured and fitted curves for a 350 µm length, 50 µm width, T-shaped resonator using SiO 2 insulating layer. Table 1 shows results for a TFEAR, quality factor Q is calculated using formula 4.

Resonance frequency and geometry

Figure 6 shows the resonance frequency as a function of the inverse length of T-shaped TFEAR. Resonance frequencies from 10 MHz to 50 MHz have been measured. We can notice the well matching between measurements and simulation results for wafer 1. Wafer 2 results show a good agreement with theory calculated from equation 1. During the fabrication of TFEAR, alignment errors could occur and some other process problems could appear for each step, which can explain curve divergences for higher frequencies (lower length). Other parameters, like inductive and resistive components of the access line, are not represented and the piezoelectric coefficients used to simulate the AlN layer are theoretical values [START_REF] Muralt | Is there a better material for thin film BAW applications than AlN?[END_REF]. Work is in progress to include these new parameters in the model and to characterize thin film used. This will improve fitting between theory, simulation and measurements.

R m L m C m C 0 R p R sub L sub R sub L sub

OSCILLATOR APPLICATION

A first simple oscillator based on TFEAR resonator was designed and tested as illustrated on Figure 7: TFEAR based oscillator. Experimental values are: A=3 and R=680 Ω. The Barkhausen oscillation condition is:

1 AB = (6)
where A and B are gains of the amplifier and back loop respectively. B is given by:

TFEAR R B R Z = + (7)
where Z TFEAR is the impedance of the resonator and R is the bridge resistance. The oscillation frequency is given by arg(A B)=0 and the gain is fixed by |A B=1|. Assuming << , the solutions can be reduced to:

2 0 1 1 1 2 2 1 m osc m m m m C' R L L C R A R ω π ⎧ ⎛ ⎞ = + ⎪ ⎜ ⎟ ⎪ ⎝ ⎠ ⎨ ⎪ = + ⎪ ⎩ (8)
Oscillator's amplifier is made using an AD8138 operational amplifier with an experimental gain of 3, and then AB > 1 with all TFEAR tested. Used TFEAR resonator is a X-shaped 325 µm length and 50 µm width. We can notice frequency oscillation around 14.0540 MHz closed to the expected performance (f osc =14.11 MHz).

CONCLUSION

Thin film acoustic resonator operating on extensional mode in the range between 10 MHZ to 50 MHz has been fabricated and characterized. Performances of the TFEAR were improved with a Q factor up to 2500 and a motional resistance lower than 400 Ω. The model, which takes into account electric substrate losses, enables to design a first oscillator. Oscillations were obtained, in the expected range, which demonstrates the working ability of these resonators. Study is in progress to confirm the potential of TFEAR to complete with quartz resonator for many applications requiring small size and CMOS compatible process.

Figure 1 :

 1 Figure 1: Schematic view of T-shaped (a) and X-shaped (b) TFEAR.

Figure 1

 1 Figure 1(a) and Figure 1(b) present the piezoelectric layer shape. The excitation voltage is applied on Al electrode deposited on the entire piezoelectric layer.

Figure 2 :

 2 Figure 2: TFEAR cross section process flow views.

Figure 3 :

 3 Figure 3: SEM view of a 75 µm length and 20 µm width TFEAR.

Figure 5 :

 5 Figure 5: Examples of measured and fitted impedance amplitude and phase.

Figure 6 :

 6 Figure 6: Resonance frequency as a function of resonator length.

  The temperature dependency of the frequency resonance is evaluated by extracting the figure of merit, Temperature Coefficient Frequency (TCF), following the relation: variation of frequency with the temperature. For the measured sample (a T-shaped with 350 µm length and 50 µm width), TCF is about -25 ppm/°C between 20 °C and 100 °C. We expect to improve this temperature dependency with SiO 2 insulator layer under bottom electrode to compensate the temperature deviation.

Figure 7 :

 7 Figure 7: TFEAR based oscillator. Experimental values are: A=3 and R=680 Ω.

  bottom

	E E =	AlN	×	t	AlN	+	Al E t	Al	+	t	Al	(2)
	ρ ρ =	AlN	×	t	AlN	+	ρ	Al	t	top Al	+	t	bottom Al	(3)
	where, AlN t ,

and ( ) top Al t , and bottom Al t are the thickness of AlN, Al top and Al bottom layers respectively. E AlN and E Al are the Young modulus, ρ AlN and ρ Al , the mass density of AlN and Al layers respectively.

Table 1 :

 1 Measured MBVD elements, resonance frequency and quality factor of a fabricated TFEAR.

	TFEAR 350 µm
	R m (Ω)	376.4
	L m (mH)	9.6
	C m (fF)	12.8
	C' 0 (pF)	1.88
	R' p (k Ω)	8.9
	f s (MHz) 14.3505
	Q	2300

ACKNOWLEDGMENT

The authors acknowledge J. Rouessard and A. Fougerat (NXP Semiconductors) for their support on characterization and ANSYS simulations, the clean room and characterization laboratory staff at NXP Semiconductors.

CONTACT

Laurie Valbin, tel: +33-145-926-647; l.valbin@esiee.fr TFEAR A R


Matthieu Urvoy

Dalila Goudia

Florent Autrusseau

Perceptual DFT watermarking with improved detection and robustness to geometrical distortions

Keywords: Watermarking, visibility, robustness, contrast sensitivity, Fourier, subjective experiment, Grubbs' test I

More than ever, the growing amount of exchanged digital contents calls for efficient and practical techniques to protect intellectual property rights. During the past two decades, watermarking techniques have been proposed to embed and detect information within these contents, with four key requirements at hand: robustness, security, capacity and invisibility. So far, researchers mostly focused on the first three, but seldom addressed the invisibility from a perceptual perspective and instead mostly relied on objective quality metrics.

In this paper, a novel DFT watermarking scheme featuring perceptually-optimal visibility versus robustness is proposed. The watermark, a noise-like square patch of coefficients, is embedded by substitution within the Fourier domain; the amplitude component adjusts the watermark strength, and the phase component holds the information. A perceptual model of the Human Visual System (HVS) based on the Contrast Sensitivity Function (CSF) and a local contrast pooling is used to determine the optimal strength at which the mark reaches the visibility threshold.

A novel blind detection method is proposed to assess the presence of the watermark. The proposed approach exhibits high robustness to various kind of attacks, including geometrical distortions. Experimental results show that the robustness of the proposed method is globally slightly better than state-of-the-art. A comparative study was conducted at the visibility threshold (from subjective data) and showed that the obtained performances are more stable across various kinds of contents.

Watermarking is one such technique and has been extensively studied for the past two decades; applied to still images, it comes down to embedding an invisible information, called watermark, that can be retrieved and matched even when the watermarked image was attacked to some degree.

Four key requirements have been driving researchers in designing watermarking algorithms: the invisibility, the robustness, the capacity and the security. Any watermarking algorithm should ideally provide the best tradeoff between these four aspects. Lately, security has been widely studied [START_REF] Cayre | Watermarking security: theory and practice[END_REF]; as for the invisibility, the robustness and the capacity, they influence each other and often must be addressed together. For instance, when the robustness is increased, the perceptual quality of the watermarked image inevitably decreases.

Yet, few works address the problem of invisibility from a perceptual perspective. Watson's visual models [START_REF] Watson | DCT quantization matrices visually optimized for individual images[END_REF], [START_REF] Watson | Visual thresholds for wavelet quantization error[END_REF], in particular, have been used to compute Just Noticeable Differences (JND) masks, thus providing a perceptual model of the visibility of the watermark [START_REF] Podilchuk | Image-adaptive watermarking using visual models[END_REF], [START_REF] Tsai | Blind Wavelet-based Image Watermarking Based on HVS and Neural Networks[END_REF] which was embedded either in the Discrete Cosine Transform (DCT) or the Discrete Wavelet Transform (DWT) domains. More recently, another Human Visual System (HVS) model was used for DWT embedding [START_REF] Chou | A perceptually tuned watermarking scheme for color images[END_REF]. Some of these methods are non-blind [START_REF] Podilchuk | Image-adaptive watermarking using visual models[END_REF] or require side information to be transmitted in order to reconstruct the JND mask at the receiver [START_REF] Tsai | Blind Wavelet-based Image Watermarking Based on HVS and Neural Networks[END_REF], [START_REF] Chou | A perceptually tuned watermarking scheme for color images[END_REF]. Some other perceptual methods are fully blind but they use heuristics instead of HVS models to derive JND masks [START_REF] Bartolini | Mask building for perceptually hiding frequency embedded watermarks[END_REF], [START_REF] Voloshynovskiy | A stochastic approach to content adaptive digital image watermarking[END_REF], [START_REF] Mairgiotis | Bayesian watermark detection and new perceptual mask based on a spatially weighted Total Variation image prior[END_REF].

Alternatively, some approaches based on statistical Objective Quality Metrics (OQMs) tune the embedding strength based on the computed quality score [START_REF] Zhang | Spread spectrum image watermarking based on perceptual quality metric[END_REF], [START_REF] Nguyen | Statistical Analysis of Image Quality Metrics for Watermark Transparency Assessment[END_REF] and do not require any side information at the receiver either.

However, two main problems emerge: these metrics often provide wrong estimates of the perceived quality; more importantly, OQMs provide a continuous quality scale and were not intended to scenarios, such as watermarking, targeting the visibility threshold. In practice, the only reliable way to assess the invisibility is to conduct a subjective experiment in which observers are asked whether they can notice the watermark or not. In addition, most comparisons between watermarking techniques found in the literature are also based on OQMs -most often the Peak Signal to Noise Ratio (PSNR) -which, as will be shown in this paper, can be misleading.

Besides invisibility, the robustness is also a key aspect in watermarking. None of the aforementioned techniques are robust to common attacks such as geometrical distortions: they cannot efficiently withstand their desynchronization effects. Rotation-Scale-Translation (RST) invariant techniques have been proposed to address this issue through the use of appropriate transform domains. In [START_REF] O'ruanaidh | Rotation, scale and translation invariant digital image watermarking[END_REF], the watermark is embedded in the Fourier-Mellin domain, to the cost of an important algorithmic complexity. Later, in [START_REF] Lin | Rotation, scale, and translation resilient watermarking for images[END_REF], log-polar mapping of the Fourier domain is used to make the watermark robust against RST attacks. Another log-polar Fourier embedding technique is proposed in [START_REF] Kang | Efficient General Print-Scanning Resilient Data Hiding Based on Uniform Log-Polar Mapping[END_REF]. Besides these log-polar mapping techniques -that commonly exhibit an important computational load -, some directly embed the watermark within the modulus of the Discrete Fourier Transform (DFT). In [START_REF] Solachidis | Circularly symmetric watermark embedding in 2-D DFT domain[END_REF], the watermark is made rotation invariant thanks to its circular shape. A similar technique is proposed in [START_REF] Poljicak | Discrete Fourier transformbased watermarking method with an optimal implementation radius[END_REF] where circular dots are embedded in the Fourier transformed coefficients.

The success of a watermarking algorithm not only depends on its embedding strategy, but also on its detection strategy. In particular, properly estimating the optimal detection threshold is crucial and thus has been extensively studied. Usually, the problem is stated in terms of two hypotheses H1 and H0, which respectively denote scenarios May 6, 2014 DRAFT with and without embedded watermark in the considered host content. In [START_REF] Piva | Threshold Selection for Correlation-Based Watermark Detection[END_REF], [START_REF] Barni | Improved wavelet-based watermarking through pixel-wise masking[END_REF], statistical models for the signal and the watermark are used to derive a theoretical threshold. Typically, these are Gaussian or binomial models; the threshold is chosen so as to limit the probability of false detections. However, the statistics of both the watermark and the host image are likely to change with various parameters (e.g. watermark size, embedded message, etc), thus impacting the parameters of the fitted models, which in turn requires for the detection threshold to be re-evaluated.

Moreover, if the statistics of both the host image and the watermark do not quite strictly fit the derived models, the obtained threshold is likely not to be optimal.

In this work, a robust, blind, substitutive and perceptual watermarking technique is proposed. Various characteristics of the HVS are used to determine and adjust the visibility level of the embedded watermark, thus resulting in an optimal invisibility versus robustness tradeoff. The proposed technique is designed to be robust against various kinds of attacks (including geometrical distortions). In addition, a subjective experiment is conducted in order to assess the invisibility of the watermark. Finally, the proposed detection strategy is efficient and adapts to the statistics of the watermark.

II. OVERVIEW OF THE CONTRIBUTION

Similarly to [START_REF] Solachidis | Circularly symmetric watermark embedding in 2-D DFT domain[END_REF], [START_REF] Poljicak | Discrete Fourier transformbased watermarking method with an optimal implementation radius[END_REF], it is proposed to embed the watermark in the Fourier domain; the magnitude is used to control the energy of the watermark while the phase is used to hold its information. For best robustness, the watermarked coefficients are grouped into two symmetrical square patches, which can be expressed as a sum of sinewaves. Once displayed, this results in a combination of sine gratings at various visual frequencies and orientations: a perceptual model is used to adjust their amplitude at the visibility threshold. The watermarked image is then obtained by inverse Fourier transform. Sec. III describes the perceptual model, then Sec. IV focuses on the embedding process. Sec. V details the detection process.

The proposed method is designed to be robust against multiple kinds of attacks. Sec. VI investigates the robustness of the proposed method in comparison to [START_REF] Solachidis | Circularly symmetric watermark embedding in 2-D DFT domain[END_REF] and [START_REF] Poljicak | Discrete Fourier transformbased watermarking method with an optimal implementation radius[END_REF]. Sec. VII reports the results of a subjective experiment assessing the visibility of the watermark Sec. VIII then re-evaluate the robustness when the strength has been adjusted according to the observers feedback. Finally, Sec. IX validates the proposed approach in the context of low quality Print & Scan (P&S).

III. REACHING THE VISIBILITY THRESHOLD

Choosing the appropriate watermark strength is a delicate but crucial step while designing a watermarking technique. Over-estimations of the strength are likely to result in visible watermarks; on the contrary, underestimations are likely to reduce the robustness performances. In this paper, it is proposed to use some properties of the HVS to automatically determine the perceptually optimal watermarking strength, at which the embedded watermark appears at the visibility threshold.

Computational models of the HVS providing estimates for the visibility of differences between an original and a distorted image have been proposed, such as in [START_REF] Barten | Evaluation of subjective image quality with the square-root integral method[END_REF], [START_REF] Daly | Visible differences predictor: an algorithm for the assessment of image fidelity[END_REF], [START_REF] Mantiuk | HDR-VDP-2[END_REF], [START_REF] Watson | Model of visual contrast gain control and pattern masking[END_REF] and [START_REF] Watson | A standard model for foveal detection of spatial contrast ModelFest experiment[END_REF]. They generally implement some embed ourier rotation is the depends In threshold Usually, and without [START_REF] Lin | Rotation, scale, and translation resilient watermarking for images[END_REF], atermark these

I sRGB Gamma expansion Γ RGB to CIE XYZ Y Foveal filter Φ 1/x Fourier transform F (f k , θ k ) CSF (f k , θ k ) Psycho- metric function Pooling f k ,θ k Ψ mark

Viewing conditions

I RGB Y mean Y local F Y C τ * local Ψ(C τ * local )
Fig. 1. The proposed HVS model estimates the probability Ψ mark that a watermark -embedded in an image I sRGB at visual frequencies f k and orientations θ k -is visible. Please refer to Section III for employed notations.

of the following elements: (1) non-linear sensitivity to the amplitude of the luminance, (2) conversion to contrast,

(3) contrast sensitivity to visual frequencies, (4) oblique effect, (5) subband decomposition into visual channels, [START_REF] Chou | A perceptually tuned watermarking scheme for color images[END_REF] masking effects, (7) psychometric probability of detection and (8) error pooling. A survey of these models can be found in [START_REF] Zhang-Ge | Numerical Observers for the Objective Quality Assessment of Medical Images[END_REF].

Practical applications, however, may not implement all of these HVS properties; visual channels, for instance, may be discarded due to their high complexity and their rather small contribution to the visibility estimates [START_REF] Watson | A standard model for foveal detection of spatial contrast ModelFest experiment[END_REF].

In a preliminary subjective experiment involving 4 expert observers, it was reported that the proposed watermark is first noticed in uniform image areas where masking effects do not occur (see Sec. VII-B for experimental details).

Therefore, the proposed HVS model is simplified and discards the perceptual channel decomposition. Moreover, excluding the masking effects from the model underestimates the visibility threshold and thus the watermark strength: the invisibility requirement is met anyway. The proposed model is illustrated in Fig. 1.

A. Modeling viewing conditions

Typically, HVS models require both the viewed image and the viewing conditions to be input. Let I sRGB (x, y)

denote an image to be watermarked, 0 ≤ x<R x , 0 ≤ y<R y , where R x and R y are respectively its horizontal and vertical resolutions. Let S x and S y denote I sRGB 's displayed size in meters. The display illumination is noted L; it is set to 280 cd.m -2 in the proposed model, a typical value for modern LCD monitors (slightly above the advised illumination level in [START_REF] Bt | General viewing conditions for subjective assessment of quality of SDTV and HDTV television pictures on flat panel displays BT Series Broadcasting service[END_REF]). Michelson's formula for contrast [START_REF] Michelson | Studies in Optics[END_REF], when applied to a sine grating of peak amplitude A peak , defines the contrast as

C(x, y)= A peak Y mean (x, y) , (1)
where Y mean is the mean illumination of the area supporting the sine grating. In typical images, the illumination varies locally; Y mean thus needs to be computed locally as well. As proposed in [START_REF] Frazor | Local luminance and contrast in natural images[END_REF], a raised cosine filter Φ of diameter one degree of visual angle is applied to Y to obtain Y mean (x, y); this provides an estimate of the average foveal illumination over the entire image domain. The Fourier domain, which is later used to embed the watermark, is not spatially localized. Prior to applying such non-localized processing steps, it is proposed to normalize Y with respect to the local luminance Y mean . The locally normalized luminance map Y local (x, y)= Y (x,y) /Ymean(x,y) may then be input to Fourier computations. Within this normalized luminance space, the local contrast of a sine grating writes as

C local = A peak . (2)

C. Contrast sensitivity

The sensitivity to contrast is influenced by numerous factors [START_REF] Barten | Contrast sensitivity of the human eye and its effect on image quality[END_REF]. Intensive research has made available a large number of computational models known as Contrast Sensitivity Functions (CSFs), describing our sensitivity to contrast levels as a function of the visual frequency: the CSF returns the inverse of the threshold contrast above which a sine grating becomes visible. In [START_REF] Barten | Contrast sensitivity of the human eye and its effect on image quality[END_REF], an accurate CSF is built by interlinking models for the optical Modulation Transfer Function, the photon noise, the lateral inhibition, and the neural noise. Some aspects of this model, though, are not theoretically supported; authors in [START_REF] Roka | Contrast sensitivity model of the human eye[END_REF] addressed these issues and incorporated better models for cone and midget ganglion cell densities. However, the resulting models are computationally intensive and do not fit practical applications such as ours.

CSFs with medium to low complexities have been made available as well. In [START_REF] Mannos | The effects of a visual fidelity criterion of the encoding of images[END_REF], the CSF proposed by Mannos and Sakrison features a single parameter, the visual frequency. However, such a model lacks important factors such as light adaptation and stimulus size. In, the CSF proposed by Daly features key factors, but was designed to be used within his Visible Differences Predictor (VDP) quality model and should be used with caution in other contexts.

In [START_REF] Barten | Formula for the contrast sensitivity of the human eye[END_REF], finally, Barten provides a simplified formula for his initial CSF, which also incorporates the oblique effect May 6, 2014 DRAFT and the influence of the surround luminance. In this paper, the proposed Fourier watermark embedding technique (as will be further discussed in section IV), modifies frequency coefficients whose orientations are oblique. For this reason, it is necessary to implement the oblique effect as our sensitivity varies with the orientation of the visual pattern. Barten's simplified CSF formula [START_REF] Barten | Formula for the contrast sensitivity of the human eye[END_REF], at binocular viewing, will thus be used in the proposed model:

CSF(f, θ) = 5200 • e -0.0016 • f 2 (1+ 100 L ) 0.08 • 1+ 144 Θ 2 (I) +0.64 • f 2 • ⇣ 1 + 3 sin 2 (2θ) ⌘ ! -0.5 • ✓ 63
L 0.83 + 1 1 -e -0.02f 2 ◆ -0.5 (3)
where f is the visual frequency in cycles per degree (cpd), L is the adaptation luminance in cd.m -2 and is assumed to be equivalent to the display illumination (see Sec. III-A), Θ 2 (I) is the square angular area of the displayed image I in square visual degrees, and θ the orientation angle.

From Eqs. ( 2) and ( 3), one may now obtain the threshold amplitude A τ peak (f, θ) of a sine grating

A τ peak (f, θ)=C τ local (f, θ)= 1 CSF(f, θ) (4)
where C τ local is the local contrast threshold.

D. Psychometric function

The psychometric function Ψ is typically used to relate the parameter of a physical stimulus to the subjective responses. When applied to contrast sensitivity, Ψ may describe the relationship between the contrast level and the probability that such contrast can be perceived [START_REF] Daly | Visible differences predictor: an algorithm for the assessment of image fidelity[END_REF]. In the proposed approach, Daly's Weibull parametrization [START_REF] Daly | Visible differences predictor: an algorithm for the assessment of image fidelity[END_REF] is used:

Ψ C τ ⇤ local =1-e -(C τ * local ) β , (5)
where

C τ ⇤ local = Clocal /C τ local
is the ratio between the locally normalized contrast and its threshold value given in Eq. ( 4). β is the slope of the psychometric curve; its value is usually obtained by fitting from experimental data. Typical values for β range from 1.3 (contrast discrimination) up to 4 (contrast detection) [START_REF] Wallis | The slope of the psychometric function and non-stationarity of thresholds in spatiotemporal contrast vision[END_REF]. Moreover, the psychometric function is typically applied locally, contrary to our model which applies it in the Fourier domain, hence globally.

For both these reasons, it is proposed to use β =2 , a rather low value which provides a large transition area between invisible and visible domains.

E. Watermark frequency pooling

In this paper, the watermark is embedded into multiple Fourier coefficients, which results into the superimposing of multiple sine gratings G k , 0 ≤ k<N, with visual frequencies f k and orientations θ k . The CSF solely provides an estimate for the visibility of a single sine grating; a summation model is thus required to estimate the combined visibility level of all embedded gratings. Numerous models have been proposed [START_REF] Watson | A standard model for foveal detection of spatial contrast ModelFest experiment[END_REF], [START_REF] Daly | Visible differences predictor: an algorithm for the assessment of image fidelity[END_REF]. Recent findings suggest [START_REF] Mantiuk | HDR-VDP-2[END_REF] DRAFT that Minkowski summation [START_REF] Watson | A standard model for foveal detection of spatial contrast ModelFest experiment[END_REF] may not be adequate [START_REF] Meese | Theory and data for area summation of contrast with and without uncertainty: evidence for a noisy energy model[END_REF]; the proposed model will thus use probability summation as in [START_REF] Daly | Visible differences predictor: an algorithm for the assessment of image fidelity[END_REF].

Let C τ ⇤ local (f k , θ k ) denote the contrast threshold ratio of G k .
As in [START_REF] Daly | Visible differences predictor: an algorithm for the assessment of image fidelity[END_REF], the probability Ψ mark that the watermark is perceived by an observer computes as

Ψ mark =1- N -1 Y k=0 ⇣ 1 -Ψ C τ ⇤ local (f k , θ k ) ⌘ , (6)
where Ψ mark is set according to Eq. ( 5): 1 -e -(1) β =0.632.

Let us now work our way back from the entire watermark to individual gratings. In the proposed approach, it is assumed that all embedded gratings share the same probability of detection

Ψ C τ ⇤ local (f k , θ k ) = ψ, ∀k ∈ [0,N[
. This probability can be derived from Ψ mark by solving Eq. ( 6):

ψ =1-(1 -Ψ mark ) 1/N . ( 7
)
The corresponding normalized contrast, which is equal to the optimal grating peak amplitude, may thus be obtained by inversion of the psychometric function (Eq. ( 5))

A peak (f k , θ k )=C τ local (f k , θ k ) •  -ln ⇣ (1 -Ψ mark ) 1/N ⌘ 1/β (8)

IV. MAGNITUDE-PHASE SUBSTITUTIVE EMBEDDING

The proposed embedding technique relies on both the amplitude and the phase Fourier components to modulate a binary watermark. The amplitude allows to adjust the watermark strength whereas the phase holds the binary information.

Let an N-bits binary sequence M 1..N be the message to embed. A binary watermark W(i, j), 0 ≤ i<M , 0 ≤ j<M, is generated from M 1..N using a PRNG and arranged into a M × M matrix.

Let F Y (u, v) denote the Fourier transform of Y local over frequency domain

Ω = (u, v):-Rx /2 ≤ u< Rx /2, -Ry /2 ≤ v< Ry /2 ,
where u and v are the horizontal and vertical frequencies. Let

Ω + W = (u, v):u w ≤ u<u w + M, v w ≤ v<v w + M } and Ω - W = n (u, v):(-u, -v) ∈ Ω + W o
be two subsets of Ω within which the watermark will be embedded; u w and v w are called watermark modulation frequencies. The watermarked spectrum is obtained by substitution as follows:

f F Y (u, v)= 8 > > > > > < > > > > > : α /2 • A ⇤ peak (u, v) • e π • W(u-uw,v-vw) , (u, v) ∈ Ω + W α /2 • A ⇤ peak (u, v) • e π • W(-u-uw,-v-vw) , (u, v) ∈ Ω - W F Y (u, v), elsewhere (9)
where α is a weighting factor that controls the watermark energy relatively to the predicted visibility level. A ⇤ peak (u, v) equals to A peak f (u, v), θ(u, v) and denotes the embedded grating's optimal amplitude obtained in Eq. ( 8) for discrete frequency (u, v); f (u, v) and θ(u, v) respectively denote its visual frequency and orientation. In this paper, the values for u w and v w are both set to 80% of the highest horizontal and vertical Fourier frequencies, thus ensuring low sensitivity to watermarked frequencies, which induces high threshold amplitudes, hence high watermark energy.

Lower values for u w and v w were also experimented (20% of the maximum frequency), with similar robustness performances.

The proposed embedding scheme given in Eq. ( 9) is a substitutive technique. Here, the choice for amplitude substitution is straightforward as it allows precise control over the energy of the watermark. Moreover, as can be seen in Eq. ( 9), phase substitution is performed in order to encode the binary watermark information: zeros are coded as null phases, ones as π phases. Eventually, the watermarked image is obtained by inverse Fourier transform, local luminance de-normalization, and transform into the original color space. watermark to be detected W(i, j). Then, outlier detection is performed to assess the presence of matching locations within the correlation matrix, thus providing a binary decision (zero-bit watermarking).

V. T WO-STEPS

A. A 2D-correlation algorithm based on template matching

Often, detection schemes are solely based on a correlation coefficient that is compared to a predefined threshold value (e.g. [START_REF] Podilchuk | Image-adaptive watermarking using visual models[END_REF], [START_REF] Solachidis | Circularly symmetric watermark embedding in 2-D DFT domain[END_REF]). Detectors based on cross-correlation, cross-covariance and their optimized variants [START_REF] Lewis | Fast normalized cross-correlation[END_REF] May 6, 2014 DRAFT generally perform better, at the cost of an increased computational load. Most often, however, one-dimensional correlation is performed [START_REF] Poljicak | Discrete Fourier transformbased watermarking method with an optimal implementation radius[END_REF], while the searched pattern is likely to be two-dimensional. Moreover, geometric distortions affect the location and shape of the watermark in the Fourier domain.

Not only does 2D cross-correlation account for the two-dimensional structure of the embedded watermark, but the displacement of its correlation peaks allows to estimate the geometrical distortions undergone by the input image [ I sRGB . Therefore, it is proposed to perform two-dimensional TM similarly to [START_REF] Kang | Efficient General Print-Scanning Resilient Data Hiding Based on Uniform Log-Polar Mapping[END_REF]. Furthermore, no additional resynchronisation watermark is required (contrary to [START_REF] Solanki | Estimating and undoing rotation for print-scan resilient data hiding[END_REF], [START_REF] Solanki | Print and Scan' Resilient Data Hiding in Images[END_REF], [START_REF] Pereira | Template based recovery of Fourier-based watermarks using log-polar and log-log maps[END_REF]); the TM is able to retrieve the watermark as long as it is located, even partially, within the searched area.

Typically, TM involves: (1) a supporting signal S, (2) a searched template T , (3) a search area Ω S and (4) a matching criterion ρ. In our case, S is derived from c F Y :

S(u, v)= 8 > > > < > > > : c F Y (u, v) , -π /2 < arg ⇣ c F Y (u, v) ⌘ ≤ π /2 -c F Y (u, v) , π /2 < arg ⇣ c F Y (u, v) ⌘ ≤-π /2 (10)
which encompasses both the Fourier energy and the phase information into a single signed signal. This stands in contrast to Kang et al. [START_REF] Kang | Efficient General Print-Scanning Resilient Data Hiding Based on Uniform Log-Polar Mapping[END_REF], who perform phase-only correlation to cope with the large dynamic range of the energy of the watermarked coefficients. In our scenario, this range depends on the CSF; therefore it is rather low, which allows us to take the energy level into account as well. The template T is derived from the binary watermark W as follows:

∀(i, j) ∈ [0,M[ 2 , T (i, j)=A ⇤ peak (u w + i, v w + j) • 2 •W(i, j) -1 (11)
which maps W's binary values from {0, 1} to {-1, 1} and accounts for the optimal amplitude A ⇤ peak of the embedded gratings (see Eq. ( 9)). The search area Ω S is centered around the original watermark modulation frequencies, such that Ω S = (u, v):u w -M<u<u w + M, v w -M<v<v w + M }. Finally, Pearson's correlation coefficient serves as matching criterion ρ, so that the final 2D correlation matrix is given by ξ(u, v)=ρ(T , S(Ω M ⇥M u,v

)), where

(u, v) ∈ Ω S and Ω M ⇥M u,v = (u 0 ,v 0 ):u ≤ u 0 <u+ M, v ≤ v 0 <v+ M .

B. Peak detection

Typically, watermark detection is granted when the correlation score exceeds some given threshold value. It is generally obtained by measuring [START_REF] Solachidis | Circularly symmetric watermark embedding in 2-D DFT domain[END_REF] (experimentally) or estimating [START_REF] Barni | Improved wavelet-based watermarking through pixel-wise masking[END_REF] (theoretically) the distribution of the values of the correlation matrix under true or false detection scenarios. On the one hand, theoretical approaches are based on statistical analysis, and thus make assumptions that may not hold in practice. On the other hand, experimental thresholds are only valid within the scope of the test signals. In any case, the obtained threshold is constant, and might not be optimal.

In contrast, the dynamic range of a correlation matrix varies with numerous parameters. To illustrate this, the proposed embedding method was used to watermark the Lena image. The watermarked image was then rotated May 6, 2014 DRAFT

(1.5 • ) and blurred (Gaussian noise, σ =2.0) to simulate an attack. Fig. 2 plots the obtained correlation matrices, for M = 32 (Fig. 2a) and M = 64 (Fig. 2b). Although the detection peaks are obvious in both cases, their amplitude, as well as the amplitude of the surrounding noise, differ significantly. In other words, the detection decision should be driven by the relative difference in amplitude between the correlation peak(s) and the surrounding noise. In Fig. 2, for the same image, the same embedding method, and against the same attack, a fixed detection threshold (e.g. 0.4) would properly detect the correlation peak from Fig. 2a and miss the one from Fig. 2b.

Instead, detection peaks may be seen as outliers. Numerous methods for outlier detection have been proposed in the context of statistical analysis [START_REF] Hodge | A Survey of Outlier Detection Methodologies[END_REF]. Grubbs' test [START_REF] Grubbs | Procedures for detecting outlying observations in samples[END_REF] is one such method, both robust and low computational, and solely requires that input correlation matrices are (approximately) normally distributed. This assumption is quite common [START_REF] Barni | Improved wavelet-based watermarking through pixel-wise masking[END_REF]; Kolmogorov-Smirnov tests further confirmed Gaussianity of the correlation matrices in practical experiments. The proposed detection method features an iterative implementation of Grubbs' test that removes one outlier value at a time, up to a predefined maximum number of outliers. Alternatively, the performances of the Extreme Studentized Deviate (ESD) test were also investigated and proved to be identical in practical experiments.

Such an approach prevents us from fixing the threshold, and adapts the detection to the observed correlation matrix.

Moreover, it can be applied to any correlation-based watermark embedding method. Last but not least, the test's significance level α G can be used to control the tradeoff between detection capacity and false alarm rate. High (resp. low) values bring higher (resp. lower) True Positive (TP) and False Positive (FP) rates. In terms of complexity, the proposed method proves to be very fast: with 1024 × 1024 images and 64 × 64 watermarks, both detection and embedding steps took approximately 0.2s to perform on a mid-2012 Mac Book Pro (2.3 GHz Intel Core i7 processor, 8Gb of RAM), of which 0.26ms was needed to apply Grubbs' test.

VI. EXPERIMENTAL RESULTS

Most experimental results were obtained from dataset D a containing nine images with resolution 512 × 768, including six natural color images, two text images and one cartoon color image. They are listed in Table I. When a larger number of images were required (e.g. threshold selection), a second dataset D b featuring 1000 natural color images was used. The results are compared to those of [START_REF] Solachidis | Circularly symmetric watermark embedding in 2-D DFT domain[END_REF] and [START_REF] Poljicak | Discrete Fourier transformbased watermarking method with an optimal implementation radius[END_REF]. Datasets, watermarked images and additional data are available online 1 . Fig. 3 shows examples of watermarked images at default embedding strength α 0 (α =1, see Eq. ( 9)).

A. Grubbs significance level

In order to determine the optimal significance level for Grubbs' test, D a 's images were watermarked with the proposed algorithm. 90 attacks, from the Stirmark benchmark [START_REF] Petitcolas | Attacks on copyright marking systems[END_REF], were then applied to both original (hypothesis H0) and watermarked (hypothesis H1) images. Detection was finally run on each resulting image with varying values for α G . This process was repeated for several watermark sizes M ∈ {16, 32, 64, 128}. The results obtained for both H0 and H1 scenarios are shown in Fig. 4. As can be seen from Fig. 4a, the TP rate slightly varies with α G , and mostly with the size of the embedded watermark. As for the FP rate, it also increases with the size of the watermark, but exponentially increases with α G . Therefore, the value of α G may be adjusted according to the desired false alarm rate. In the next experiments, we set M to 64 -in order to provide high detection performances -and α G to 1e -4 -in order to ensure a low FP rate (null in experimented images) -.

B. False alarms and true detections

To ensure optimal detection performances, the gap between the detection scores in watermarked images (hypothesis H1) and those in un-marked images (hypothesis H0) need to be as large as possible. Detection scores were collected on dataset D b (1000 images). Fig. 5 plots the histograms (solid lines) of the detection scores in both scenarios H0 and H1, along with the fitted Gaussian models (colored areas). As can be seen, H0 and H1 distributions are clearly disjoint in the proposed technique, similarly to [START_REF] Solachidis | Circularly symmetric watermark embedding in 2-D DFT domain[END_REF]. Conversely, Fig. 5 shows that the distributions of H0 and H1 are much closer in [START_REF] Poljicak | Discrete Fourier transformbased watermarking method with an optimal implementation radius[END_REF], thus resulting in respectively higher and lower FP and TP rates. FP rat e (%) 10 -7 10 -6 10 -5 10 -4 10 -3 10 -2 0.05

• • • • • • • • • • • • • • • • • • • • • • • • Mark size M = 16 M = 32 M = 64 M = 128
• • • • • • • • • • • • • • • • • • • • • • • • Mark size M = 16 M = 32 M = 64 M = 128 (

C. Robustness to attacks

The robustness of the proposed watermarking algorithm to various attacks was measured with the Stirmark benchmark [START_REF] Petitcolas | Attacks on copyright marking systems[END_REF]. All 90 Stirmark attacks were applied on dataset D a , thus resulting in 810 attacked images. The robustness is reported in terms of maximum of correlation against experimented attacks; yet, the proposed detection scheme depends on the output of Grubbs' test (absence/presence of outliers). Therefore, the robustness of the proposed algorithm is reported in terms of percentage of images with a positive detection (detection rate). The obtained results are compared with [START_REF] Solachidis | Circularly symmetric watermark embedding in 2-D DFT domain[END_REF] and [START_REF] Poljicak | Discrete Fourier transformbased watermarking method with an optimal implementation radius[END_REF], where detection rates correspond to the percentage of images for which correlation is higher than the predefined thresholds.

The obtained results are plotted in Fig. 6; details for each attack are available online. For each attack, the dark bar(s) correspond to the best performing algorithm(s) and the light bar(s) to the least performing algorithm(s).

Overall, it appears that the proposed technique performs best with an average detection rate of 62.5%, while [START_REF] Solachidis | Circularly symmetric watermark embedding in 2-D DFT domain[END_REF] and [START_REF] Poljicak | Discrete Fourier transformbased watermarking method with an optimal implementation radius[END_REF] respectively reach 52.1% and 53.1%. It should be noted that [START_REF] Poljicak | Discrete Fourier transformbased watermarking method with an optimal implementation radius[END_REF] performs well against most kinds of attacks, although seldom performing best out of the three algorithms. The proposed technique is especially robust to geometric distortions such as rotations (despite the fact that [START_REF] Solachidis | Circularly symmetric watermark embedding in 2-D DFT domain[END_REF] and [START_REF] Poljicak | Discrete Fourier transformbased watermarking method with an optimal implementation radius[END_REF] embed a circular watermark contrary to the proposed approach) and shearing. Still, only [START_REF] Poljicak | Discrete Fourier transformbased watermarking method with an optimal implementation radius[END_REF] withstands severe rotations and cropping. In addition, the proposed technique is nearly as robust as [START_REF] Solachidis | Circularly symmetric watermark embedding in 2-D DFT domain[END_REF] to scaling and filtering operations, and nearly as robust as [START_REF] Poljicak | Discrete Fourier transformbased watermarking method with an optimal implementation radius[END_REF] to cropping. Moreover, it is significantly more robust to random bending, which randomly adds local geometrical distortions, as can be induced by severe P&S effects. On the downside, it appears that the proposed technique is not robust to severe JPEG compression, although this can be explained by the fact that the chosen watermark modulation frequencies are high and thus very likely to be strongly affected by compression. Better robustness against coding artifacts could be achieved by simply shifting the embedded watermark towards lower frequencies.

Further information can be found in Table II which lists average detection rates amongst groups of attacks.

VII. SUBJECTIVE EXPERIMENT

Previously, we assessed the robustness of the proposed watermarking scheme at default strength α 0 =1and compared it to [START_REF] Solachidis | Circularly symmetric watermark embedding in 2-D DFT domain[END_REF], [START_REF] Poljicak | Discrete Fourier transformbased watermarking method with an optimal implementation radius[END_REF] at suggested strengths. Yet, one may wonder whether these strengths are equivalent in terms of visibility, hence whether previous results are truly comparable. For this reason, a subjective experiment was conducted to determine the optimal strengths at which the watermarks appear at visibility threshold. Section VIII will re-assess the robustness of the compared algorithms at the obtained strengths. Det ect ion rat e (%) 0 20 40 60 80 100 [START_REF] Poljicak | Discrete Fourier transformbased watermarking method with an optimal implementation radius[END_REF] 0 20 40 60 80 100 [START_REF] Solachidis | Circularly symmetric watermark embedding in 2-D DFT domain[END_REF] Fig. 6. Robustness to Stirmark attacks at default embedding strength (α 0 ): detection rates on dataset Da

A. Apparatus and methodology

Observers were seated in a standardized room [START_REF]Methodology for the subjective assessment of the quality of television pictures[END_REF] Eye One Pro luminance meter: gamma correction was set to 2.2 and white point to 6600 K. Room illumination was set to 30 cd.m -2 behind the screen, hence 15% of the perceived screen brightness. Finally, the viewing distance was set to six times the height of the presented images (as in [START_REF] Chou | A perceptually tuned watermarking scheme for color images[END_REF]).

The Two Alternative Forced Choice (2AFC) protocol was used as it is best suited to estimate visibility thresholds.

37 naive observers were recruited to participate to the experiment; screening tests ensured that they had a perfect visual acuity (Snellen chart) and no color deficiencies (Ishihara plates). The 2AFC methodology [START_REF] Cox | Digital Watermarking and Steganography[END_REF] (chap. 8, p.

258) was used to assess the visibility of the embedded watermarks. Nine anchor images were first presented for the observers to familiarize themselves with the sources and the type of degradations. All 243 images were then shown randomly to each observer. Observers were showed each pair of images (the watermarked image and its original version) during 10 seconds, at which point images were hidden. Observers were forced to vote for either one or the other image they thought was containing the watermark, before they could resume the experiment.

B. Stimuli set

Source images were taken from dataset D a . The nine images were then watermarked with three different algorithms (proposed, [START_REF] Solachidis | Circularly symmetric watermark embedding in 2-D DFT domain[END_REF] and [START_REF] Poljicak | Discrete Fourier transformbased watermarking method with an optimal implementation radius[END_REF]). This watermarking was repeated for various embedding strengths

α i = w i • α 0 , 1 ≤ i ≤ 9
, where w i are weighting coefficients, and α 0 is the default strength. α 0 = α =1according to Eq. ( 9) in the proposed algorithm; α 0 =0.3 in [START_REF] Solachidis | Circularly symmetric watermark embedding in 2-D DFT domain[END_REF]; finally, α 0 is image-adaptive in [START_REF] Poljicak | Discrete Fourier transformbased watermarking method with an optimal implementation radius[END_REF].

In 2AFC experiments, it is especially important to make sure that presented stimuli span the entire visibility range, therefore from completely invisible to severely distorted: the obtained results can thus be best fitted to the expected psychometric curve. Proposed 1-7 0.13 0.25 0.50 0.75 1.00 1.50 2.00 4.00 6.00 8-9 0.25 0.50 0.75 1.00 1.25 1.38 1.50 2.00 4.00 [START_REF] Solachidis | Circularly symmetric watermark embedding in 2-D DFT domain[END_REF] 1-7 0.06 0.13 0.19 0.25 0.38 0.50 1.00 2.00 3.00 8-9 0.03 0.06 0.09 0.13 0.15 0.17 0.25 0.50 1.00 [START_REF] Poljicak | Discrete Fourier transformbased watermarking method with an optimal implementation radius[END_REF] 1-7 0.25 0.37 0.50 0.75 1.00 2.00 5.00 10.00 20.00 8-9 0.50 0.75 1.00 2.00 3.00 5.00 7.00 10.00 20.00

C. Analysis

Let V(i, j) denote the vote of the j th observer at the i th image: it is equal to 1 if the watermarked image was correctly identified, and to 0 otherwise. For each of the 243 watermarked images, the percentage of observers who correctly identified the watermarked image R(i)= 1 /Nobs • P Nobs j=1 V(i, j) is called the 2AFC detection rate, where N obs = 37 is the number of observers. A clustering-based analysis revealed 6 observers whose votes were inconsistent, these were removed from subsequent analysis, thus reducing N obs to 31. Fig. 7 plots the 2AFC detection rates R(i) obtained for image ED from dataset D a (circles) against the strength weighting coefficients w i . Typically, the visibility threshold is assumed to be obtained for a 2AFC score of 75% [START_REF] Cox | Digital Watermarking and Steganography[END_REF]. Experimental data are fitted to a sigmoidal Weibull curve to obtain the corresponding psychometric curve, from which the threshold value and therefore optimal watermarking strength α ⇤ can be obtained. The same fitting process is repeated for each image and each algorithm, therefore resulting in 27 perceptually May 6, 2014 DRAFT optimized strengths α ⇤ = w ⇤ • α 0 . Table IV lists the obtained weights w ⇤ . Fig. 8 plots w ⇤ 's distributions for each of the three algorithms; for plotting purposes, the histogram was computed on the binary logarithm of w ⇤ : log 2 (w ⇤ )=0 therefore corresponds to α = α 0 , log 2 (w ⇤ )=1to α =2• α 0 and log 2 (w ⇤ )=-1 to α =0.5 • α 0 . These strengths α ⇤ , in the scope of experimented images, optimize the robustness versus visibility tradeoff. The further away the default strength α 0 is from its optimum α ⇤ , the more (in)visible is the watermark. When the default strength α 0 is greater than its optimum α ⇤ (i.e. log 2 (w ⇤ ) < 0), the watermark is likely to become visible. Conversely, when α 0 is less than α ⇤ (i.e. log 2 (w ⇤ ) > 0), the watermark is invisible but the strength is not maximized, hence likely resulting in a loss of robustness. within the [0.46; 3.64] range (Table V). In addition, both of these algorithms are more likely to under-estimate the strength than to over-estimate it, which is preferable as it keeps the watermark invisible. Above all, these results tend to show that the proposed HVS model provides good estimates for the visibility threshold. While [START_REF] Poljicak | Discrete Fourier transformbased watermarking method with an optimal implementation radius[END_REF] does not feature any kind of psychophysical model, the embedded watermark also nears the visibility threshold.

VIII. PERFORMANCES AT THE VISIBILITY THRESHOLD

Here, the algorithms are re-evaluated with the obtained perceptually optimized strengths α ⇤ . Images from dataset D a are thus watermarked using α ⇤ instead of α 0 .

A. Robustness to attacks

The robustness benchmark performed in section VI-C was thus rerun at α ⇤ . The obtained results are plotted in Fig. 9; again, the darkest bar(s) correspond to the best performing algorithm(s) and the light bar(s) to the least performing algorithm(s). Table VI summarizes average detection rates amongst groups of attacks.

In average, the three algorithms now feature similar robustness percentages: [START_REF] Solachidis | Circularly symmetric watermark embedding in 2-D DFT domain[END_REF] and [START_REF] Poljicak | Discrete Fourier transformbased watermarking method with an optimal implementation radius[END_REF] tie at 65.3%, and the proposed method reaches 63.2%. In the proposed method, the robustness to individual attacks is very similar to the scores obtained at default strength (see Fig. 6). Still, the robustness to filtering, shearing and random bending moderately improves with optimal strengths. In [START_REF] Poljicak | Discrete Fourier transformbased watermarking method with an optimal implementation radius[END_REF], the robustness is significantly improved for all kinds of attacks; it even slightly outperforms our method against JPEG, rotation, and the rotation & scaling. As a consequence, the introduction of a perceptual model into [START_REF] Poljicak | Discrete Fourier transformbased watermarking method with an optimal implementation radius[END_REF] would be very likely to significantly enhance its robustness instead of targeting a given PSNR.

Surprisingly in [START_REF] Solachidis | Circularly symmetric watermark embedding in 2-D DFT domain[END_REF], there is no overall loss in robustness despite the strong reduction in the embedding strength (α ⇤ =0.24 • α 0 in average), but on the contrary a drastic improvement. This can be explained by the fact that the normalized correlation value in [START_REF] Solachidis | Circularly symmetric watermark embedding in 2-D DFT domain[END_REF] is inversely proportional to the embedding strength. The detection threshold from [START_REF] Solachidis | Circularly symmetric watermark embedding in 2-D DFT domain[END_REF] being fixed, the FP rate is likely to rise.

Besides high detection capabilities, providing stable performances over a rather large range of embedding strengths is also essential. In the proposed technique, variations in the embedding strength do not significantly affect the robustness thanks to the efficiency of the perceptual model that accurately set the watermark strength with respect to its visibility threshold. Conversely in [START_REF] Solachidis | Circularly symmetric watermark embedding in 2-D DFT domain[END_REF] and [START_REF] Poljicak | Discrete Fourier transformbased watermarking method with an optimal implementation radius[END_REF], small variations in the embedding strength strongly affect their robustness, in other words, unstable performances.

May 6, 2014 DRAFT remain stable across experimented w i ; in [START_REF] Poljicak | Discrete Fourier transformbased watermarking method with an optimal implementation radius[END_REF], they are dramatically reduced (from 78% down to 5%) at low w i ;

surprisingly in [START_REF] Solachidis | Circularly symmetric watermark embedding in 2-D DFT domain[END_REF], they increase (from 52% up to 88%) with decreasing w i .

PSNR (dB) 2AFC det ect ion rat e This effect was already observed in Sec. VIII-A (see Fig. 9). We further investigated this unexpected behavior and studied the Receiver Operating Characteristics (ROCs) of [START_REF] Solachidis | Circularly symmetric watermark embedding in 2-D DFT domain[END_REF]. The ROC curves obtained for various strengths w • α 0 , with w ranging from 0.24 (the average w ⇤ across dataset D a in [START_REF] Solachidis | Circularly symmetric watermark embedding in 2-D DFT domain[END_REF]) to 0.14 (the minimum w ⇤ ). While the ROC curve remains near-optimal for w =0 .24, it quickly declines at lower strengths. At w =0 .14, the TP rate starts decreasing at high FP rates (10 -1 ) and drops to 0.1 for an FP rate of 10 -16 . Neither in the proposed technique nor in [START_REF] Poljicak | Discrete Fourier transformbased watermarking method with an optimal implementation radius[END_REF] was this observed. C. Content-dependency: some aspects Fig. 12 plots, for each algorithm, the average robustness to Stirmark attacks in each image belonging to dataset D a , at default strength (dashed line) and optimal strength (solid line). In addition, the bars in Fig. 12 plot the absolute difference in robustness between default and optimal strengths, i.e. the gain or loss in robustness when switching from α 0 to α ⇤ . It can be seen that the performances are consistent across all images in the proposed approach. Conversely in [START_REF] Solachidis | Circularly symmetric watermark embedding in 2-D DFT domain[END_REF] and [START_REF] Poljicak | Discrete Fourier transformbased watermarking method with an optimal implementation radius[END_REF], performances are much more inconsistent. This can be explained by the fact that, thanks to the accurate estimation of the visibility threshold provided by the proposed perceptual model, the watermark strength is properly adapted to the visual contents whilst this is not the case in [START_REF] Solachidis | Circularly symmetric watermark embedding in 2-D DFT domain[END_REF] and [START_REF] Poljicak | Discrete Fourier transformbased watermarking method with an optimal implementation radius[END_REF].

Robust ness (%) 0 20 40 60 80

• • • • • • • • • P roposed

IX. ROBUSTNESS TO PRINT &SCAN

In order to validate the proposed scheme in a more realistic scenario, the robustness to Print & Scan (P&S) of the proposed method was evaluated on seven 512 × 512 standard images (Baboon , Barbara, Boats, Fruits, Lena, Monarch and Peppers denoted as D i s , 1 ≤ i ≤ 7, in Table VII) and nine images of the D a database (denoted D j a , 1 ≤ j ≤ 9). As demonstrated in Sec. VII-C, default (α 0 ) and optimal (α ⇤ ) strengths are close; here, images were watermarked at α 0 . To ensure best reproducibility, each watermarked image was first printed five times at 600 ppi (with a 80% downscaling) on a Dell 2335dn laser printer. Printed images were then scanned at [START_REF] Mantiuk | HDR-VDP-2[END_REF] DRAFT 75, 150, 200 and 300 ppi with a Lexmark CX410de scanner, thus leading to 320 scanned images (images were placed approximately straight on the scanner glass). For each watermarked image and scanning resolution, Table VII lists the number of images (out of five) for which the watermark was detected. As can be seen, the proposed approach is not robust to very low scanning quality (75 ppi), but is robust to slightly higher resolutions 150 to 300 ppi. The detection fails only in text images wilk (see Fig. 3b) and dsc. In comparison, [START_REF] Poljicak | Discrete Fourier transformbased watermarking method with an optimal implementation radius[END_REF] is robust to a combination of print (600 ppi) and scan (150 ppi). The technique in [START_REF] Solanki | Print and Scan' Resilient Data Hiding in Images[END_REF] is not robust to scanning resolutions below 600 ppi, while the technique in [START_REF] Kang | Efficient General Print-Scanning Resilient Data Hiding Based on Uniform Log-Polar Mapping[END_REF] can withstand resolutions as low as 100 ppi.

An informal 2AFC subjective test was run on the 16 prints with 8 observers. The results show that the watermark remains under the visibility threshold (53.1% of detection).

X. CONCLUSION

This paper proposes a new watermarking method. The watermark, a square patch of coefficients, is embedded within the Fourier domain by substitution of both the magnitude (energy) and the phase (information). The watermark strength is perceptually optimized. The detection features both template matching and outlier detection, the latter being applied to the obtained correlation matrix. The decision is positive if at least one outlier is detected, and negative otherwise.

The proposed method was extensively compared to two competing algorithms from the literature. A subjective experiment was conducted in order to determine the perceptually optimized watermarking strengths (i.e. at the visibility threshold). The performances of both the proposed and the compared algorithms were evaluated twice: at default then perceptually optimal strengths. Experimental results showed that:

1) the proposed perceptual model accurately sets the watermark to its visibility threshold and is stable across all kinds of experimented contents;

2) the template matching accurately locates the watermark for slight geometrical distortions;

3) Grubbs' test for outlier performs very well both in terms of True Positives and False Positives;

May 6, 2014 DRAFT 4) the proposed method is robust against Print & Scan and shows state-of-the-art performances.

ACKNOWLEDGMENT

This work was partially supported by project OSEO A1202004 R, from which the author MU was funded.

35 Fig. 2 .

 352 Fig. 2. Correlation matrices ξ(u, v): how the dynamic range of the detection peaks and the surrounding noise varies with the watermark size M .

  BLIND WATERMARK DETECTION Let [ I sRGB be a supposedly watermarked image; it is assumed to have the same resolution and dimensions as the original image I sRGB , otherwise it is re-scaled. Let Y local ( b I) denote its locally normalized luminance map (see Section III-B). The proposed detection algorithm is blind, thus neither requires the original image I sRGB nor the original Fourier coefficients nor the CSF weighting coefficients. It is performed in two steps. At first, Template Matching (TM) is used to compute a 2D correlation matrix between the Fourier transform of Y local ( b I) -denoted c F Y -and the

Fig. 3 .

 3 Fig. 3. Examples of atypical images (cartoon and text): watermarked images.

Fig. 4 .

 4 Fig. 4. True Positive and False Positive detection rates against Stirmark attacks on Da: influence of the watermark size (M ) and α G

Fig. 5 .

 5 Fig. 5. H0 (original) versus H1 (marked): comparison with [15] and [16] on dataset D b . Solid lines represent the experimental histograms, and colored areas represent the fitted Gaussian models.

Fig. 7 Fig. 7 .

 77 Fig. 7. Proposed method: 2AFC analysis on image ED -dataset Da

Fig. 8 .

 8 Fig. 8. Distribution of perceptually optimized weights w * = α * /α 0 on Da.

Fig. 8

 8 Fig.8shows that algorithm[START_REF] Solachidis | Circularly symmetric watermark embedding in 2-D DFT domain[END_REF] overestimates the watermarking strength by factors ranging from 2 to 8: the embedded watermark is always visible at default strength. The other two algorithms perform much better in this regard: the distribution of their w ⇤ are nearly centered around the default position; misestimation factors remain

Fig. 9 .

 9 Fig. 9. Robustness to Stirmark attacks at visibility threshold (α * ): detection rates on dataset Da

Fig. 10 .

 10 Fig. 10. 2AFC visibility percentage versus PSNR for Da's image k21

Fig. 11 .

 11 Fig.11. Receiver Operating Characteristic curves of[START_REF] Solachidis | Circularly symmetric watermark embedding in 2-D DFT domain[END_REF] for various strengths.

Fig. 12 .

 12 Fig.[START_REF] O'ruanaidh | Rotation, scale and translation invariant digital image watermarking[END_REF]. Average robustness per image for dataset Da at default strength α 0 (dashed lines) and optimal strength α * (solid line). The bars plot the corresponding difference in robustness.

a

  Out of five printed & scanned copies.

  The viewing distance d is generally expressed as a multiplicative factor of the active image area's height, such that d = d • S y , where d is the normalized viewing distance.

	B. From pixel values to perceived contrast
	Contrast sensitivity models generally take physical luminance in cd.m -2 as an input; digital images on the other
	hand are gamma encoded for display purposes. The proposed model assumes that a typical monitor (γ =2 .2)
	is used for display; the standard RGB (sRGB) color space [26] will thus be used. Gamma expansion (Γ) is first
	applied to I sRGB in order to transform the original sRGB pixel values into linear RGB values: I RGB = Γ (I sRGB ). In
	a second step, I RGB is converted into the CIE XYZ color space, with D65 white point. Let I XYZ denote the obtained
	image; its component Y is proportional to the physical luminance and will thus be used for contrast computation
	and watermark embedding.
	May 6, 2014	DRAFT

TABLE II STIRMARK

 II BENCHMARK AT DEFAULT EMBEDDING STRENGTH: AV E R AG E D E T E C T I O N R AT E S a (%) AMONGST GROUPS OF ATTACKS ONDATASET Da. In each column, the bold value designates the best performing algorithm.

	Alg.	Overall	Rows & cols	removal	Filtering &	sharpening	Cropping	JPEG	Linear	Ratio	Rotation	Rotation &	scale	Scale	Shearing	Random bend
	Prop. 62.5 100 88.9 29.6 27.8 100 100 51.4 50.0 88.9 94.5 55.6
	[16] 53.1 77.8 71.1 33.3 51.9 59.3 76.4 45.8 48.6 64.8 40.7 44.4
	[15] 52.1 100 100 0.0 100 92.6 100 11.8 11.1 100 42.6 11.1

a

  Table III lists the values for w i , whose suitability was assessed by 4 expert observers who performed the proposed experiment prior to naive observers. In total, 243 images were included in the subjective experiment.

	May 6, 2014	DRAFT

TABLE III DATASET

 III Da: WATERMARK STRENGTH WEIGHTING COEFFICIENTS w i Algorithm Img. w 1 w 2 w 3 w 4 w 5 w 6 w 7 w 8 w 9

TABLE IV DATASET

 IV Da: PERCEPTUALLY OPTIMAL STRENGTH WEIGHTING COEFFICIENTS w *

			Algorithm ED k03 k07 k15 k20 k21 k23 dsc wilk
			Proposed 0.62 2.10 2.98 1.46 0.46 1.96 3.00 0.84 0.55
			[15] 0.32 0.48 0.30 0.20 0.18 0.14 0.22 0.14 0.16
			[16] 0.88 1.08 0.72 1.01 1.04 0.49 0.96 3.64 1.30
	Frequency	01234	P roposed
	Frequency	01234	[16]

TABLE VII ROBUSTNESS

 VII TO PRINT &SCAN: NUMBER OF DETECTED IMAGES a IN FIVE STANDARD 512 × 512 IMAGES AT SEVERAL SCANNING

RESOLUTIONS

Image D 1

[START_REF] Mantiuk | HDR-VDP-2[END_REF] DRAFT

http://www.irccyn.ec-nantes.fr/ ∼ autrusse/DFTWmking/ May 6,

DRAFT

B. Objective quality

Among the three tested algorithms, three different strategies are observed. In the proposed method, a perceptual model automatically adjusts the watermark strength at the visibility threshold. In [START_REF] Solachidis | Circularly symmetric watermark embedding in 2-D DFT domain[END_REF], the embedding strength is fixed and is thus independent of the image content. In [START_REF] Poljicak | Discrete Fourier transformbased watermarking method with an optimal implementation radius[END_REF], the embedding strength is progressively increased, so as to reach a target PSNR of 40 dB. Although it is widely accepted that the PSNR does not reflect image quality, it is still a recurrent quality metric when assessing the visibility of a watermark. . At a constant PSNR of ∼ 37.9 dB (see the vertical dimension line A in Fig. 10), the proposed watermark is invisible (2AFC percentage of 27.3%), while the watermark from [START_REF] Solachidis | Circularly symmetric watermark embedding in 2-D DFT domain[END_REF] is highly visible (2AFC percentage of 93.9%). Focusing now on data points close to the visibility threshold (see the horizontal dimension line B in Fig. 10), it appears that the obtained PSNRs significantly differ between [START_REF] Poljicak | Discrete Fourier transformbased watermarking method with an optimal implementation radius[END_REF] (43.2 db) and the proposed method (29.3 dB), while the watermark in [START_REF] Poljicak | Discrete Fourier transformbased watermarking method with an optimal implementation radius[END_REF] is slightly above the visibility threshold (81.8%) and the proposed watermark is slightly under the threshold (69.7%). In addition, the robustness of the proposed method is higher (62%) than the one of [START_REF] Poljicak | Discrete Fourier transformbased watermarking method with an optimal implementation radius[END_REF] (56%). Therefore, assessing the quality of experimented algorithms, from their PSNRs only, would conflict with the subjective results (ground truth), which in turn would lead to erroneous conclusions.

Concerning the robustness performances, one can see that experimented algorithms behave differently with varying w i . As it was already seen in section VIII-A, slight variations in embedding strength have a significant impact on the robustness in [START_REF] Poljicak | Discrete Fourier transformbased watermarking method with an optimal implementation radius[END_REF] and [START_REF] Solachidis | Circularly symmetric watermark embedding in 2-D DFT domain[END_REF], but have little impact on the proposed approach. In the proposed method, performances May 6, 2014 DRAFT


Ronan Hamon

Pierre Borgnat

Patrick Flandrin

Céline Robardet
email: celine.robardet@insa-lyon.fr

NONNEGATIVE MATRIX FACTORIZATION TO FIND FEATURES IN TEMPORAL NETWORKS

Keywords: nonnegative matrix factorization, temporal networks, Fourier analysis, dynamic graphs, multidimensional scaling

Temporal networks describe a large variety of systems having a temporal evolution. Characterization and visualization of their evolution are often an issue especially when the amount of data becomes huge. We propose here an approach based on the duality between graphs and signals. Temporal networks are represented at each time instant by a collection of signals, whose spectral analysis reveals connection between frequency features and structure of the network. We use nonnegative matrix factorization (NMF) to find these frequency features and track them over time. Transforming back these features into subgraphs reveals the underlying structures which form a decomposition of the temporal network.

INTRODUCTION

Many systems associated to networks, whether physical, biological or social, can be described by graphs [START_REF] Newman | Networks: An Introduction[END_REF]. The study of the structure of these graphs helps us to describe, understand and predict the behavior of systems. Often, these systems are not frozen but have a time evolution, e.g. associated connections which appear and disappear over time, or apparition or disappearance of nodes. We study here temporal networks where the edges are changing with times, keeping the same given set of nodes (possibly isolated).

Temporal networks have been introduced and studied in different fields under various names. Holmes and Saramäki [START_REF] Holme | Temporal networks[END_REF] proposed a review about temporal networks which includes discussion about the different types of temporal networks, the measures of the structure and the associated models. Casteigts et al. [START_REF] Casteigts | Time-varying graphs and dynamic networks[END_REF] also reviewed some works about dynamic networks with a graph theory approach. Evolving graphs were also considered as models for social networks [START_REF] Grindrod | Social networks: Evolving graphs with memory dependent edges[END_REF] or networks with information spreading [START_REF] Clementi | Information spreading in stationary markovian evolving graphs[END_REF]. As for visualization, Xu et al. [START_REF] Xu | Visualizing the temporal evolution of dynamic networks[END_REF] proposed a method based on multidimensional scaling (MDS). The temporal networks considered in this work consist in a series over time of snapshots of graphs, which are static at a given instant of time.

We propose a new approach based on the duality between graphs and signals to study temporal networks. Duality between graphs and signals consists in the transformation from (static) graphs to signals and inversely, in order to study either signals thanks to graph theory or graphs with signal processing. Campanharo et al. [START_REF] Campanharo | Duality between time series and networks[END_REF] proposed a This work is supported by the programs ARC 5 and ARC 6 of the région Rhône-Alpes method where the graphs represent a Markov chain coding the signals. Haraguchi et al. [START_REF] Haraguchi | Transformation from complex networks to time series using classical multidimensional scaling[END_REF] and later Shimada et al. [START_REF] Shimada | From networks to time series[END_REF] proposed a method to transform graphs into time series using classical multidimensional scaling (CMDS). We previously extended this approach in order to highlight relations between the structure of a graph and the associated frequency patterns, as well as to visually track temporal networks by following such patterns [START_REF] Hamon | Transformation de graphes dynamiques en signaux non stationnaires[END_REF]. This method has been exploited to build approximation of the graph over time with an application to the study of a bike sharing system in Lyon [START_REF] Hamon | Tracking of a dynamic graph using a signal theory approach : application to the study of a bike sharing system[END_REF] [START_REF] Hamon | Networks as signals, with an application to bike sharing system[END_REF]. The novelty here is to decompose a whole temporal network, first by obtaining frequency patterns by duality and then by using nonnegative matrix factorization (NMF) [START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF] to decompose these local spectra of the temporal network into features associated with weights evolving over time.

Section 2 recalls the duality mapping from graphs to signals. In Section 3 we extend the method to temporal networks using NMF and explain how the features are obtained. Section 4 gives an example on a controlled model of synthetic temporal network. A conclusion is given in Section 5.

FROM NETWORKS TO SIGNALS

Let G be a simple undirected and unweighted graph with n nodes. We note (Aij)i,j=1,..,n its adjacency matrix.

Principle

Shimada et al. [START_REF] Shimada | From networks to time series[END_REF] proposed a method to transform a graph into a collection of signals of n points indexed by the vertices of the graph by using classical multidimensional scaling (CMDS).

This transformation consists in applying CMDS to a matrix distance between vertices of a graph, noted ∆ = (δij)i,j=1,..,n and defined for two nodes i, j of G by

δij =    0 if i = j 1 if aij = 1 and i = j w > 1 if aij = 0 and i = j
Following [START_REF] Shimada | From networks to time series[END_REF], we choose w = 1.1.

Multidimensional scaling (MDS) [START_REF] Borg | Modern multidimensional scaling: Theory and applications[END_REF] is a set of mathematical techniques used to represent measurements of dissimilarity among pairs of objects as distances between points in a multidimensional space whose dimension is low. Classical MDS is a particular case of metric MDS where the dissimilarities are assumed to be Euclidean distances. The matrix X of coordinates in the low-dimensional, transformed space can be computed analytically. Starting with the distance matrix ∆, we first compute a double centering of the matrix whose terms are squared :

B = -1 2 J∆ (2) J with J = In -1 n 1n1 T n
where In is the identity matrix and 1n1 T n a n × n matrix of ones. The CMDS solution is given by X = Q+Λ + with Λ+ a diagonal matrix whose terms are the eigenvalues of the matrix B sorted in an increasing order and Q+ is the matrix of the corresponding eigenvectors. The resulting signals are the components (or columns) of the matrix X. The j-th signal is noted X (j) .

The inverse transformation from signals to graph is trivial: according to the principle of the CMDS, the distance matrix D between points in the Euclidean space is equal to the distance matrix ∆. The adjacency matrix A derives clearly from ∆.

Relabeling of a graph

Shimada et al. [START_REF] Shimada | From networks to time series[END_REF] showed that ring lattices are transformed into periodic signals: each component is a monochromatic oscillation whose frequency depends on the index of the component. More generally, it is relevant to use spectral analysis to link frequency features of signals with graph properties. Spectral analysis is nonetheless closely related to the indexation of signals and so to the labeling of the graph. Finding a good labeling of vertices is therefore essential to avoid sudden variations of signals. That means that it is necessary to have close labeling between neighbor vertices, which are defined closer in the distance matrix than unlinked vertices. This problem can be related to another graph labeling problem called cyclic bandwidth sum problem [START_REF] Hao | Cyclic bandwidth sum of graphs[END_REF]. We proposed in [START_REF] Ronan Hamon | Relabeling nodes according to the structure of the graph[END_REF] a heuristic that will be used throughout this paper but not explained further. 1. Procedure of the analysis of a temporal network. Description of each step is given in Section 3.

Analysis of signals

Let us consider a collection of K signals indexed by n vertices. We use the amplitude spectrum of the Fourier transform. From the magnitude of each frequency for each component k ∈ {1, . . . , K}, this amplitude reads as:

S(k, f ) = |F X (k) (f )| (1)
estimated, for positive frequencies, on F = n 2 + 1 bins, F being the Fourier transform.

Fig. 1 shows that specific graph properties, such as being regular in degrees (Watts-Strogatz model [START_REF] Newman | Networks: An Introduction[END_REF]) or being structured in communities, have specific frequency patterns. A random structure (Erdös-Rényi model [START_REF] Newman | Networks: An Introduction[END_REF]) is also visible in the spectrum of the collection of signals as the magnitudes are spread out over all the frequencies.

TEMPORAL NETWORK AND NMF

The duality between networks and signals defined in the previous section enables us to link frequency features with graph properties. We propose in this section to extend the transformation towards temporal network and use nonnegative matrix factorization (NMF) to catch these features and follow their level of activation over time. We finally exhibit a method to build subgraphs from these features in order to visualize the graph properties highlighted by NMF. Table 1 summarizes the steps of our procedure, and the elements we detailed in this section.

Transformation of temporal networks

We associate to a temporal network a dynamic graph defined as a succession of graphs. We note Gt the graph at time t for t = 0, . . . , T -1 where T is the duration of the considered evolution. Gt consists of all the nodes and edges activated at time t. We then define the transformation of such a temporal network into signals as the transformation of Gt for all t = 0, . . . , T -1. We note Xt for t = 0, . . . , T -1 the collection of signals obtained at time t.

We assume that the number of components K and so the number of frequencies F is the same at each time step. This assumption can be easily achieved by fixing the number of components as the maximal number of nodes and zero-padding the missing components. We can now define for each component k

= 1, . . . , K St(k, f ) = 2|F X (k) t (f )| for f = 0, . . . , F -1 as previously.
St represents T matrices of dimension K × F . For t fixed, it is possible to reshape the matrix St to a vector of size (KF ) by successively adding end-to-end the columns of the matrix St. For all t = 0, . . . , T -1, we can form the matrix V of dimension (KF )×T where V (t) t-th column of V , is equal to the vectorized matrix St.

Identification of features

Nonnegative matrix factorization (NMF) [START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF] is used to decompose St(k, f ) into features and levels of activation over time. NMF can be written as the following problem: given a nonnegative matrix V of dimension A × T , find a factorization V ≈ W H where W and H are two nonnegative matrices of dimensions respectively A × Q and Q × T , with Q the number of features (usually small). W is described as the features matrix and H the corresponding levels of activation over time. Formally the NMF problem almost comes down to solve:

min W ,H A-1 a=0 T -1 t=0 d([V ]at|[W H]at) (2)
subject to W and H nonnegative matrices. d(x|y) is a scalar cost function.

Févotte et al. [START_REF] Févotte | Algorithms for nonnegative matrix factorization with the β-divergence[END_REF] proposed an algorithm to find a solution of the NMF where the cost function is the β-divergence, a parametrized function with a simple parameter β which encompasses the Euclidean distances (β = 2), the generalized Kullback-Leibler divergence (β = 1) and the Itakura-Saito divergence (β = 0) as special cases. Varying β and the number of features Q has not been here the subject of specific studies. Here β is chosen equal to 1 because of the spectral origin of V . We will adapt Q to the expected number of features in the example.

We use NMF to find patterns in spectra of the collections of signals obtained from the transformation of the temporal network. At each time t, the t-th column of V (t) represents the vectorized representation of St(k, f ). V of dimension A × T with A = KF is taken as input of the algorithm.

Reconstruction of signals from features

The reconstruction of signals from features is tricky considering that the feature represents only the magnitudes of the spectrum. To obtain a consistent collection of signals, it is necessary to provide some additional information about the phase. The strategy used here is based on the assumption that at some specific time step, the considered feature is the only one activated in the network or is dominant. That means that the level of activation of the considered feature is high while levels of activation of other features are low. That leads to select the time step which maximizes the following criterion: for

q = 0, • • • , Q and t = 0, • • • , T -1, tq = arg max t=0,••• ,T -1 (Hqt - p∈{0,••• ,Q-1}\{q} Hpt) (3)
Once properly identified, the phase φ of the signals representing the graph at step tq is consistent with the magnitudes and so the most activated feature.

Reconstruction of graph from features

The transformation from a feature to a graph is complicated because the reconstructed signals are not directly derived from a transformation from graph to signals but have been modified. The distances between points are neither equal to 1 nor to w but have a distribution whose width can be variable. Under the assumption that if two points have a small distance, then the proximity between these points in the network is high, a reconstruction can be set up such that only the links associated to the smallest distances are retained. A good threshold would be the one which maximizes the separation between e ∈ E p e / ∈ E p e ∈ Et-1 1 -10 -3 0.9 e / ∈ Et-1 0.1 10 -3

Table 2. Probability to have an active link in E according to the presence of this link in Et-1, the active edges in the temporal network at time t, and Ep the set of edges in the prescribed graph.

the distances which denotes linked pairs of nodes, and the distances which denotes unlinked pairs of nodes. Nonetheless it might be in practice very tough to select such a threshold, except manually. We propose here to use the density of the graph at the best time step as described previously, to select the number of edges. In order to highlight visually the structure, we select a higher number of links (1.5 times the density of the original graph).

EXAMPLE

This section illustrates the previous sections on an example of a simulated temporal network which contains specific structures. The aim of this example is to show that the method is able to separate features corresponding to these structures, even in presence of noise or features coupled between them. The analysis of this network follows the steps of the procedure described in Table 1. The simulated temporal network is generated using the following procedure: starting from a graph without links, the algorithm adds or removes links at each time instant according to a probability depending on a prescribed graph. This probability is given in Table 2 and depends on whether an edge is in the prescribed graph and whether it exists in the graph at the previous instant. The algorithm runs during 30 time steps for a prescribed graph. Four different graphs with a constant number of nodes equal to n = 49 have been used and drawn from the following models: The graph at each time instant is transformed into a collection of signals and frequency analyses are performed to obtain the matrix V with the vectorized spectra at each time at each column. The NMF algorithm is applied to this matrix using β = 1, the number of features is Q = 3 and the number of iterations is 3000. Fig. 3 displays on the left the obtained features reshaped in a matrix form to show magnitude of each frequency for each component. Connections between frequency patterns highlighted in Fig. 1 and these features can be made. Feature 3 looks like the pattern of a Watts-Strogatz model of degree k = 2 (Fig. 1a) with high values on the diagonal. Feature 2 is close to the pattern of a graph with communities (Fig. 1b) where the magnitudes of low frequencies for the first components are very high. As for feature 1, the magnitudes spread out over all frequencies are related to those found for the Erdös-Rényi model (Fig. 1c). These three structures correspond to the three models used to build the temporal network.

Levels of activation over time of features are displayed on the right of Fig. 3, as well as the four intervals corresponding to the four prescribed models used to build the temporal network. It reveals that the features are activated at the expected time steps. Feature 3 associated to a regular structure is fully activated when the model (CYC) is active, and activated with less intensity when the prescribed model is (BAR). For others time steps the feature is close to 0. As to feature 2, linked with structure in communities, it is activated when models (COM) and (BAR) are actives, composed of cliques. Finally feature 1, corresponding to a random structure, is mainly activated when model (ER) is active and for model (COM), denoting the random structure inside communities which are not totally built as cliques. It is also activated at the beginning of each interval, describing the random transition between two prescribed graphs. Construction of subgraphs from features enables us to visualize the structure of the temporal network when the features are activated. Fig. 4 displays these subgraphs corresponding to the features. The layout takes into account the links between nodes and enables us to have an overview of the structure of these subgraphs. These structures correspond to the frequency features by similarity with illustrations given in Fig. 1. They also correspond to the models we used to build the graph. NMF is hence able to separate, in a temporal network with three different structures coupled between them, these subgraphs. This gives indication about their time instants of activation. The shape of the distributions of weights used to build the subgraphs is displayed in insert in order to highlight the difficulty to choose a good threshold.

CONCLUSION

We have proposed a novel method to track the structure of temporal networks over time using NMF and duality between graphs and signals. At each time, the temporal network is represented by a graph which is transformed into a collection of signals. NMF is used to discover features in the amplitude spectra of these signals; these features, whose levels of activation vary over time, represent a specific structure of the graph. The effectiveness of the method was demonstrated on a simulated temporal networks containing three types of structures.

These results provide insights in the characterization of temporal networks, but also call for further studies, in particular concerning the parameters in the NMF in this context. The extension to a continuous model for temporal networks, where some continuity in time of the transformation into signals, is also worth considering.

Fig. 1 .

 1 Fig. 1. (left) Snapshot of the network. (right) Frequency analysis of the resulting signals. (a) The regularity of the structure of is apparent as high values on the diagonal of S(k, f ). (b) Structure in communities implies very high magnitudes at low frequencies for the first components. (c) Random structure leads to a spread of magnitudes over all frequencies.

#

  Definition of the steps 1 Transformation of Gt into a collection of signals Xt 2 Frequency analysis of signals: V t = |F Xt(f )| 3 NMF: V ≈ W H 4 Identification of features W k and levels of activation H k 5 Transformation into a spectrum Ŷ = W k .e iφ 6 Inverse Fourier transformation: Y = F -1 Ŷ 7 Transformation from the signals Y to a graph H Table

4. 1 .Fig. 2 .

 12 Fig. 2. Four snapshots of the models used for the generation of the temporal network.

1 .

 1 [0, 30] -Erdös-Rényi model p = 0.5 (ER) 2. [30, 60] -Two cliques of variable size (COM) 3. [60, 90] -Cycle graph (CYC) 4. [90, 120] -Barbell model (BAR): two cliques of size 15 connected with a path of length 19 These models represent three different structures which can be found in networks: random linkage of nodes in model 1, structure in communities in model 2 and regularity in model 3. Frequency patterns related to signals describing graphs with these structures have been highlighted in Fig. 1. The model 4 is a mixture between structure in communities and regularity.

4. 2 .

 2 Steps 1-2-3-4: Transformation into signals and identification of features with NMF

Fig. 3 .

 3 Fig. 3. (left) Features found using NMF, reshaped in a matrix form to show magnitude of each frequency for each component. Connections can be made with magnitudes found in Fig. 1: feature 1 with (c), feature 2 with (b) and feature 3 with (a). (right) Levels of activation over time are displayed as well as the four intervals corresponding to the four prescribed models used to build the temporal network.

4. 3 .Fig. 4 .

 34 Steps 5-6-7: Transformation of features into static graphs Graphs built from features. (insert) Shape of the distribution of distances. The vertical line indicates the retained threshold. Random structure for feature 1, structure in communities for feature 2 and regular structure for feature 3 are visible.