Quantitative Analysis of LISA Pathfinder Test Mass Noise

Quantitative Analysis of LISA Pathfinder Test Mass Noise

Luigi Ferraioli luigi@science.unitn.it University of Trento and INFN, via Sommarive 14, 38123 Povo (Trento), Italy    Martin Hewitson Albert-Einstein-Institut, Max-Planck-Institut fuer Gravitationsphysik und Universität Hannover, Callinstr. 38, 30167 Hannover, Germany    Giuseppe Congedo University of Trento and INFN, via Sommarive 14, 38123 Povo (Trento), Italy    Miquel Nofrarias Institut de Ciències de l’Espai, (CSIC-IEEC), Facultat de Ciències, Campus UAB, Torre C-5, 08193 Bellaterra, Spain    Mauro Hueller University of Trento and INFN, via Sommarive 14, 38123 Povo (Trento), Italy    Michele Armano SRE-OD ESAC, European Space Agency, Camino bajo del Castillo s/n, Urbanización Villafranca del Castillo, Villanueva de la Cañada, 28692 Madrid, Spain    Stefano Vitale University of Trento and INFN, via Sommarive 14, 38123 Povo (Trento), Italy
Abstract

LISA Pathfinder (LPF) is a mission aiming to test the critical technology for the forthcoming space-based gravitational wave detectors. The main scientific objective of the LPF mission is to demonstrate test-masses free-falling with residual accelerations below at mHz. Reaching such an ambitious target will require a significant amount of system optimisation and characterisation, which will in turn require accurate and quantitative noise analysis procedures. In this paper we discuss two main problems associated with the analysis of the data from LPF: i) Excess noise detection and ii) Noise parameter identification. The mission is focused on the low frequency region ( mHz) of the available signal spectrum. In such a region the signal is dominated by the force noise acting on test masses. At the same time, the mission duration is limited to days and typical data segments will be hours in length. Considering those constraints, noise analysis is expected to deal with a limited amount of non-Gaussian data, since the spectrum statistics will be far from Gaussian and the lowest available frequency is limited by the data length. In this paper we analyze the details of the expected statistics for spectral data and develop two suitable excess noise estimators. One is based on the statistical properties of the integrated spectrum, the other is based on Kolmogorov-Smirnov test. The sensitivity of the estimators is discussed theoretically for independent data, then the algorithms are tested on LPF synthetic data. The test on realistic LPF data allows the effect of spectral data correlations on the efficiency of the different noise excess estimators to be highlighted. It also reveals the versatility of the Kolmogorov-Smirnov approach, which can be adapted to provide reasonable results on correlated data from a modified version of the standard equations for the inversion of the test statistic. Closely related to excess noise detection, the problem of noise parameter identification in non-Gaussian data is approached in two ways. One procedure is based on maximum likelihood estimator and another is based on the Kolmogorov-Smirnov goodness of fit estimator. Both approaches provide unbiased and accurate results for noise parameter estimation and demonstrate superior performance with respect to standard weighted least-squares and Huber’s norm. We also discuss the advantages of the Kolmogorov-Smirnov formalism for the estimation of confidence intervals of parameter values in correlated data.

LISA Pathfinder; LISA Technology Package; LTP; Noise Analysis
pacs:
04.80.Nn, 95.75.-z, 07.05.Kf

I Introduction

LISA Pathfinder (LPF), a European Space Agency mission, will be used to characterize and analyze all possible sources of disturbance which perturb free-falling test masses from their geodesic motion. The system is composed of a single spacecraft (SC) enclosing a scientific payload, the LISA Technology Package (LTP), which is composed of two test masses (TMs) whose position is sensed by an interferometer. The spacecraft cannot simultaneously follow both masses, therefore the trajectory of only one test mass is used as a drag-free reference along the measurement axis. In order to prevent the trajectories of the test-masses from diverging, the second test mass is capacitively actuated to follow the first (free-falling) TM. In the main science operating mode, the first interferometer channel measures the displacement of the SC relative to the free-falling TM. The second interferometer channel (the differential channel) measures the relative displacement between the two TMs.

LPF is a controlled system which can only be fully assessed during flight operation, therefore a considerable number of experiments will be devoted to the identification of the details of the dynamics of the system. A dynamical model of LPF is built in advance on the basis of physical considerations and from the results of test campaigns. The dynamical model is parametric so that it can be updated on the basis of the experiments that will be conducted during mission operations. The overall aim of the process is to reach the best free-fall quality (below at frequencies around mHz) in a step-by-step procedure in which the result of the previous experiment is used to adjust the detailed configuration of the following experiments LTPVitale2011 (); LTPMartin2011 (); LTPPaul2011 (); ArmanoCQG2009 ().

Such a demanding program requires daily analysis of the instrument signals constrained by two major factors: i) the amount of available data is tightly constrained by LTP mission duration ( days), the telemetry bandwidth, and the length of each data segment (typically h); ii) the scientific interest is mainly focused on the analysis of noise sources which act directly on the TMs since that should provide a baseline reference for the forthcoming space-based gravitational waves observatories LISAYellowBook (); Cornish2005 (); DECIGO (); ASTROD (). The direct forces on the TMs are expected to dominate the instrument output in the frequency range mHz. Sample power spectra are typically calculated with Welch’s averaging periodogram method (WOSA) Percival (). In order to keep enough frequency resolution at low frequencies, the sample power spectra can not be averaged more than few times (we average times in the present paper), this results in highly non-Gaussian data for which we are developing dedicated techniques. In particular the paper aims to propose a solution for two major data analysis challenges encountered in LPF: i) Different measurements of the same physical quantity can exhibit different noise content if they are performed under slightly different environmental conditions. The objective of LISA Pathfinder data analysis during operations will be to discover such differences, understand their origin and adjust spacecraft physical parameters accordingly. Such a problem requires reliable excess noise detection procedures which have to be based on solid statistical considerations; ii) Along with the demonstration of unprecedented test-mass free-fall, LPF will provide a model for the expected test mass force noise for future space-based gravitational wave detectors. In order to do this, we need to be able to match an analytical model to a noisy power spectral density measurement. The quality of the match must be statistically quantified. Both data analysis problems deal with sample spectra and the corresponding statistical properties. Section II reports on the properties of different experimental procedures for the detection of noise excess. In particular we considered two cases where the noise excess is evaluated with respect to reference data or a reference model. The accuracy of the methods is theoretically analyzed for the case of broadband and band-limited excess noise. In Section III the problem of noise parameter estimation for non-Gaussian data is explored and an algorithm based on maximum likelihood is derived. In parallel, the application of the Kolmogorov-Smirnov formalism to the construction of a goodness of fit estimator is discussed. Section IV reports briefly on the extension of the analysis procedures to the case of non-stationary noise and time-frequency investigations. In Section V we provide an application of the developed algorithms to synthetic data with LPF-like qualities. This allows us to shed light on the effects of data correlations on the accuracy of the developed excess noise estimators. The analysis procedures are available as MATLAB tools in the framework of the LTPDA Toolbox HewitsonCQG2009 (); matlab (); toolbox (); an object oriented MATLAB Toolbox for advanced data analysis.

Ii Quantitative detection of noise power variations

The problem of the detection of noise power variation in consecutive measurements can be formulated in two different ways. i) Two different measurements of the power are compared; ii) The different measurements of the power are compared with a reference model. The problem in the first case is of general character and can be applied to a wide range of experiments, the second case, on the other hand, assumes that a reference model for the noise power is known and data must be compared against the given model in order to understand if the system is performing under known conditions. The latter case is likely to be the scenario for LPF operations. The quantity typically used for the detection of noise power variations is the total energy content of the data series. It is defined as . Unfortunately, , provides a poor estimator for two reasons: i) As soon as the data series departs from zero mean Gaussian white noise, the statistic of becomes ill-defined and the definition of a confidence interval becomes cumbersome; ii) provides global information as it is not sensitive to noise changes in a given frequency band. While the first problem could be overcome (without little difficulty) by a numerical identification of the expected statistic, the second problem suggests that a spectral based estimator would provide supplementary information which could be fundamental in discriminating different noise sources.

ii.1 Detection of noise variations with a model

In the case that a reference model is available, the detection of excess noise in the spectral domain can be effectively implemented with a test on the normalized Welch’s overlapped segment averaging (WOSA) spectrum defined in equation (24). In particular, the test can pursue two different philosophies of which one aims to test a global scalar indicator of the properties of the data and another aims to test the details of the statistical distribution of the data.

A sensible estimator for the first approach is provided by the integral of the normalized spectrum which, in the discrete case, can be written:

(1)

In the simplifying assumption of independent spectral data, the statistic of each element of is described by a gamma distribution as defined by equation (23) with and . The sum over the different values at the frequencies is still a gamma distribution with and . The expectation value for is easily obtained as .

The natural estimator for the second approach is provided by the empirical cumulative distribution function (ECDF). The ECDF for a data series can be defined as , where is the number of observations reporting a value less than or equal to . denotes the values taken by the data, in our case with . The ECDF can be tested against the theoretical expectation provided by equation (23). If the model well represents the given sample spectrum then is distributed according to the expected gamma distribution. Alternatively, if the sample spectrum contains a excess noise with respect to the reference model, then the distribution of the normalized WOSA spectrum will be different from the one reported in equation (23). The hypothesis that the two distributions are equal can be tested if a ’distance’ between the ECDF and the theoretical reference, , is defined as:

(2)

with assuming values on the interval and . Kolmogorov found that the statistical properties of

(3)

are independent from the specific distributions under test Kolmogorov1933 (); Feller1948 (). This property qualifies as an excellent candidate for the construction of a general test for cumulative distribution functions, the limiting statistic for was identified by Kolmogorov himself and then inverted by Smirnov Smirnov1939 (); Miller1956 () who provided an analytical expression for the calculation of as a function of the significance level. General details about the Kolmogorov-Smirnov test (KS-test) are reported in appendix B.

It is interesting to calculate the expected sensitivity for the two estimators. is expected to be distributed as a gamma distribution (at least when the model and the data are in agreement), the corresponding cumulative distribution function (CDF) is provided by the incomplete gamma function AbramoStegun (). The CDF assumes values in the interval as it defines the probability associated with the observation . The inverse of the CDF provides the critical values, , associated with a given probability, . The confidence range for a probability, , for a gamma distributed variable can then be defined by the boundaries and , with . We assert that the measured sample spectrum is compatible with the reference model if for the given probability, , or significance level, .

If the noise excess is provided by a scale factor, , which affects the noise on the complete band of frequencies, the expected value for changes to . Therefore the detection threshold for is fixed by the interval . In other words, the can detect a noise difference with respect to the reference model only if or if . If is non-zero only in a restricted band of frequencies , then the expectation value for changes to . In this case can detect the noise difference only if or . is the number of frequency points in the interval .

In the case of the Kolmogorov-Smirnov (KS) estimator, the ECDF of the WOSA normalized spectrum is compared with the theoretical expectation . In the case that the noise excess is simply a constant scale factor, , over the whole frequency band, the expected CDF for the normalized WOSA spectrum is . The expected value for the KS test variable is then written as:

(4)

Once a significance value is provided, the corresponding critical value, , can be calculated from the equations for the inversion of the limiting CDF for Smirnov1939 (); Miller1956 (). The two distributions in equation (4) are incompatible at the given significance level if . This defines the detection threshold for . If is non-zero only in a restricted band of frequencies , then the expected distribution for the normalized WOSA estimator is difficult to calculate. Nevertheless the detection threshold for the KS estimator can be calculated numerically from synthetic data.

Figure 1: Non-detection ranges for and KS estimators. Noise excess is assumed as a constant multiplicative factor extending along the whole band of frequencies. Darkened intervals define the thresholds for detection, if is larger or smaller than the shaded values it can be detected with a confidence of . refers to WOSA averages. is the number of frequency points in the spectrum.
Figure 2: Non-detection ranges for and KS estimators in the case that noise excess coefficient is different from zero in a restricted band . Data are presented as a function of the ratio where is the number of frequency points in the interval and is the total number of frequency points in the spectrum. refers to WOSA averages. Confidence level for excess detection is fixed to .

In figure 1, the ranges of non-detectability of the KS and estimators are reported as a function of for three different values of the WOSA averages . Data refer to the case that the scale factor, , extends over the complete frequency band. As can be clearly seen, the estimator always has a better sensitivity than the KS estimator. The sensitivity for the case of a band limited excess noise is reported in figure 2 as a function of the ratio . As in the previous case, the estimator provides a better sensitivity with respect to the KS estimator. The difference is particularly relevant for , below such values the sensitivity of the KS estimator becomes poor. Both in figure 1 and in figure 2 the confidence level for detection is fixed at .

It is worth noting that KS algorithm can be used directly on time series to quantitatively gauging departures from a given distribution (e.g. Gaussian). Once the ECDF for the data is calculated, it can be compared with the expected distribution by calculating from equation 2. Since the distribution of coefficients is known it is straightforward to set a confidence threshold. The procedure can be applied even in presence of correlations thanks to the generalizations discussed in section V.2 and in appendix B.

ii.2 Detection of noise variations without a model

In the case that an excess of noise has to be detected by comparing different measurements the Parseval’s theorem Percival () () suggests that the sum of the elements of a sample spectrum in a given frequency band could provide a sensitive estimator for noise power variations. Such an estimator would be loosely equivalent to the discussed above, except that its statistic is hard to determine in a typical experimental situation. The statistic of at each frequency, , is defined by equation (19). Therefore, in the case of non-white noise, its parameters depend on and the statistic is different at different frequencies. Thus the statistic of is not easily known.

An interesting alternative to is provided by the KS estimator. Given two data series and of length and respectively. The hypothesis that their ECDFs have the same limiting cumulative distribution function can be tested if a distance in the ECDFs space is defined as:

(5)

with defined on the interval and Feller1948 (). Also in this case, the statistical properties of

(6)

are independent from the distributions of and . The same equations used in the case of the comparison with a given model can now be used for the inversion of the cumulative statistic of 111As formulated by equation (6), the test is searching for differences in the noise content of the data series. It does not provide information about which series has the largest noise content. Such information can be recovered graphically with a distribution plot (for example, a quantile-quantile plot). Alternatively, a pure excess noise test for a data series, , with respect to a reference series, , can be formulated by substituting with . The formalism used for the calculation of and is the same, therefore the procedures described in this paper can be applied in either case without any significant modification.. Considering the simplifying assumption of independent spectral data, the sensitivity of the KS estimator at a given significance level can be calculated in analogy to what was discussed in the previous paragraph.

Figure 3: Non detection ranges for the KS estimator for the comparison of two sample spectra. A) Noise excess coefficient on the full frequency band. is the number of WOSA averages. is the number of frequency points in the sample spectrum. B) only in a restricted interval of frequencies . is the number of points in the interval.

In figure 3 the calculated interval for non-detection is reported for the case that the excess, , is extending over the whole frequency band and for the case that in .

Iii Noise model identification

Closely related to the problem of excess noise detection, the problem of noise model identification is one of the principal scientific objectives of LPF mission. Of particular interest is the identification of a model for the force noise acting on the TMs, which can be used as a baseline for the forthcoming space-based gravitational wave observatories.

The main constrains on the identification of force noise on the TMs in LPF are fixed by i) the limited data series length, typically h; ii) the frequency range in which force noise on the TMs is dominating the signal is below mHz; iii) force noise data are not directly accessible since the system measures and reports TM displacement, requiring that force noise on TMs be reconstructed by a numerical procedure LISA7Anneke (); LISA7Luigi ().

The result of the combination of the first two constraints is that the number of segment averages in the WOSA procedure for sample spectrum estimation should be taken as low as possible so as to have a reasonable number of frequency points in the range mHz. As a consequence, the distribution of the WOSA spectrum strongly departs from a Gaussian distribution 222The distribution of the WOSA estimator for the sample spectrum is a gamma distribution, as reported in equation (23). Such a distribution tends to a Gaussian distribution as the number of averages increases. The difference between the ‘true’ distribution and the corresponding Gaussian can be quantified by the maximum distance among cumulative distribution functions, . As an example for , for , for , for and for . It is worth noting that since is a difference in the space of the values of the CDF, it can assume values in . As can be seen, the distribution of the WOSA spectrum reasonably approaches a Gaussian only for as large as ., meaning that the classical least-squares minimization procedure for parameter estimation is not well conditioned and a full maximum likelihood procedure is required.

iii.1 Maximum Likelihood Parameter Estimation

If we replace in equation (24) with a parametric model for the spectrum, , the normalized WOSA spectrum becomes parametric, , and can be used for the estimation of noise model parameters . In this section we develop the likelihood formalism for the simple case that the noise model is a function of a single parameter, . This allows us, not only to find a sensible goodness of fit estimator for that can be used also in the case of multiple parameters, but also to place the excess noise estimator, , in a more solid theoretical framework.

Indicating with the ‘true’ value for the parameter, can be rewritten as:

(7)

Here . Assuming that correctly reproduces the expected value for the spectrum, the distribution of is reported in equation (23) with and . The distribution of the samples is then:

(8)

Under the simplifying assumption that the values of are independent for different , the likelihood function can be written as:

(9)

Here is a constant term required to have a finite probability from the probability distribution function . are observed samples corresponding to . It is typically more convenient to work with the natural logarithm of the likelihood function .

(10)

is the total number of frequency samples. Taking the first derivative with respect to and equating to we find the maximum likelihood estimator for the parameter :

(11)

The value of at which corresponds to the maximum likelihood estimation for the parameter. The estimator is unbiased, since, remembering that , it can be verified that . For the practical purpose we can find the zero crossing of the reduced estimator

(12)

which is crossing zero at the same value of , since when . It is worth noting that is practically the excess noise estimator with the expectation value subtracted. It is then worth noting that can be used not only in the simple case of one parameter but it can also be applied in the general case of a model since the condition when is always satisfied.

iii.2 Kolmogorov-Smirnov Parameter estimation

As discussed in section II.1, the KS estimator can be used as an alternative to a maximum likelihood procedure for parameter estimation.

Thanks to the statistical properties of , the closer is to , the better the distribution of is described by equation (23) with and . Therefore, the Kolmogorov-Smirnov distance parameter provides an effective goodness-of-fit estimator:

(13)

is the ECDF for the current estimate, is the limiting distribution for .

KS estimation for is obtained by the minimization of with respect to . A confidence range for the parameter estimation can be readily defined from the non-rejection region of the KS-test at a given significance level. In practice, having defined a significance level, the corresponding critical values of the KS statistic, , can be calculated with standard equations Miller1956 () or Monte Carlo simulations in the case of correlated data. The values of for which provide the boundary for the confidence range at the given significance. The KS statistic can also be successfully applied to multi-parameter identification, since the convergence of to is always verified when .

Iv Analysis of non-stationary noise

The implementation of noise analysis procedures was so far discussed in the context of stationary or pseudo-stationary noise 333The term pseudo-stationary indicates a data series which is affected by slight non-stationarity of the kind that can be removed with mean subtraction or standard polynomial fit trend removal.. In the case of truly non-stationary noise the spectral content of a time series is investigated by time-frequency analysis techniques which include the spectrogram and the wavelet transform. The spectrogram is estimated by the square modulus of the short-time Fourier transform of the data Oppenheim1999 (). It provides a direct extension of the PSD formalism to non-stationary time series. Given a data series of samples, a fraction of length is windowed and the Fourier transform is applied. Then the data window is time shifted and the process is repeated. The calculation of the spectrogram is based on data windowing and the application of the Fourier transform, therefore the considerations noted in the previous sections for stationary noise can be applied directly to the spectrogram analysis of non-stationary noise.

Since the short-time Fourier transform has the same resolution across the time-frequency plane, it is often preferable to resort to the wavelet transform. Wavelet transform is a decomposition of the time series over time-frequency elements that are obtained by scaling and translating a mother function :

(14)

The function has zero average and the wavelet elements are normalized to 1. The wavelet transform of a function is then defined as:

(15)

In the discrete case, the result of the wavelet transform on a time series is an array of coefficients where is the time index and is the scale index which is associated to a given frequency band MallatWavelet (). In the case of uncorrelated Gaussian noise, the distribution of the coefficients is still Gaussian with a certain amount of correlation introduced by the convolution-like transform operation MallatWavelet (). In such a favorable situation the extension of the method described in the stationary case it appears straightforward. In particular Kolmogorov-Smirnov procedures are excellent candidates since their robustness to correlation and the possibility of extending KS distance definition to a two-dimensional space Peacock1983 (); Fasano1987 ().

V Application to LPF synthetic data

In this section the different procedures for excess noise detection and parameter estimation are applied to synthetic LPF data. This provides not only an interesting framework for testing their accuracy and precision, but it also helps clarifying the role of correlations among spectral data.

v.1 Synthetic data and noise projection

LPF provides two output channels along the principal measurement axis which are sensing the displacement of the SC relative to the free falling TM and the relative displacement between the TMs.

From the knowledge of the displacement signals and a linear model for the system dynamics, an effective force-per-unit-mass, , acting on the TMs can be extracted. is the combination of the ‘true’ force-per-unit-mass acting on TMs and a projected interferometer readout noise. Details of the calculation are reported in appendix D.

Following the same scheme, a given prediction for the different input noise sources can be projected to a prediction for the power spectrum of , as reported in figure 4. Here the model for the spectrum of force noise acting on TMs and the model for the spectrum of readout noise are projected into a model for the power spectrum of (TMs + Readout in the figure notation), which represents this paper’s baseline for the spectrum of .

Figure 4: Projection of the spectrum of . ‘Readout’ is the projection of the readout noise to . ‘TMs’ is the projection of the force noise on the test masses to . ‘TMs + Readout’ is the complete noise projection for , it represents the baseline noise level assumed in the present paper. ‘LPF worst case’ refers to a worst case scenario for and ‘LPF Spec.’ corresponds to the mission specifications.

In figure 4 we also report the project specifications for LPF and the expected noise spectrum for in a worst case scenario. In our baseline we assumed a reduced force noise on the TMs compared to the worst case but choose to keep the worst case for the readout noise. This was done in order to represent one of the possible scenarios (not the best one) that can be experienced during the mission.

The model assumed for the force noise on the TMs is characterized by a low frequency behavior and a flat part for mHz. The model can be written as . is an adjustable parameter which assumes values for the worst case scenario and for our baseline model. It is worth noting that is projected (together with the readout noise model) through LPF dynamics in order to obtain the expected noise spectrum for which we indicate with . with corresponds to ‘TMs + Readout’ in figure 4, with corresponds, instead, to the ‘LPF worst case’.

v.2 Excess noise detection

A change of the noise level on () produces a variation of the energy content of . Such variation, which may be ‘improperly’ identified as excess noise, can be detected with the procedures defined in section II. In particular we tested the detection of excess noise between two data series and between a data series and a reference model. Synthetic data were produced according to the following procedure:

  1. Different models for are produced changing around the reference value . Readout noise level is kept fixed.

  2. Corresponding noise time series for are generated using the procedure reported in FerraioliPRD2010 (). The time series are h long and have a sampling frequency of Hz.

  3. Sample spectra are calculated for each series with the WOSA algorithm. We chose a Blackman-Harris data window, segment overlap and number of segments averages .

  4. The analysis is restricted to the frequency interval mHz since, as can be seen from figure 4, it represents the region in which the force noise on the TMs dominates .

Spectral data are tested for excess noise. We used the KS algorithm (equation (6)) in the case of the test of two data series. The data series for are compared against the reference series with . In the case of the test of a data series against a model, both the KS algorithm (equation (3)) and the IR algorithm (equation (1)) are used. The reference model is the projected for . The results of the tests are summarized in table 1. KS critical values and IR confidence intervals are calculated for a significance level which corresponds to a confidence.

Each value of corresponds to a value of the in-band energy content of in the analyzed frequency band. We report in table 1 the relative change in energy corresponding to a relative change in . The quantity plays the same role of the parameter in figures 1 and 3, even though the two quantities are not completely comparable since in figure 1 assumes independence of the data. Spectral data are correlated among different frequency values because of two effects Percival (); Thomson1977 (): i) Data windowing which corresponds to a convolution in the frequency domain of the window function with the sample spectrum; ii) WOSA overlapped segment averaging. The first effect is unavoidable, the second effect, instead, can be attenuated by a proper choice of the segment overlap. It can be demonstrated 444If is the number of averaging segments, the length of each segment and the shift factor, which is indicating the number of data points two consecutive data segments are shifted by, then the correction to the variance of the averaged process is proportional to Percival () . Here are window elements shifted by a factor . With our choice of , a segment overlap of and a total data length of h, we find . that for a Blackmann-Harris window the effect is negligible with overlap.

KS vs. data KS vs. model IR
Conf. Int. MC Conf. Int.
Table 1: Detection of noise differences in the frequency band [0.1, 10] mHz. The symbol indicates compatibility between tested objects. The symbol is instead used for indicating test rejection. is an adjustable parameter which assumes value for our baseline model. Different values of correspond to different values of the in-band energy content of . We reported here relative values with respect to reference. Details on , and MC confidence interval for IR estimator can be found in appendix B and C respectively.

Since the standard statistics for the estimators ( for KS and confidence interval for IR) are calculated in section II with the assumption of data independence, the standard critical values, , for the KS estimators and the confidence intervals for the IR estimator can be applied only if the correlations among data are negligible. If this is not the case, the effective statistic of the estimators can be numerically calculated with Monte Carlo simulations. The corresponding results of a Monte Carlo simulation with realizations of the reference data series are indicated in table 1 with the suffix MC. The symbol is used to indicate that the spectral data for the corresponding is compatible with the reference. On the contrary the symbol indicates a rejection.

Observing the test results reported in columns and for the KS algorithm and the results in columns and for the IR algorithm it is readily seen that correlations among data play a role. In the case of the KS test, the comparison of the data with determines a rejection for . On the other hand, the same value is accepted when Monte Carlo result is used. In the case of the IR estimator, the comparison with the standard confidence interval leads to the rejection in correspondence with and . Such values are instead considered compatible by the Monte Carlo confidence interval. These results provide us with clear information that the presence of correlations among data has affected the tests statistics and therefore the standard equations, assuming independence among data, are not usable in this situation.

Looking at the results for the KS test between two data series, we discover that and are practically equal. In fact the results for the two corresponding columns of table 1 (columns and ) are in perfect agreement. This is the practical result of one of the most interesting properties of the KS test. The test is based on equation (6). It states that the two empirical cumulative distributions under test have the same limiting CDF. Since the sample spectra in our test are calculated following the same WOSA procedure, the degree of correlation among different frequency points is the same, therefore the test statistic is not spoiled.

Comparing and of columns and we see that the effect of correlation is to increase the maximum expected spread between ECDF and limiting CDF. Therefore the effect of data correlation is to distort the expected statistic for . The ‘distortion’ of statistic can be taken into account if an effective value for the parameter is introduced. In the case of the comparison of the ECDF for correlated data against a theoretical CDF the application of the standard values for , where , being the number of data elements, leads to a statistically unfair test. We then discovered that test fairness can be recovered if an effective value for is used rather than the standard . In particular, for spectral data produced with the WOSA method, Blackmann-Harris window, averages on overlapped segments, we obtained with for a significance level . It is worth noting that the value of is independent from the number of data points considered, and from the spectral shape, provided that the different shapes have reasonably comparable smoothness on a frequency interval comparable with the width of the first lobe of the data window. As an example, the value of is valid for LPF-like data and for white-noise equivalently. The requirement on the smoothness of the spectrum is connected to the expression of the correlations introduced by data windowing. It can be demonstrated that if the spectrum can be assumed constant in a frequency range of the order of the width of the first lobe of data window Percival () then the correlations are independent from the particular shape of the spectrum and are determined only by the window function. For such class of spectra we expect the same value for once the required significance level is fixed.

v.3 Noise model identification

The problem of noise parameter identification is strictly connected to the problem of excess noise detection by a comparison of a data series with a reference model. While in excess noise detection different data series are compared with a given reference, in parameter estimation different realizations of a parametric model are compared with a dataset in order to find the best fit. Due to this, the same algorithms (i.e. KS and IR) can be applied to the solution of the two problems. Precision and accuracy in the estimation of the parameter, , controlling the excess force noise on the TMs is tested with a Monte Carlo simulation over realizations of the same process. The data series reproduce corresponding to the reference value . Sample spectra are calculated with the procedure described above and the analysis is restricted to the frequency range mHz. For each realization, data are compared with a bank of models obtained by the projection of TMs force noise and readout noise into for different values of around the reference value. The KS and IR estimators are calculated for each element of the model bank, in particular it is simpler to analyze the results in terms of . Both KS and are expected to have a minimum corresponding to the best estimate for the parameter . The two methods proposed here are compared with the performance of a classical weighted least-squares method, which works by minimizing the mean squared error , and a Huber’s norm estimator (details are reported in appendix E). The results of the analysis are reported in figure 5, where the histograms of the best fit vales over realizations are reported for the four procedures. In the same figure we also show the evolution along the model grid of the four estimators for a particular set of data from the available .

Figure 5: KS, IR, MSE and Huber’s norm accuracy and precision in the estimation of the force-per-unit-of-mass noise parameter . Corresponding histograms report the result of MC realizations of the parameter search on a grid of template models. We also show (black dots) the evolution along the model grid of the four estimators for a particular set of data taken from the available.

The distributions for the best fit parameter are reasonably symmetric for all the estimators; mean values and sample standard deviations are respectively , , , , , , and . From the analysis of the Monte Carlo results, it is readily seen that both the KS and the estimators provide equivalently precise and accurate results. On the other hand, the MSE algorithm provides a poor estimation, both from the accuracy and from the precision point of view. The best estimation for the parameter is , which is strongly biased with respect to the reference value of . Also, the distribution of the parameter values is wider than those obtained from the KS and estimators. Huber’s norm estimator, with the chosen parameter E, performs better than MSE but the result is still far from the true value.

It is worth discussing the effect of correlations among spectral data on the estimation of the parameter . As can be observed from the results of the MC simulation, the accuracy of the estimators is not affected by correlations; the results for KS and estimators are practically indistinguishable. Some problems with can arise from data correlations when the confidence interval for a single estimation is required. As discussed in the previous paragraph, correlations modify the statistic of the IR estimator making it impossible to easily calculate confidence intervals from the standard equations. The statistic of the KS estimator is affected by correlations too, but we have demonstrated that accurate critical values can be recovered if an effective value for is used, . is a shape parameter depending only on the correlations in the spectral data and on the required significance level. It can be calculated for a specific spectrum (e.g., white noise) and effectively extended to a wide family of spectral shapes. Using the value reported in appendix B, we obtain for . The intersection of such a value with the curve of as a function of reported as black dots in figure 5A provides a confidence interval for the single estimation. In this specific case such an interval is .

Vi Conclusions

The problem of excess noise detection and noise parameter estimation for non-Gaussian data is analyzed in the framework of the LISA Pathfinder mission. Excess noise detection can be approached in two ways. In one way, the noise content of a data series is compared with a reference data series, in the other way the noise content of the data is compared with a reference model. In the first case, simple estimators like the total energy content in a data series are not suitable for formulating quantitative statements on a solid statistical basis. As an alternative a Kolmogorov-Smirnov (KS) estimator is proposed and successfully applied to LPF synthetic data. The KS estimator has the advantage of being independent of the statistical properties of the data under test. It is demonstrated in the paper that such a convenient property makes the estimator robust to correlations among spectral data. Two different estimators are investigated for the problem of comparing a data series to a model for excess noise detection. One estimator (IR) is based on the integral of the normalized WOSA spectrum, the other is a KS estimator for the comparison of an empirical cumulative distribution function with a limiting theoretical function. Despite the fact that IR estimator proves to have better sensitivity on independent data, the versatility of the KS estimator is highly advantageous on correlated data. Since the statistics of the estimators are distorted by the presence of correlations among test data, the standard and simple procedures for the determination of the confidence intervals, based on the inversion of the limiting distribution function, provide inaccurate results. While, in the case of the IR estimator, the problem can be overcome only with dedicated Monte Carlo simulations, we demonstrated that the introduction of a shape parameter allows us to use standard equations to calculate proper boundaries for KS confidence intervals. The shape parameter depends on the required significance level and on the data correlations. It does not depend on the number of data considered for the test. Since correlations among spectral data are mainly introduced by the windowing process, the shape parameter is fixed for a wide family of spectral shapes. For example, synthetic LPF data share the same shape parameter with white noise.

Closely related to noise excess detection, the problem of noise parameter identification is analyzed with a maximum likelihood approach which, in the particular case of linear dependence on a single parameter, provides an algorithm analogous to the IR estimator used for excess noise detection. A KS algorithm was proposed as an alternative to the IR algorithm and the accuracy and precision of both were tested with a Monte Carlo simulation on LPF synthetic data. Both the IR and KS estimators were demonstrated to give equivalently good results, even though the capability of the KS to handle data correlation is a clear advantage for the definition of a confidence interval for the estimated noise parameter.

Data analysis procedures introduced in this paper are easily extended to the vast context of time-frequency analysis of non-stationary noise. KS algorithm can be applied effectively both to spectrogram and wavelets coefficients provided that correlations among data are taken into account. KS algorithm has the advantage that can be generalized to two dimensions thus allowing to extend the analysis to the time-frequency plane.

Appendix A Statistic of the sample spectrum

In the case of a discrete, real-valued stationary process , the continuous spectral density function is defined as:

(16)

Here is the autocovariance sequence of the process , is the sampling period and is the frequency expressed in Hz. is defined in the range and is known as the Nyquist frequency. In the case of a finite representation of the discrete process , the approximation to the spectral density function is provided by the sample spectrum

(17)

in this case is also defined on the interval . If the sample spectrum is calculated on the grid of Fourier frequencies (, ) then it corresponds to the squared modulus of the discrete Fourier transform of the data sequence . In practical applications only the positive frequency part of the spectrum is considered, and the one sided sample spectrum is defined as with . In the rest of the paper the one sided sample spectrum will be simply named the sample spectrum.

If the data series is Gaussian distributed and the elements are independent then the Fourier transform produces a complex series whose elements are approximately independent and their real and imaginary parts are Gaussian distributed. The term is then the sum of two independent variables distributed as a where is the number of degrees-of-freedom of the distribution ( in this case). If the correlations among the elements of the data series are non-vanishing, the statistical properties of the sample spectrum, in the simplifying assumption of independent elements, can be calculated from the case of a distributed variable multiplied by a constant . The characteristic function for is

(18)

Here indicates the expected value. The inverse Fourier transform of provides the probability density function for :

(19)

This is a gamma distribution, , with and . is the gamma function. In the case of the sample spectrum at a given frequency, and 555It is easy to verify that, in the case of and , the probability density function is correctly represents the statistic of the sample spectrum. as expected for the spectrum Percival (). since in the present case.. It is useful, for the statistical analysis of the spectrum, to introduce the normalized sample spectrum

(20)

which, at each frequency , is distributed as .

The sample spectrum in equation (17) can also be seen as a special case of the power spectral density, , when the infinite data series is chopped by a square data window of length . This operation introduces a considerable amount of spectral leakage because of the convolution with the frequency response of the square data window Percival (); Harris1978 (). Therefore it is common practice to multiply the time series with a more performant data window, which increases the accuracy of the sample spectrum in the case of processes with a high dynamic range. The application of a data window introduces correlations among different elements of the sample spectrum. Such correlations affect the statistics of the spectrum, resulting in a change in the probability distribution of the sample spectrum. An analytical treatment of the spectrum statistics under such conditions is cumbersome, and it is easier to numerically evaluate (with a Monte Carlo simulation) the statistics of the sample spectrum for the case of interest.

In order to improve the variance properties of the sample spectrum, Welch’s overlapped segment averaging (WOSA) method is applied Percival (). The data series is divided in overlapping windowed segments. The estimates of the sample spectrum of each segment are then averaged. The practice of averaging overlapping segments can modify the expected statistics of the spectrum since the data in different segments can be correlated. Then, even in the simplifying assumption of vanishing spectral correlations from windowing and overlapping, the averaging process changes the statistics of the estimated sample spectrum. In the assumption of vanishing window and overlap correlations, the statistic of the WOSA spectrum corresponds to the average of gamma distributed variables

(21)

The critical function for the sum is , where is the critical function for . Thanks to equation (18)

(22)

where and . The inverse Fourier transform of provides the probability density function for the WOSA spectrum

(23)

with and . Again it is useful to define a normalized WOSA spectrum as

(24)

which is gamma distributed (equation (23)) with and .

Appendix B Kolmogorov - Smirnov Test

Kolmogorov - Smirnov is a well known test for distributions Kolmogorov1933 (); Smirnov1939 (); Feller1948 (); Miller1956 (); Fisher1983 (); Wilk1968 (). An empirical cumulative distribution (ECDF) is tested against a continuous theoretical model or, alternatively, two ECDFs are tested with the hypothesis that they share the same limiting cumulative distribution function. Indicating with the probability density function associated with a given random process , the corresponding cumulative distribution function (CDF) is defined as:

(25)

Given a particular realization of the random process :

(26)

ECDF is written as where is the number of observations which is smaller or equal to .

Given two data series and , with and not necessarily equal, we can test if the two series are two particular realizations of the same random variable by analysis of their ECDFs. Under the hypothesis that the two data series comes from the same distribution function, Kolmogorov has demonstrated that the maximum distance between the two ECDFs has a limiting distribution which is independent from the statistical properties of the corresponding random variable. If the test is performed against a theoretical distribution, the distance is defined as:

(27)

In such a case . In alternative, If the test is performed between two ECDFs, and:

(28)

The test is defined as follows:

  1. In the case of the test on a single data series, the null hypothesis is that the data are realizations of a random variable which is distributed according to the given probability distribution. In the case of two data series, the null hypothesis is that the two data series are realizations of the same random variable, which means their ECDFs should tend to the same CDF. The test rejects or accepts the null hypothesis on the basis of the analysis of .

  2. A significance level is defined as the probability that the test rejects the null hypothesis when it is indeed true.

  3. The test can be formulated in terms of critical values. The critical value for the test is the value of corresponding to the significance level. Then if , the null hypothesis is rejected.

KS theory was formulated for independent data sets and the available equations for critical values are valid only if this condition is satisfied Miller1956 (). The test can be formulated also in the presence of data correlation but the distortion to statistic introduced should be taken into account. This is possible if an effective value for is introduced as , with a shape parameter depending only on the data correlations and the required significance level. Alternatively, realistic critical values can be calculated with dedicated Monte Carlo (MC) simulations Weiss1978 (). The advantage of the shape parameter is that it depends only on correlations and the required significance level. Therefore it can be determined for a specific spectrum (e.g., white noise) and shared among a wide family of spectral shapes. Once is known it can be used to calculate critical values for correlated data using standard equations reported in the literature Miller1956 ().

Focusing on the particular problem, we performed a Monte Carlo estimation of for WOSA spectra representative of LPF. The number of frequency data considered is , corresponding values are in the case of the test against a theoretical distribution and in the case of the test between two ECDFs. The results are summarized in table 2. In the same table we report the values of that are required to obtain proper critical values from the standard equations in the case of the test against a theoretical distribution.

Table 2: Table of KS critical values for correlated spectral data. Critical values are calculated with a Monte Carlo simulation for different values of the significance level . Datasets used are representative of the data analyzed in the current paper. We reported the values for testing an ECDF against a theoretical CDF () and the values for testing two ECDFs for the same limiting CDF (). We also reported the values of the shape parameter that can be used to calculate proper critical values from standard equations in the case of correlated data. refers to the significance level whereas is the corresponding confidence level for the test. .

Appendix C IR Test

The statistics of the IR excess noise estimator are numerically estimated with a Monte Carlo simulation on the spectral data used for the present paper. Data series are h long synthetic reconstructions of the force-per-unit-mass expected on the TMs. WOSA spectra are calculated with averages on overlapped segments which were windowed with the Blackman-Harris window. Results are reported in table 3.

Table 3: Table of confidence bounds for the IR estimator on correlated spectral data. The values are calculated with a Monte Carlo simulation for different values of the significance level . Data sets used are representative of the data analyzed in the current paper. The number of frequency points is , corresponding to the number of available spectral data in the range mHz, for a time series h long and number of averages for the WOSA estimator .

Appendix D Conversion of displacement noise

LPF can be considered as a three body controlled system composed of the two TMs and the spacecraft (SC). Its equations of motion along the measurement axis can be written as:

(29)

Here:

  • and are TMs coordinates along the sensitive axis. They are relative coordinates in the SC reference frame.

  • is the absolute SC coordinate along the sensitive axis.

  • , and are the masses of the two TMs and of the SC.

  • and are the parasitic stiffnesses coupling the TMs and the SC. The TMs are coupled to the spacecraft through the parasitic stiffness thus producing an oscillator like equation of motion. The spacecraft at the same time experiences reaction forces given by and .

  • , and are the forces acting on TMs and SC respectively.

  • and are control forces on the second TM and the SC respectively. Since is an internal force to the system the SC experiences a reaction force .

  • Dots over symbols represent time derivatives.

In the main LPF science operation mode, one TM (indicated here as ) is in free-fall and provides the reference for the other TM () and the SC. In order to avoid unwanted drifting, both and the SC are controlled to follow . It is worth noting that the system is, by construction, symmetric and the role of the two TMs can be inverted. In order to avoid confusion, we indicate with the free-fall reference and with the actuated TM.

Moving to the Laplace domain, substituting for the SC dynamics and substituting for the differential coordinate of with respect to , the equation (29) can be rewritten as:

(30)

Here:

  • and are output displacement signals as provided by the interferometer readout system. is the displacement between the SC and the . is the displacement of relative to .

  • and are transfer functions of the control systems on and SC. The force applied by the controllers is calculated on the basis of the output displacement, therefore and