Adam Cheminet
email: adam.cheminet@cea.fr

Yasar Ostovan

Valentina Valori

Christophe Cuvier

Fançois Daviaud

Paul Debue

Bérengère Dubrulle

Jean-Marc Foucaut

Jean-Philippe Laval

Optimization of regularized B-spline smoothing for turbulent Lagrangian trajectories

Keywords: Turbulence, Lagrangian measurements, signal processing, denoising techniques, regularized B-splines

The denoising of Lagrangian trajectories based on regularized B-spline is investigated. The aim is to find systematic criteria for optimization of algorithms used in 4D-PTV in order to optimize the quality of 4D-PTV measurements of turbulent flows as well as high-order of turbulence statistics. We introduce and adapt to this context two innovative tuning strategies which are commonly used in the Tikhonov regularization of inverse problems based on L-curve shape and Normalized Cumulative Periodogram (NCP). The corresponding strategies are tested on synthetic Lagrangian trajectories computed from Direct Numerical Simulation with additional white Gaussian noise. Error-based quantities like Signal-to-Noise Ratio as well as statistical Lagrangian quantities are investigated to compare the different strategies. We then apply the algorithm to experimental data from a 4D-PTV Lagrangian measurements in a turbulent Von Kármán flow. We show the ability of those strategies to optimize the quality of the signal compared to conventional methods. Moreover, the strategies are more adaptable to real experimental noise.

Introduction

Turbulent flows are inherently time dependent, multi-scale and three-dimensional. Corresponding measurements need to resolve all the timescales and lengthscales of the flow from the largest energy injection scales down to the smallest dissipative scales. State-of-the art methods, termed 4D-Particle Tracking Velocimetry (4D-PTV) [START_REF] Schanz | Shake-The-Box: Lagrangian particle tracking at high particle image densities[END_REF], strive to achieve this. These optical techniques are based on the displacement of passive seeded tracers in the flow, lighted by a laser volume and imaged in successive frames separated by a small time interval. Reconstruction of the flow field can then be obtained through the displacement of each individual tracer in-between two consecutive frames. Historically, 3D-Particle Tracking Velocimetry (3D-PTV) first emerged in the 1990's [START_REF] Maas | Particle tracking velocimetry in three-dimensional flows. part 1. photogrammetric determination of particle coordinates[END_REF], [START_REF] Malik | Particle tracking velocimetry in three-dimensional flows[END_REF]) but could only give access to very low density Lagrangian flow fields by reconstructing through triangulation then by tracking individually each particle. In 2006, [START_REF] Elsinga | Tomographic particle image velocimetry[END_REF]) introduced tomographic reconstruction followed by cross-correlation on a voxel based representation of the volume to directly obtain a dense Eulerian flow field. However those methods were still impaired by the impact of ghost particles at relatively high particle densities. To overcome these obstacles, and in the case of time-resolved measurements, [START_REF] Schanz | Shake-The-Box: Lagrangian particle tracking at high particle image densities[END_REF] showed that combining tomographic reconstruction and the tracking of reconstructed particles could drastically increase the accessible seeding density thus improving the overall spatial resolution. This results in a hybrid Lagrangian/Eulerian measurement where the Eulerian flow is accessed by interpolation of the Lagrangian flow field. Numerous hybrid Lagrangian/Eulerian techniques [START_REF] Yang | Lagrangian Particle Image Velocimetry[END_REF], [START_REF] Lasinger | 3D Fluid Flow Estimation with Integrated Particle Reconstruction[END_REF], [START_REF] Cornic | Double-frame tomographic PTV at high seeding densities[END_REF]) are emerging nowadays aiming to achieve higher seeding density thus enabling high spatial resolution of the measurement.

In practice, the 4D-PTV measurement consists in noisy particle trajectories. One can say that noise is inherent to any measurement system. The sources of noise in 4D-PTV are multiple. They can be internal to the measurement system, such as the CCD camera background noise, the pixelation of the cameras, the pixel sensitivity, laser illumination temporal and spatial inhomogeneity as well as external like vibrations. They can also come from experimental conditions such as the camera setup, imaging conditions and seeding density. This noise creates uncertainties which will propagate non-linearly through the measurement chain, from the calibration and self-calibration method to the method of 3D particle reconstruction (Bhattacharya et al. 2006). This uncertainty in the particle position is detrimental to the determination of the particle velocity and acceleration using differentiation [START_REF] Feng | Errors in particle tracking velocimetry with high-speed cameras[END_REF] as well as any direct computation of turbulent quantities based on higher order derivatives. This issue is well known in experimental studies of Lagrangian turbulence and is not restricted to 4D-PTV. A thorough review of Lagrangian experimental methods can be found in [START_REF] Bourgoin | Lagrangian Methods in Experimental Fluid Mechanics[END_REF]. Statistically speaking, this noise is often seen and presumed to be white and uncorrelated with the particle trajectory [START_REF] Machicoane | Estimating two-point statistics from derivatives of a signal containing noise: Application to auto-correlation functions of turbulent Lagrangian tracks[END_REF]) and will mostly impact the high frequencies of the signal, such as in Acoustic Doppler Velocimeter (García et al. 2016). Strategies to evade the noise issue exist and depend on what the experimenter is trying to measure. For instance, [START_REF] Machicoane | Estimating two-point statistics from derivatives of a signal containing noise: Application to auto-correlation functions of turbulent Lagrangian tracks[END_REF] designed a method to directly access statistical properties such as moments and correlation functions of signal derivatives, without using filtering schemes. However, in the case of 4D-PTV, the aim is to access both the instantaneous Lagrangian and Eulerian velocity field at any given measurement snapshot. Though this is still an ongoing research and techniques are rapidly evolving, state-of-the art methods compute the Eulerian flow field directly from a spatial interpolation and data assimilation of the Lagrangian flow field [START_REF] Gesemann | From Noisy Particle Tracks to Velocity, Acceleration and Pressure Fields using B-splines and Penalties[END_REF], [START_REF] Schneiders | Dense velocity reconstruction from tomographic PTV with material derivatives[END_REF]. The interpolation process is done at every measurement snapshot separately using the particles ' instantaneous positions, velocities and accelerations. Any measurement noise from the temporal evolution of Lagrangian field will pollute and reduce the spatial and time resolution of the Eulerian flow field. Removing this temporal noise as much as possible in the whole Lagrangian data (positions, velocity & acceleration) becomes necessary.

A classic approach consists in oversampling the particle trajectory so as to distinguish high fre-quency noise from the turbulent signal using finite-impulse response low-pass filters as pointed out in [START_REF] Lawson | Bias in particle tracking acceleration measurement[END_REF]). An ideal filter has two essential characteristics : a cutoff frequency which can be tuned depending on a desired frequency, and a strong negative cutoff slope that removes the noise while keeping the real signal as much as possible and that does not add any non-physical high frequencies which might lead to biases in statistics of higher-order derivatives. To this effect, several filters exist and have already been used in the literature. Filtering can be done through Least-square Polynomial fitting as in [START_REF] Voth | Measurement of particle accelerations in fully developed turbulence[END_REF], [START_REF] Del Castello | Lagrangian Acceleration of Passive Tracers in Statistically Steady Rotating Turbulence[END_REF], or fitting a moving cubic spline [START_REF] Lüthi | Lagrangian measurement of vorticity dynamics in turbulent flow[END_REF]. Those filters enter the category of Savitsky-Golay filters. A very common filter is the Gaussian filter which consists in doing a convolution of the signal with a Gaussian kernel [START_REF] Mordant | Experimental lagrangian acceleration probability density function measurement[END_REF], [START_REF] Ouellette | A quantitative study of three-dimensional Lagrangian particle tracking algorithms[END_REF], [START_REF] Xu | Tracking Lagrangian trajectories in position-velocity space[END_REF], [START_REF] Biferale | Lagrangian structure functions in turbulence: A quantitative comparison between experiment and direct numerical simulation[END_REF], [START_REF] Stelzenmuller | A Lagrangian study of inhomogeneous turbulence[END_REF] allowing for an easy computation of the filtered signal and its derivatives. The Gaussian filter has a stronger negative cutoff slope than Savitsky-Golay filters, thus eliminating high frequencies more efficiently [START_REF] Stelzenmuller | A Lagrangian study of inhomogeneous turbulence[END_REF], [START_REF] Schafer | On the frequency-domain properties of Savitzky-Golay filters[END_REF]. However, border issues are still present when using the convolution of kernel, limiting the use of such filters to relative long trajectories thus losing some valuable information. Finally, [START_REF] Berg | Experimental investigation of Lagrangian structure functions in turbulence[END_REF] chose a binomial filtering instead of Gaussian filtering, without clearly stating the benefit of such a choice. Whatever the filter shape, choosing the proper time scale and frequency cutoff is a key question in signal processing. Indeed, this choice will decide the trade-off between how much noise is removed and how much actual physical signal is kept. Too much smoothing will remove high frequency turbulent events but not enough smoothing will create errors due to noise. This question is either overlooked or studied rigorously. Several methods have been used but hardly ever following the same guidelines. The ideal filter length is often looked for in a range where a small difference in the filter scale leaves no impact on high-order statistical quantities such as the acceleration variance [START_REF] Lawson | Bias in particle tracking acceleration measurement[END_REF]). However, this range is not always observed [START_REF] Berg | Experimental investigation of Lagrangian structure functions in turbulence[END_REF]). On the contrary, two distinct behaviors towards the filter lengthscale are sometimes observed depending on the statistical quantity under study. [START_REF] Berg | Experimental investigation of Lagrangian structure functions in turbulence[END_REF] focuses on the acceleration PDF while [START_REF] Stelzenmuller | A Lagrangian study of inhomogeneous turbulence[END_REF] focuses on the wall normal acceleration variance. Each behavior corresponds to each extreme filtering situation : keeping too much noise or removing too much signal. The optimal filtering is thus found in-between those two ranges of filter lengthscale though [START_REF] Berg | Experimental investigation of Lagrangian structure functions in turbulence[END_REF] admit that their choice is somehow arbitrary with no rigorous technique to find the optimal filtering. In a context of evolving 4D-PTV techniques, a simple filtering algorithm for Lagrangian trajectories (TrackFit) was designed by [START_REF] Gesemann | From Noisy Particle Tracks to Velocity, Acceleration and Pressure Fields using B-splines and Penalties[END_REF]) using regularized B-Spline and the minimization of a weighted functional. Such minimization is performed via a tuning parameter which controls the level of smoothness of the solution. This parameter is directly linked to a cutoff frequency. The best cutoff frequency is chosen to be the frequency at which the Signalto-Noise ratio is equal to one. In TrackFit, the tuning is performed by selecting the frequency laying at the crossroads of the signal (i.e. position) and noise spectrum, either visually or by a corner detection [START_REF] Gesemann | From Noisy Particle Tracks to Velocity, Acceleration and Pressure Fields using B-splines and Penalties[END_REF]. [START_REF] Lawson | Bias in particle tracking acceleration measurement[END_REF] showed that this simple procedure could lead to significant gains in the accuracy of acceleration statistics compared to classical Gaussian filtering when combined with an additional independent measurement. In the quest for the optimal signal and in the context of detection of extreme and complicated flow events, one needs to optimize the quality of the experimental data with respect to high order statistics as well as instantaneous quantities. Therefore, the following questions come up : is this frequency the optimal filtering frequency ? Are the proposed strategies to find this frequency adapted to experimental conditions ? Are there other systematic and rigorous ways of finding this frequency ? If so, how do those methods compare to one another ?

The goal of this paper is to provide a tentative answer to these questions by investigating the performances of two alternative strategies, namely L-curve and NCP, and by comparing their performances with respect to several performance criteria : Error-based quantities, statistical Lagrangian quantities such as acceleration PDFs, velocity spectra, geometrical properties and velocity structure functions. We start by describing the mathematical background of the filter, as well as different tuning strategies. Then, we use Direct Numerical Simulation (DNS) of Homogeneous Isotropic Turbulence (HIT) seeded with convected Lagrangian tracers to perform tests on the filtering strategies with respect to the addition of a well-controlled Gaussian noise. Finally, we apply the algorithms to experimental data from a 4D-PTV Lagrangian measurements in a turbulent Von Kármán flow [START_REF] Ostovan | 4D Particle Tracking Velocimetry measurements in a Von Karman turbulence experiment[END_REF]) and compare them via several performance criteria.

Smoothing method and tuning strategies

Let X P (t i ) = (X(t i ), Y (t i ), Z(t i ))
be the coordinates of the measured P particle track at N discrete time instants such as t i = t 1 + (i -1)∆t. The measurement Nyquist frequency is defined as f N = 1 2∆t . For each coordinate, f (t i ), a Lagrangian trajectory smoothing technique is required to damper the noise level. In the sequel, we focus on techniques based on B-Spline projection in the "least square sense".

B-Spline curve

An n th degree B-spline curve S is a piecewise polynomial function made of B-spline basis functions of order n, b i,n (t), which are defined on a sequence of control points or knots α i such as

α 0 ≤ α 1 ≤ α 2 ≤ ... ≤ α m (1)
The B-spline basis functions are defined recursively using the De Boor formula [START_REF] De | A Practical Guide to Splines[END_REF] :

b i,0 (t) = 1 if α i ≤ t ≤ α i+1 0 otherwise b i,n (t) = t -α i α i+n -α i b i,n-1 (t) + α i+n+1 -t α i+n+1 -α i+1 b i+1,n-1 (t) (2)
Overall, the B-spline curve can be written as :

S(t) = m-n i=1 c i .b i,n (t) (3)
where c i are weighting coefficients.

Regularized B-spline

To fit a B-spline curve S(t) to a noisy signal f (t j ) j∈{1,N } and find the optimum control coefficient c i which reduces the noise, one needs to build a distance functional that has to be minimized :

F dist (c) = N j=1 f (t j ) - m-n i=1 c i b i,n (t j ) 2 (4)
The smoothness of the B-spline curve is controlled by the number and location α i of the knots with respect to the data points (t j ) as well as the polynomial order (n). A wide variety of methods are to be found in the literature for an optimal knot arrangement for a given smoothness criterion.

A different projection strategy (O'Sullivan 1986) & (O'Sullivan 1988) uses a large number of knots compared to the number of data points and has an additional regularization term F regu (c) in the functional which will penalize higher-order curve derivatives to avoid overfitting. [START_REF] Eilers | Flexible smoothing with B-splines and penalties[END_REF] use k th order finite difference of the adjacent B-Splines coefficients (∆ k c i ) as regularization. Strictly speaking, the regularization is not on the k th derivative of the B-spline fitted curve but on the k th derivative of its control coefficients (c i ). However, regularizing the control coefficient derivatives will penalize the curve derivatives which is the aim of the smoothing technique. In the case of the second derivative, [START_REF] Eilers | Flexible smoothing with B-splines and penalties[END_REF] showed the strong link between the two implementations.

The TrackFit method [START_REF] Gesemann | From Noisy Particle Tracks to Velocity, Acceleration and Pressure Fields using B-splines and Penalties[END_REF]) is very similar to [START_REF] Eilers | Flexible smoothing with B-splines and penalties[END_REF]. It uses third order B-spline functions from an analytical description, placed at and centered on every data points. One can say that they are used like Radial Basis Functions (RBF). This is equivalent to using third order (n = 3) B-spline curve defined as equation (2) and using the data points as the control knots of the B-spline curve (t i = α i , m = N ).

In TrackFit, penalization is performed on the particle acceleration derivative also called jerk. Since B-spline are third order, this is achieved using ∆ 3 c i as regularization term. The differentiating scheme used here is a first order forward finite difference. The regularization functional is :

F regu (c) = N -3 i=1 ∆ 3 c i 2 (5)
Smoothing is then achieved through minimization with respect to c of the following functional :

F (c) = F dist (c) + λF regu (c), ( 6
)
where λ is a given smoothing parameter, weighting the strength of the regularization on the jerk compared to the error on the data which needs to be tuned depending on the Signal-tonoise Ratio (SNR) and the expected cutoff frequency f c . Using generalized normal equations [START_REF] Eilers | Flexible smoothing with B-splines and penalties[END_REF], minimization is thus achieved when the functional gradient ∇F is set to 0, which gives :

B t y = B t B + λD t k D k c (7)
with B ij = b j,n (t i ) and D k , the k th order differentiation matrix (in our case, k = 3) and y i = f (t i ). To avoid border issues as much as possible, we put three additional equally spaced knots at the two extremities of the time series resulting in a knots series (t i ) i∈{-2,N +3} . B and D matrices are therefore a N by N + 6 -n matrix and a N -n + 3 by N + 6 -n matrix respectively, with n = 3 in our case.

Solving eq. 7, one then gets the c i as a function of λ, thereby getting the smoothed coordinate via eq. 3. In this procedure, we see that the choice of λ determines the smoothing properties. In Section 2.4, we will discuss the TrackFit choice, as well as two alternative choices based on properties of Tikhonov regularization. In the next section, it will be worth investigating the link between the smoothing parameter λ and the signal cutoff frequency f c .

Smoothing parameter and cutoff frequency

The Lagrangian smoothing acts as a low-pass filter. Like any such filter, it has a transfer function, H(f ) defined as :

H(f ) = F (f ) F (f ) (8)
where F (f ) is the Fourier transform of the unfiltered signal and F (f ) is the Fourier transform of the filtered signal. Knowing H(f ) allows to define a cutoff frequency f c . Since H(f ) depends on λ, so does f c . We can provide both numerical and analytical estimates of such dependence, by applying the smoothing technique to the case of a wave function for varying frequencies. Assuming that the filter response to a sinusoidal signal f (t) = Asin(kt) of period T is a sinusoidal signal f (t) = A sin(kt), one can model and write the distance F dist and regularization F regu functionals of the filter optimization problem.

F dist = 1 T T 0 (f (t) -f (t)) 2 dt = 1 2 (1 -) 2 A 2 F regu = 1 T T 0 ( f (t)∆t 3 ) 2 dt = 1 2 2 A 2 k 6 ∆t 6 (9)
Though the regularization is done on the derivative of the adjacent B-spline coefficients, and not the derivative of the filtered signal, a first rough approximation allows us to equal the two.

The filter full functional F ( ) can be now expressed as a function of , the damping parameter.

Setting the functional gradient to zero ( dF d = 0) allows us to obtain the analytical transfer function:

H(f ) = 1 1 + ( f fc ) 6 , f c = 1 πλ 1/6 f N (10)
f c can be defined as the frequency cutoff of the filter and provides a link between f c and λ.

In signal processing, it is also common to define the -3 dB frequency cutoff (f -3 dB ) as the effective frequency cutoff of the filter. f -3 dB can be obtained by setting H(f ) = 1/ √ 2 in eq. 10, leading to :

f -3dB = ( √ 2 -1) 1/6 f c (11)
To check such formula, the filter transfer function can be estimated numerically. To do so, we filtered a noiseless signal made of a sinusoidal wave of varying frequencies and for different regularization parameters. Fig. 1a shows the transfer function for different smoothing parameters. Fig. 1a clearly depicts the -6 slope behavior of the low-pass filter. Fig. 1b shows the link between the analytical description of f c given in eq. 11 and the measured -3dB frequency cutoff. It shows that eq. 11 is reliable and offers a good estimate of the relationship between the regularization coefficient and the cutoff frequency.

Two main characteristics of the filter can be inferred from this analysis. First, a clear relation between tuning parameter and cutoff frequency allows us to accurately tune the smoothing intensity so as to remove all the frequencies above a given frequency. Second, the filter exhibits a strong -6 slope suitable to efficiently denoise the signal. It is thus reasonable to think that the best cutoff frequency should be equal to or in the close vicinity of the frequency at which the Signal-to-Noise Ratio (SNR) is equal to 1. Depending on the definition of SNR, this can also be seen as the frequency where the low-frequency physical signal spectrum intersects with the noise spectrum [START_REF] Gesemann | From Noisy Particle Tracks to Velocity, Acceleration and Pressure Fields using B-splines and Penalties[END_REF]). This frequency is termed f SN R=1 thereafter.

2.4. Frequency cutoff finding strategies 2.4.1. Spectrum shape derived methods : To find the f SN R=1 , it is worth looking at the shape of the trajectories position spectrum. Indeed, as mentioned in [START_REF] Gesemann | From Noisy Particle Tracks to Velocity, Acceleration and Pressure Fields using B-splines and Penalties[END_REF], this spectrum can be divided in two parts : a low-frequency signal-dominated part and a highfrequency noise-dominated part. For particle trajectories, the signal spectrum usually exhibits a negative power law. In contrast, the noise spectrum is usually assumed to be flat (white noise). Therefore the f SN R=1 can be found in the vicinity of the corner or kink of the spectrum.

Two methods can be inferred : the spectrum corner can be either found by computing the maximum curvature of the spectrum curve in logarithmic scaling or by finding the intersection of two straight lines which locally approximate the two parts of the spectrum on both sides of its corner. Those two strategies lead to regularization parameters respectively termed λ c and λ d .

To compute λ c , the raw signal spectrum is first smoothed, interpolated on a finer grid and then its curvature is computed. The frequency yielding the maximum curvature is chosen as the cutoff frequency leading to a regularization parameter using eq. 10 as suggested in [START_REF] Gesemann | From Noisy Particle Tracks to Velocity, Acceleration and Pressure Fields using B-splines and Penalties[END_REF]). To compute λ d , the raw spectrum is fitted by two straight lines using least-square fitting. The cutoff frequency is defined as the frequency where both lines cross (noted f cross ), yielding a regularization parameter from eq.10. This criterion is very sensitive to the spectrum shape which does not always have a constant power law (see Fig. 12). The length of each segment is also user-related, thus making this criterion dependent on the insight of the experimenter.

2.4.2. Methods deriving from Tikhonov regularization strategies : It is interesting to realize the similarities between our regularization problem and Tikhonov regularization [START_REF] Tikhonov | Solutions of Ill-Posed Problems[END_REF]. In Tikhonov regularization, solutions of ill-posed problems are searched by minimizing the following functional :

y -Kx 2 2 + λ x 2 2 (12)
In this problem, properly choosing the regularization parameter is also crucial. Several techniques were developed to circumvent this issue. In this article, we will focus on two of them, namely the L-curve criterion [START_REF] Hansen | The use of the L-curve in the regularization of discrete ill-posed problems[END_REF]) and the Normalized Cumulative Periodogram (NCP) criterion [START_REF] Hansen | Exploiting residual information in the parameter choice for discrete ill-posed problems[END_REF]). The L-curve criterion uses the characteristic L-shape of the parametric curve

L(λ) = (log( y -Kx λ 2 ), log( x λ 2 )) where x λ = arg min x y -Kx 2 2 + λ x 2 2 .
The optimal regularization is searched in the vicinity of its corner (Fig. 2a). To locate the curve corner, we can use the L-curve curvature θ(λ), resulting in a criteria for the choice of λ such as :

λ L = arg max λ {|θ(λ)|} (13)
The NCP criterion uses the spectral behavior of the residual vector of the filter r λ = x f iltered(λ) -x data . Given the vector r λ (j) j∈[0,N ] , and its power spectrum

P λ (k) k∈[0,N/2] , its NCP is computed as : NCP r λ (k) = k n=0 |P λ (n)| N/2 n=0 |P λ (n)| (14)
NCP r λ visualizes the spectral behavior of r λ . If a signal is made of white noise, its NCP should be close to the identity line (see Fig. 2b). A signal dominated by low frequency (respectively high frequency) has its NCP above (and respectively below) the identity line. Ideally, the filter residual should correspond to the noise. Its NCP should therefore be the NCP of a noisy signal. Thus, the point is to choose the regularization parameter that gives an NCP of the residual vector that is as close as possible to an estimated measurement noise (NCP N oise ). A criterion based on the principles of a Kolmogorov-Smirnov test can be inferred such as :

λ N CP = arg min λ max k |NCP r λ (k) -NCP N oise (k)| (15)
We will call λ N CP W N the criterion based on the assumption that the noise is a white noise, thus where NCP N oise (k) is the x = y line.

However in our case, the TrackFit filter is a low pass-filter, which leaves unaltered low frequencies of the signal and damps the high frequencies above its cutoff frequency which depends on λ, namely, f c (λ). Therefore it is logical to compare the spectral behavior of r λ to the spectrum of the high-pass filtered estimated noise using the TrackFit transfer function. This high-pass filtered estimated noise can be seen as the best estimated residual TrackFit can produce for each λ. If E N oise is the spectrum of the estimated noise then the spectrum of the best estimated residual is

E EN,λ = (1 -H λ (f )) E N
λ N CP EN = arg min λ max k |NCP r λ (k) -NCP EN,λ (k)| (16)
Both L-curve and NCP criteria (λ N CP W N eq.15 and λ N CP EN eq.16) will be investigated in synthetic tests and experimental data in sections 3 and 4.

Tests on Numerical Simulations

In this section, we will test and compare the above strategies on synthetic tracks from convected particles in a Direct Numerical Simulation (DNS) in a triply periodic box. Gaussian noise is added on the particle trajectories so as to have both noisy trajectories and their references to quantify the filter behavior and compare the TrackFit tuning strategies.

For the numerical simulation, we use a well resolved Direct Numerical Simulation of a forced Taylor Green Turbulent flow at R λ = 62. A pseudo-spectral scheme on a 768 3 periodic cubic box of length 2π was used to solve the Navier-Stokes equations. The forcing is introduced through the Taylor-Green modes. Ideal tracers were randomly introduced in the flow and convected along the Eulerian flow field by means of a classical 4 th order Runge-Kutta scheme and linear interpolation. Tracers' positions were saved at a given frequency f DN S such that f DN S τ η ≈ 1.2, τ η being the Kolmogorov time. Unless otherwise specified, we added white noise on the Lagrangian track positions using a Gaussian distribution of µ N oise = 0 mean and σ N oise = 0.1η standard deviation where η is the Kolmogorov lengthscale. The smoothing algorithm was tested on trajectories of length N = 200. Although more than 2.10 6 had been introduced in the DNS, the study in this section was only done on 50 000 trajectories which were found sufficient for the purposes of our statistical analyses.

First, we will describe the computation of the different regularization parameters namely λ c , λ d , λ L , λ N CP W N and λ N CP EN . Then, we will use error statistical analysis, position, velocity and position SNR spectral behavior as well as Lagrangian Structure Functions to evaluate their respective performances.

Criteria computation

Fig. 3a shows the Power Spectral Density (PSD) of the particle X positions from the DNS, the noisy data, as well as the PSD of the added noise. From this spectrum, we can define f SN R=1 , yielding a regularization parameter λ SN R=1 following eq. 10. The f SN R=1 frequency cannot be directly computed from the experimental data since we do not know the real noiseless physical spectrum nor do we know the noise spectrum. Fig. 3b shows the computation of the λ c parameter from the corner detection (corresponding to the red star). The λ d criterion was computed using the yellow segment of the raw spectrum and corresponds to the intersection of both dotted lines. The exploration of different regularization parameters allows for the computation of the L-curve (Fig. 4a) as well as its maximum curvature yielding λ L . The NCP of the residual vector r λ for different λ is shown in Fig. 4b. Since the noise is a white Gaussian noise, we estimated the noise spectrum as a horizontal line whose intensity E N oise was estimated from the flat tail of the spectrum. λ N CP W N is estimated from the comparison between the NCPs and the identity line. The different regularisation parameters as well as the corresponding cutoff frequencies are summed up in table 1.

λ SN R=1 λ c λ d λ L λ N CP W N λ N

Statistical and spectral error analysis

The statistical analysis of the filter normalized errors can be done first by computing the mean and the standard deviation of the errors for different regularization parameters. The filter errors are based on the particle position, velocity and acceleration, respectively normalized by η, η/τ η and η/τ 2 η . To compute the statistics, the first and last 10 points of each trajectory were excluded to avoid border issues. One could argue that with increasing lambdas, ie increasing filter width, 10 points may not be enough to account for border issues. However, tests (not included in this article) showed clearly that the error statistics were but only slightly impacted by this choice at regularization parameters well above the optimum one.

The mean normalized errors based on position, velocity and acceleration were all found to be O(10 -5 ) which indicates that TrackFit does not create any biases since errors are centered on 0. Fig. 5a shows the normalized error standard deviations based on either the particles positions, velocities or accelerations with respect to the regularization parameter. Each curve has a "v" shape curve : for low λ, there is scarcely any smoothing, the error on the position is of the order of the noise standard deviation (0.1η in our case), and it is much higher when considering errors on velocity and acceleration. For high λ, high signal frequencies are far too damped which leads to high errors. For each curve, a minimum can be found. The curve minima for the standard deviations of the errors are all located in the close vicinity of λ SN R=1 . The λ c criterion seems to severely under-smooth the signal and leaves too much noise unfiltered, while in contrast, the λ N CP W N criterion over-smooths the signal. Since λ N CP EN and λ d are the closest to λ SN R=1 , their errors are the lowest. However λ d seems to have higher dispersion of acceleration error than λ N CP EN . The λ L criterion has higher position error dispersion than λ d and λ N CP EN but relatively comparable dispersion error on the velocity and acceleration to λ N CP EN . To better understand the filter behavior and the error dispersion, the Probability Density Functions (PDF) of the acceleration errors was computed in Fig. 5b. An ideal filter would have the narrowest error PDF centered around 0. In our case, one sees that the error PDFs which have a high value in 0 have wide tails, as in contrast, error PDFs with lower tails have a low value at 0. This behavior is directly linked to how TrackFit reconstructs the signal gradients. To be able to reconstruct strong gradients, the regularization parameter need not be too high, thus resulting in low PDF tails, with low errors in areas of strong acceleration. However, this will lead to under-filtering in regions of the low gradients and an increase of low-level acceleration errors which have a considerable statistical weight as PDF of Lagrangian acceleration are centered on 0. This means that a trade-off has to be found between extreme high gradients reconstruction and low gradients area of the particle trajectories. The λ N CP W N criterion is the one which has the highest PDF value at 0. However, this choice corresponds to an over-filtering of the acceleration leading to wider tails of the PDF. Additionally, one sees that the λ c criterion creates a Gaussian-like distribution of the acceleration error. This behavior is due to the unfiltered Gaussian noise which is still left in the signal. The best regularization criterion for the acceleration is the one closest to the λ SN R=1 which is λ N CP EN representing a trade-off between high gradients reconstruction and low gradients smoothing.

A spectral analysis of the trajectories and velocities will enable us to further understand the ability of the filter to recover the turbulent small scales. Fig. 6 shows the PSD of positions and velocities. Since the trajectories are non periodic signals, spectra are periodized using a point symmetry of the signal [START_REF] Foucaut | PIV optimization for the study of turbulent flow using spectral analysis[END_REF]. Fig. 7a shows the impact of the noise on the particle positions. The choice of the regularization parameter does not lead to very different results except around f SN R=1 . The corner detection strategy, as expected, does not filter enough. On the other hand, the λ N CP W N leads to an over-filtering around f SN R=1 .

To distinguish between the other strategies, the regularization choice can be better evaluated from the velocity spectrum. With λ SN R=1 , the smoothed velocity spectrum can recover the correct energy until f ≈ 4 Hz < f SN R=1 where the spectrum falls below the reference one while still being very similar. λ N CP EN and λ SN R=1 lead to similar results when considering the spectrum. It is interesting to realize that the λ d strategy leads to an energy spectrum which is the closest to the DNS reference down to f ≈ 10 Hz. However, we have seen that it does not mean that this strategy minimizes the velocity errors (Fig. 5).

To further our understanding, we need to investigate the error spectrum behavior and find how it compares to the real signal spectrum behavior. The idea is to understand, in the frequency space, how the real physical energy level compares to the error energy levels made by the filtering operation. Similar to the Signal-to-Noise ratio (SNR) terminology used above, we can define a Signal-to-Error ratio (SER) spectrum computed by normalizing the DNS reference spectrum by an error spectrum. For a given frequency, if the SER is above one, it means that error levels are below the real signal. On the contrary, if the SER is below one, it means the error levels of the filter are above the level of the real noise. This will help us quantify, frequency by frequency, how the filter is able to remove noise, and how far in the frequency space, the filter has the ability to denoise effectively and with precision. Furthermore, the error considered can be the error on the particles position, velocity or acceleration. In Fig. 7, we show the filter's position and velocity SER. One sees that the filtering method allows us to access frequencies above the f SN R=1 since SER = 1 is reached for f ≈ 6 Hz for the best regularization strategy both for position and velocity. The best strategy is the one that leads to the largest SER ≥ 1 on a wider range of frequencies. The λ N CP W N deteriorates the SER for the lowest frequencies by over-smoothing the signal. The λ c is an excellent choice for the small frequencies but the SER degrades rapidly due to the under-smoothing even well before the f SN R=1 . The L-curve criterion has a slightly lower SER than λ d and λ N CP EN for frequencies between 2 and 5 Hz. However, its SER is the highest for frequencies larger than 5 Hz. The SERs for λ d and λ N CP EN are the same and above the others. Starting from f SN R=1 , the SER for λ N CP EN is slightly above the SER for λ d and the SER is kept above 1 until f ≈ 6 Hz. Those detailed trends are the same for the SER based on position and velocity. Overall this shows that the λ N CP EN is the regularization strategy that has the best SER for both position and velocity.

Analysis of the trajectories statistical geometric properties

This section is devoted to the assessment of the Lagrangian smoothing method on the trajectories geometry. Following [START_REF] Braun | Geometry of particle paths in turbulent flows[END_REF] we computed the curvature κ and the norm of the torsion θ for each trajectory.

κ(t) = u ∧ a |u 3 | , θ(t) = |u.(a ∧ ȧ)| κ 2 u 6 (17)
where u, a and ȧ are the particles velocity, acceleration and acceleration derivative or jerk respectively and u = ||u||, and ∧ the cross-product operator between two vectors. For the DNS and noisy DNS, velocities and acceleration were computed using a second order centred difference scheme. For all signals, the jerk was computed by a 2 nd order central finite difference scheme from the positions.

The PDFs of curvature and torsion are shown in Fig. 8 where 2.10 5 trajectories were used to allow for a good statistical convergence of the results. The curvature PDF exhibits the same power laws as in the numerical analysis of [START_REF] Braun | Geometry of particle paths in turbulent flows[END_REF], [START_REF] Scagliarini | Geometric properties of particle trajectories in turbulent flows[END_REF]) as well as the experimental results of [START_REF] Xu | Curvature of Lagrangian Trajectories in Turbulence[END_REF]). The PDF of curvature behaves as κ for low curvatures and as κ -5/2 for high curvatures. Those robust power laws were shown to be independent from turbulence and linked to the relatively independent random Gaussian-like nature of velocity variables in turbulence [START_REF] Xu | Curvature of Lagrangian Trajectories in Turbulence[END_REF]). This is visible in Fig. 8a since the raw unfiltered noisy data also exhibit the -1 and -5/2 power laws. However, the position of the PDF peak is strongly linked to the turbulence and Re λ T , the Reynolds number based on the Taylor microscale. [START_REF] Xu | Curvature of Lagrangian Trajectories in Turbulence[END_REF] suggests that the peak location may be due to the intermittency of the acceleration and to the velocity and acceleration correlation in a fully turbulent field. Fig. 8a shows that all five regularization criteria (λ c , λ d , λ L , λ N CP W N , λ N CP EN ) are able to reproduce the PDF of the DNS, thus nicely resolving the curvature of trajectories. Only the PDF of curvature using λ c seems to be shifted to the right towards the PDF of the noisy data.

The computation of the torsion has always challenged experimenter since it includes the derivative of the acceleration, the particle jerk, which is extremely sensitive to experimental noise thus requiring a good time resolution of the measurements to distinguish noise from physical signal. Since the DNS used in those synthetic tests is well resolved, it seemed relevant to compute the torsion that can be extracted from the smoothing algorithm (see Fig. 8b). The absolute torsion PDF shows the same behavior as numerical results based on the trajectories of an ABC flow [START_REF] Braun | Geometry of particle paths in turbulent flows[END_REF]. A flat slope is seen for low torsion values and a -3 slope for higher torsion. Regarding the regularization parameters, all three criteria are capable of capturing this behavior. A closer inspection shows that the λ N CP EN criterion is the one criterion which follows the best the DNS torsion PDF. λ c seems to create higher torsion values than the DNS resulting in a slight shift of the PDF towards higher torsion values similarly to the noisy data PDF. The other parameters (λ N CP W N , λ L and λ d ) lead to similar PDF. This is a promising result in that TrackFit has the ability to recover the particle jerk. Ironically enough, one of its operating principles is to regularize the particles jerk which can lead to an almost null jerk if the regularization parameter is set too high. This preliminary work needs further investigations to quantify the quality of the jerk measurement but this lies beyond the scope of our article.

Lagrangian Velocity Structure function analysis

The aim of this section is to understand the impact of the trajectory smoothing algorithm and its parameterization on the Lagrangian velocity structure functions which are high order moments of velocity increments.

Lagrangian velocity structure functions of order p are defined as :

S p (τ ) = |e. (v(t + τ ) -v(t)) | p (18)
The . here refers to ensemble averaging and averaging over all possible directions of the axes, i.e. the random unit vector e [START_REF] Berg | Experimental investigation of Lagrangian structure functions in turbulence[END_REF]. They are of particular interest when trying to understand the statistical behavior of small scales in turbulence. Moreover, their scaling laws have been hotly debated because the Kolmogorov theory of 1941 (A.N. [START_REF] Kolmogorov | The Local Structure of Turbulence in Incompressible Viscous Fluid for Very Large Reynolds[END_REF] predicted in the inertial range a S p (τ ) ∼ ( τ ) p/2 scaling, with the dissipation rate, thus predicting statistical similarity of velocity increments in the inertial range for homogeneous isotropic turbulence. However, several studies have revealed deviations from the statistical self-similarity [START_REF] Mordant | Measurement of Lagrangian Velocity in Fully Developed Turbulence[END_REF]. This phenomenon is associated with Lagrangian intermittency which is found to be larger than in the Eulerian framework [START_REF] Xu | High Order Lagrangian Velocity Statistics in Turbulence[END_REF]).

We will focus in this section on the small scale behavior of the Lagrangian structure functions since the low-pass behavior of TrackFit should not impact the larger inertial scales of the DNS. In Fig. 9, we show the 2 nd (left top) and 6 th (left bottom) order structure functions as well as their respective local scaling exponents. The aim is to understand the impact of the regularization strategy on the structure functions. At small scales, the noise on the velocity measurement creates a well-known bias in the Lagrangian statistics [START_REF] Berg | Experimental investigation of Lagrangian structure functions in turbulence[END_REF]).

The impact of smoothing is thus mostly visible at small scale (τ ≤ 2τ η ). Furthermore, the smoothest signal has the lowest structure function values at small scales. For small scales, τ < τ η , all regularization strategies recover the expected S p (τ ) ∼ τ p scaling stemming from regularity. For the second order structure function, there are very little differences between the regularization strategies and all seem to recover well the DNS reference. The only slight difference appears in the S 2 scaling exponent (Fig. 9b). The choice of λ c leads to a lower scaling exponent, which comes from the fact that a signal smoothed with λ c is more noisy than its DNS reference. At small scales, the corresponding structure function is thus higher than the DNS. On the contrary, the λ N CP W N leads to a higher exponent than for the DNS. The λ N CP EN seems to be the one which best recovers the scaling exponent of the DNS. However, differences are more pronounced for the 6 th order structure functions as the impact of the smoothing increases with the order of the structure functions. What is surprising is the fact that S 6 for the noisy data recovers the DNS reference very well up to τ ∼ 0.4τ η .

Similarly, the λ c criterion is the one which is the closest to the DNS reference. This behavior is counter-intuitive since we saw that the choice of λ c is not able to filter the noise as well as the other criteria. We believe that this behavior is due to the fact that the higher the order is, the more sensitive the structure functions are to high values of velocity increments which may be smoothed out in TrackFit at small scales but are present in the noisy data.

This means that the seemingly good behavior of λ c is in fact a by-product of the noise which is still in the signal. The trend might be good but is due to the measurement noise. This shows how difficult it is to use such statistical criteria as performance criteria since seemingly good trends actually can come mainly from a lack of measurement precision.

Signal-to-Noise ratio analysis

Up to now, only one Signal-to-Noise ratio was considered in our numerical tests (σ N oise = 0.1η). In a real experiment, the noise level is often fixed in pixel size. Indeed, once the optical setup is set (calibration, magnification, pixel size, particle image size, particle density, defocusing effects, laser power, temporal and spatial illumination for instance), the level of noise on the particle trajectory will mostly be fixed compared to the pixel size. However, the values of position error relatively to turbulent scales increase with the Reynolds number. The higher the Reynolds number is, the smaller η is compared to the pixel size resulting in a higher SNR relatively to η.

To understand how the TrackFit filter behaves with higher SNRs, we computed mean normalized accelerations errors as a function of λ for different standard deviation noise (σ N oise ) ranging from 0.01η to 30η in Fig. 10. Each grey curve corresponding to the error for a given SNR has a minimum. The best strategy is the one which is able to capture each curve minimum. In yellow are the errors when using λ SN R=1 . One sees that they correspond to the error curves minima, except for very small SNRs where they are in the close vicinity of the minima. For the sake of simplicity, we only show here the acceleration errors for two regularization strategy which were shown to give good results in the previous tests : λ d and λ N CP EN . Fig. 10 shows that the value of the regularization parameter for the λ d strategy is always lower than the optimum and the gap between λ d and the minimum tends to grow as the SNR decreases. Though λ N CP EN is slightly less than the optimum one for the first two low SNRs, it is clearly the most capable of finding the optimum regularization parameter for a wide range of SNR, including very high SNRs.

Overall, the λ N CP EN strategy is the best strategy to minimize mean acceleration errors, and also position and velocity errors since Fig. 5 showed that the optima of the three quantities were found for the same regularization parameter.

Experimental results

Experimental set-up

Time resolved measurement of 4D-PTV "Shake-The-Box" [START_REF] Schanz | Shake-The-Box: Lagrangian particle tracking at high particle image densities[END_REF] were performed on a turbulent Von Kármán flow in water [START_REF] Ostovan | 4D Particle Tracking Velocimetry measurements in a Von Karman turbulence experiment[END_REF]. Four high speed cameras are imaging a 45 × 40 × 6 mm 3 volume lighted by a 30 mJ/pulse high speed laser Nd-YLG laser. The measurement volume is located at the center of a cylindrical water tank of radius R = 0.1 m and height H = 0.47 m. It is worth mentioning that all four cameras are in the same Y = 0 plane, the Z axis being the reconstructed dimension. The Y axis corresponds to the Von Kármán rotation axis. Energy is being injected in the fluid by two counter-rotating impellers located respectively on the top and at the bottom of the water tank. For the present analysis, the frequency is set at F = 0.1 Hz, resulting in a Reynolds number based on the tank radius Re = 2πF R 2 /ν ≈ 6300. The Kolmogorov length (η = 0.3 mm) and time scale (τ η = 92 ms) were computed from an estimate of the global mean dissipation rate in the water tank from torque measurements. The measurement frequency f m = 200 Hz was set to resolve the Kolmogorov time scale. The Taylor length scale is λ T ≈ 5.3 mm which leads to a Reynolds number based on the Taylor lengthscale Re λ T ≈ 80. For a thorough description of the experimental setup, a complete description of the experimental set-up is provided in [START_REF] Debue | Experimental approach to the problem of the Navier-Stokes singularities[END_REF].

Glass hollowsphere particles of 10 µm averaged diameter from DANTEC were inserted in the water tank resulting in a Kolmogorov based Stokes number of St τη = 6.57 10 -5 . The particle image density was estimated at 0.047 particles per pixels resulting in an average of 40 000 particles tracked per time step in the measurement volume (Fig. 11). Particle reconstruction and tracking was obtained using the STB algorithm of Davis10 software. The parameters used for the analysis are detailed in [START_REF] Ostovan | 4D Particle Tracking Velocimetry measurements in a Von Karman turbulence experiment[END_REF].

Regularization criteria : computation and remarks

In this section, we discuss the results of the aformentioned regularization criteria on the experimental data. Their values are summed up in table 2. Fig. 12 shows the spectrum of the trajectories raw positions. It was computed using N tracks = 59 222 tracks of N length = 100 time steps long. Looking at the Lagrangian spectrum gives us a lot of information concerning the temporal resolution of our measurements and the level of error in the particle position. The level of the spatial noise is seen by looking at the level of the noise section of the spectrum, at high frequencies. The temporal resolution is given by the frequency at which the noise appears in the spectrum. The position error can be estimated in between 0.01 to 0.02η. The temporal resolution of our measurement is about τ η .

Several additional remarks can be made regarding our results. First of all, one observes a clear difference between the X, Y dimensions and the third dimension Z. The reconstructed Z dimension seems to be more noisy by almost one order of magnitude than the X and Y dimensions (for f ≥ 10 Hz). We believe that this is due to the linear imaging configuration of our camera setup. The Z axis is the axis which is the most co-linear to all the cameras optical axis. This dimension is thus the most difficult to reconstruct during the 3D particle reconstruction process and therefore the dimension on which the position errors will be the highest. This is why cross-like configurations are preferred whenever possible [START_REF] Scarano | Tomographic PIV : principles and practice[END_REF]. The Y component is also less noisy than the X one since the X spectrum is slightly above the Y spectrum in the frequency range affected by the measurement noise. Furthermore, the frequencies at which the noise appears in the signal is not the same for the two largest dimensions (X, Y ) with respect to the Z dimension. Thus, the smoothing cutoff frequencies have to be separately chosen for each direction based on each direction noise spectral behavior. Finally, a striking feature of the spectra is the fact that the high frequency noise is not a pure white noise. The reason for this is still unknown but is certainly to be found in a combination of the experimental setup, the particle reconstruction and the tracking algorithm since the Z spectrum is even less flat than the X and Y spectra. We believe that this feature could be a footprint of the Wiener filter used for predicting the position of a tracked particle in the next time-step [START_REF] Schanz | Shake-The-Box: Lagrangian particle tracking at high particle image densities[END_REF]. The slightly correlated nature of this noise could also be due to the illumination inhomogeneity and two head-pulsed laser misalignement causing the particle position error to be correlated in space and time.

The f cross frequency, defined in 2.4, is therefore difficult to find especially for the computation of λ d , which requires the experimenter to approximate the noise. Fig. 12 shows that neither the signal part nor the noise part of the spectrum can easily be approximated by a straight line. We therefore only computed the λ c criterion based on the curvature of the spectrum for each direction. We also computed the L-curve and extracted their related curvature maxima (see Fig. 13a) yielding for each direction λ L,x , λ L,y and λ L,z . We recover the fact that the Y direction is the less affected by the noise with λ L,y ≤ λ L,x λ L,z .

We also computed the two NCP criteria : N CP W N , which is based on the assumption of a pure white noise and N CP EN which is based on a comparison of a current residual spectrum and a model of the best residual spectrum. The model of the best residual spectrum was chosen to be made of two distinct parts :

E λ residu (f ) = E signal (f ) if f ≥ f λ E signal (f λ ) (1 -H λ (f )) otherwise (19)
The λ N CP EN is the λ for which the NCP residue is closer to the NCP of its associated model residue E λ residu (f ). Fig. 13b shows for the Z dimension the tested NCPs of the residual vectors of TrackFit as well as N CP W N , N CP EN and their associated models. Table 2 summarizes the results for all the different regularization parameters. On the Z dimension, one sees that λ N CP EN λ N CP W N . Choosing to estimate the noise as white noise in the case of strongly non-uniform noise may lead to a wrong and under-estimated smoothing parameter value. This will be further examined in the following section.

X Y Z λ f c [Hz] λ f c [Hz] λ f c [

Results

Experimental results are brought together in this section. We will first look at positions and velocity spectra of the particles tracks. We focused on the Z direction, since we showed that this dimension is the more impacted by the measurement noise.

The spectra of Z positions are shown in Fig. 14a for the different regularization parameters. An overall -6 slope is observed for all spectra, as previously found on numerical data in section 2.3. There is almost no quantitative difference between the regularization parameters, except for λ N CP W N which as expected, has a higher cutoff frequency and does not perfectly filter the noise as visible in the small bump around f = 10 Hz. The Z velocity spectrum (Fig. 14b) helps to discriminate between the different parameters. The noise amplification on derivatives of a noisy signal is a well known issue, and depends on the derivative scheme used. It is clearly visible in the spectrum of raw unfiltered data where a 2nd order centered difference scheme was used to compute the derivatives from the positions. The velocities of TrackFit were computed from the derivation of equation 3 resulting in a linear combination of polynomial functions of order n -1. A noise-induced bump is visible in the λ N CP W N velocity spectrum as well as for λ c around f ≈ 10 Hz. This was expected for λ N CP W N which totally overestimated the cutoff frequency. The spectrum behavior of the noise clearly differs from a pure white noise, thus rendering the NCP estimate and comparison incorrect. For λ c , the bump around f ≈ 10 Hz shows that the geometrical criteria selecting the cutoff frequency was not conservative enough regarding the noise estimation since the signal remains noisy. On the contrary, the λ L and λ N CP EN do not exhibit such a bump and their respective velocity spectra do not seem to be affected by the measurement noise. Since we do not have a reference signal or another independent measurement for comparison such as in [START_REF] Lawson | Bias in particle tracking acceleration measurement[END_REF], it is difficult to state that one is better than the other. However, one could argue that λ L may smooth the signal more than λ N CP EN . Though, this trend is found both in the numerical tests in 3 as well as in this experimental dataset, we do not expect this to be always the case for any given experimental dataset. Therefore, we strongly advise experimenters to apply such treatments to their data set in order to best select the optimal filter frequency cut-off. In any case, the optimum is definitely around this value. Further investigation could also be made on the impact of this unfiltered noise on the interpolation of the Lagrangian velocity field on an Eulerian mesh, but this lies well beyond the scope of our present research.

The particle Lagrangian acceleration statistics are crucial quantities in Lagrangian turbulence and often used as a criteria to select the filtering width. We computed the particle acceleration PDF (Fig. 15.) for the three space coordinates of the acceleration and for the 4 regularization strategies considered in section 4. To avoid any biases due to the possibility (which is now known in the 4D-PTV community for the Shake-The-Box algorithm) of ghost tracks, we only considered relatively long tracks of more than 50 time steps. We also removed the first and last 10 points of the tracks to avoid any border issues. For a good statistical convergence, we computed the statistics over 40 runs of 3226 time instants each. This resulted in more than 7.10 9 data points. The PDF is non gaussian, as expected in Lagrangian turbulence [START_REF] Mordant | Experimental lagrangian acceleration probability density function measurement[END_REF]), [START_REF] Voth | Measurement of particle accelerations in fully developed turbulence[END_REF], [START_REF] Berg | Experimental investigation of Lagrangian structure functions in turbulence[END_REF]. Fig. 15 shows the PDF is well converged up to 10 -6 . The PDF of the Z acceleration is also quite symmetric for a range of normalized acceleration between -18 and +18. The tails of the PDF on this component are not symmetric, and in this range, the choice of regularization parameter can be seen. Since our measurement volume is at the center of the Von Kármán flow, there should not be any preferential flow direction.

The asymmetry of the PDF in its tails is certainly due to the measurement noise, which is an order of magnitude higher than in the X and Y direction and the PDF done on X and Y show symmetric tails (not show here).

As excepted from the previous investigations, the λ c and λ N CP W N are the ones which allow for higher accelerations with wider PDF tails. λ N CP EN and λ L regularization parameters produce smoother results. If one is interested in locating extreme events in the turbulence, less regularization will definitely allow to capture higher accelerations. However,it will also allow for smaller frequencies to have higher energy values, as seen in the velocity spectra. This will result in higher averaged errors as seen in Fig. 5 Geometrical properties of the Lagrangian trajectories such as PDF of curvature and torsion are shown in Fig. 16). All the behaviors observed in the numerical tests are recovered in this experimental application. The robust +1 and -5/2 slopes are recovered in the PDF of curvatures. The shift of the PDF towards lower curvature due to the smoothing procedure is also observed. All the different regularization parameters exhibit almost the same PDF of curvature. The PDF of torsion (Fig. 16b) is flat for low torsion values and exhibits a -3 slope for high torsion values. The transition between these two states is also affected by the smoothing procedure. One finds the same trend as observed in section 3.3 : the signals which are the more regularized (λ N CP EN and λ L ) have slightly lower torsion values. Overall, such a statistical criterion shows very little difference between the regularization parameters.

Velocity structure functions were also computed and we show (Fig. 17a &17c) the 2nd and 6 th order velocity structure functions as well as their scaling exponents (Fig. 17b &17d). The scaling exponents of the DNS are shown for comparison. Since we are only looking at relatively short Lagrangian tracks of 100 time steps, we do not observe any plateau as expected for S 2 and S 6 indicating that the inertial range was not covered in this Lagrangian exploration of our turbulent Von Kármán flow. All four regularization parameters are able to recover and regularize the small time scale velocities since all the curves superimpose in Fig. 17a and17c. Furthermore, the scaling exponents recover the expected trends of a regular differentiable function on which a Taylor expansion can be done. For the scaling exponents, differences arise for the 6 th order : λ L and λ N CP EN have slightly higher scaling exponents than λ c and λ N CP W N for τ ≤ τ K . This is coherent with the trend observed on the DNS where higher values of smoothing power lead to higher scaling exponents. The comparison with the DNS shows that the experimental structure functions have higher scaling exponents for -0.75 ≤ log 10 (τ /τ K ) ≤ 0.4. One can think of numerous reasons that might explain this difference. First, the flows are rigorously not the same since the flow geometry and the mechanical properties like Re λ T differ from one another. Furthermore, it is well known that PTV measurements can suffer from multiple biases like inhomogeneity of seeding, errors in the case of high seeding and high gradients or anisotropy of the measurement volume to name a few. In our case, the measurement laser volume has a smaller Z expansion than in the X and Y dimension. Furthermore, to avoid border issues, the structure functions were computed on trajectories larger than 50 time steps. The particles that have high Z velocities tend to come out of the volume and belong to short trajectories and thus are scarcely represented in the statistics.

Conclusion

We compared three types of optimal cutoff frequency finding strategies for the smoothing of noisy Lagrangian trajectories. Those strategies are based on the shape of the spectrum (corner detection), the shape of the regularization L curve problem (corner detection) and the comparison between the spectral behaviors of the residual vector and the estimated noise. Comparisons were done using DNS assuming a Gaussian white noise on the trajectories and also on experimental data of a Von Kármán flow at Re ≈ 6300.

Numerical simulations, assuming a Gaussian white noise, gave us an in-depth understanding of the types of errors that occur with this smoothing technique by computing mean absolute error and standard deviation error on positions, velocities and accelerations. We can conclude from this research that an optimal regularization parameter does exist for each physical quantity while the dispersion of the error is minimum for positions, velocities and accelerations for λ close to f SN R=1 . We show that this behavior is linked to the way low or high gradients of the signal are reconstructed by looking at PDFs of accelerations and acceleration errors. Low regularization parameters will faithfully reconstruct high gradients while adding non-physical fluctuations in low-gradients areas. On the contrary, higher regularization parameters will smooth out high gradients and better predict low gradients. We showed that regularization parameters based on the shape of the position spectrum (especially λ c ) tend to under-smooth the signal, thus recovering high accelerations but adding artificial fluctuations to low gradients. In comparison, we show that λ N CP W N over-smooths the signal. λ d , λ L and λ N CP EN have quite similar performances and a combined analysis of error criteria and SER spectrum is needed. λ d is able to recover the acceleration variance quite accurately but it fails to reconstruct low acceleration. It results in a lower SER for frequencies higher than f SN R=1 . λ L is less accurate for high accelerations but reconstructs low intensity gradients well thus resulting in the highest SER at high frequencies (f ≥ f SN R=1 ). λ N CP EN seems to be a good compromise between those two tendencies. Those results are confirmed by the analysis of the statistical geometrical properties of the Lagrangian trajectories where λ N CP EN is the parameter which recovers the best behavior of the torsion PDF behavior. The investigation of the Lagrangian velocity structure functions showed that λ N CP EN recovers the closest the scaling exponent of the DNS for the second order structure function. SNR analysis showed that for a wide range of SNR, λ N CP EN is able to find the regularization parameter leading to the least mean acceleration errors.

When confronted to real experimental data, the first finding is that the experimental noise can be quite different from a white Gaussian noise, especially for the smallest reconstructed dimension. At this stage, this is an ongoing research, but it could be extremely stimulating to understand the origin of this phenomenon in order to improve the measurement accuracy. This does undoubtedly impact the possible strategies to distinguish the real signal from the noise. Indeed, we show that the non-Gaussian noise renders the λ W N strategy ineffective. Furthermore, the geometrical properties of the position spectrum change, making λ d very hard to set objectively. We show that, both in the case of Gaussian noise and experimental noise, λ c under-smooths the signal. It does furthermore impact the statistical Lagrangian geometrical properties and this unfiltered noise could be detrimental to any interpolation into the Eulerian flow field. λ L and λ N CP EN have very similar performances and behavior for velocity spectrum, acceleration PDF, PDF of curvature and torsion and also velocity structure functions. This robustness combined with the performances obtained in the numerical simulation strengthens the idea that those parameters constitute the best possible tuning of the TrackFit method. The L-curve criterion, though purely based on signal processing analysis, is quite robust and offers the experimenter a good systematic cutoff frequency finding strategy but leads to smoother results than with λ N CP EN .

Overall, this study showed that TrackFit is a very powerful tool to recover denoised Lagrangian trajectories. Furthermore, the smoothing tuning strategies can be based on a desired accuracy and on the physical behavior of the measured quantities the experimenter wishes to focus on. For instance, extreme events with rapid change of direction can more readily be recovered with lower constraints on the regularization. If the position spectrum slope has a clear negative slope, then the λ d strategy can be applied and is still close to the optimum regularization parameter for mean errors. As shown in acceleration PDFs, higher acceleration values can be recovered. However, this is done at the detriment of lower frequency accuracy and may lead to more noise in low acceleration areas when interpolating the field on an Eulerian mesh. If the experimenter wishes to use TrackFit as a first step before interpolation to obtain the Eulerian flow field, or to reconstruct the pressure field using the Lagrangian acceleration, the idea is to choose a regularization parameter which has the best performances overall on all possible scales, thus λ L and λ N CP EN are the best candidates.

Ultimately, this smoothing technique requires to find a trade-off between low and high acceleration reconstruction. This can be detrimental to the measurement of intermittency which is a crucial property of turbulence. This choice is required because the same smoothing level is applied to all tracks and to the whole track length. To go further, an idea could be to tune the regularization parameter according to an estimate of the intensity of local velocity gradients. This would allow a stronger regularization parameter in low acceleration areas, and a lower regularization when high frequencies with high amplitudes are detected. This promising idea may lead to significant gains for the measurement of turbulence.

Figure 1 :

 1 Figure 1: a : Filter transfer function for different regularization parameters λ. The analytical model from eq. 10 is shown here for λ = 100 b : Link between the analytical f c from eq. 11 and the effective -3 dB frequency cutoff measured from the numerical transfer function.

  Figure 2: a : Representation of an L-curve b : NCP for a white noise signal as well as its low-pass and high-pass filtered version. The dotted line represents the x = y identity line.

  Figure 3: a : Power Spectral Density of particle positions for the DNS, noisy data and the added Gaussian noise. b : PSD of particle x-positions showing the spectrum corner detection from the maximum spectrum curvature (the red star) corresponding to the λ c parameter. The λ d parameter is computed from the intersection of the dash lines. The yellow segment was chosen for the line fitting resulting in the dash line.

  Figure 4: a : L-curve of the Lagrangian smoothing algorithm. The black star represents λ L computed from (13). b : NCP of the residual vector for different regularization parameters. In full red, the NCP closest to the x = y line corresponds to λ N CP W N . In full black, the residual NCP is associated to the λ N CP EN criterion with the corresponding NCP EN,λ in dotted black line. The black arrows represents the direction of increasing λ.

  Figure 5: a : Normalized standard deviation errors based on position, velocity and acceleration errors. Position errors were normalized by η, velocity errors by η/τ η and acceleration errors by η/τ 2 η . λ c and λ d are respectively represented by the and star marker, λ L is indicated by the marker, λ N CP W N and λ N CP EN are represented respectively by the and square marker. The vertical dotted line indicates λ SN R=1 . The color code associated with each regularization strategy is respected throughout the article. b : Acceleration error Probability Density Functions for different regularization parameters.

Figure 6 :

 6 Figure 6: PSD of X positions (a) and X velocities (b) for the reference DNS, the noisy data and the results of TrackFit using the different strategies.

Figure 7 :

 7 Figure 7: SER spectra based on X positions (a) and X velocities (b) for TrackFit results using the different tuning strategies.

  Figure 8: a : PDF of the trajectories curvature b : PDF of the norm of the trajectories torsion.

  Figure 9: a : Lagrangian velocity structure function of order 2, S 2 , as a function of nondimensional time increment b : Local scaling exponent d[log(S 2 )] d[log(τ /τη)] as a function of nondimensional time increment c : Lagrangian velocity structure function of order 6, S 6 , as a function of non-dimensional time increment d : Local scaling exponent d[log(S 6 )] d[log(τ /τη)] as a function of non-dimensional time increment.

Figure 10 :

 10 Figure 10: In grey, mean normalized acceleration errors as a function of λ for different SNR, ie σ N oise (0.01η, 0.05η, 0.1η, 0.5η, 1η, 5η, 10η, 30η). The yellow stars represent λ SN R=1 . The purple and red lines respectively correspond to λ d andλ N CP EN .

Figure 11 :

 11 Figure 11: Visualisation of Lagrangian trajectories color coded by the velocity norm.

Figure 12 :

 12 Figure 12: Raw particle positions spectrum for X, Y and Z.

  Figure 13: a : L-curve curvature for X, Y and Z dimensions b : Normalized cumulative periodogram of residual vectors for the Z coordinates.

  Figure 14: a : Power density Spectrum of Particle Z positions b : Power density Spectrum of Particle Z Velocities.

Figure 15 :

 15 Figure 15: Particles' Z acceleration PDF normalized by its standard deviation for different regularization parameters.

  Figure 16: a : PDF of the trajectories curvature b : PDF of the norm of the trajectories torsion.

  Figure 17: a : Lagrangian velocity structure function of order 2, S 2 , as a function of nondimensional time increment b : Local scaling exponent d[log(S 2 )] d[log(τ /τη)] as a function of nondimensional time increment c : Lagrangian velocity structure function of order 6, S 6 , as a function of non-dimensional time increment d : Local scaling exponent d[log(S 6 )] d[log(τ /τη)] as a function of non-dimensional time increment.

Table 1 :

 1

	CP EN

Summary of the different regularization parameters, their dimensionless values and associated cutoff frequencies.

Table 2 :

 2 Summary of the different regularization parameters, their values λ and associated cutoff frequencies f c .

	Hz]

Acknowledgments

This work was carried out within the framework of ELSAT2020 project supported by the European Community, the French Ministry for Higher Education and Research, and the Hauts de France Regional Council in connexion with CNRS Research Foundation on Ground Transport and Mobility. The experimental data were collected through the ANR EXPLOIT grant agreement no. ANR-16-CE06-0006-01 which is a collaboration between CEA/SPEC and LMFL. This research was granted access to the HPC resources of IDRIS under the allocation 021741 made by GENCI (Grand Equipement National de Calcul Intensif).