Xiaoran Ning

Ivan W Selesnick

Laurent Duval
email: laurent.duval@ifpen.fr

Chromatogram baseline estimation and denoising using sparsity (BEADS) $

Keywords: baseline correction, baseline drift, sparse derivative, asymmetric penalty, low-pass filtering, convex optimization

published or not. The documents may come

Chromatogram baseline estimation and denoising using sparsity (BEADS)

Introduction

Several sources of uncertainties affect the quality and the performance of gas and liquid chromatography analysis [START_REF] Smit | Noise and detection limits in signal-integrating analytical methods[END_REF][START_REF] Barwick | Sources of uncertainty in gas chromatography and high-performance liquid chromatography[END_REF]. As with many other analytical chemistry methods (including infrared or Raman spectra [START_REF] Brown | Comprehensive Chemometrics: Chemical and Biochemical Data Analysis[END_REF]), chromatogram measurements are often considered as a combination of peaks, background and noise [START_REF] Mcnulty | The effect of different baseline estimators on the limit of quantification in chromatography[END_REF]. The two latter terms are sometimes merged under different denominations: drift noise, baseline wander, or spectral continuum. For instance in [START_REF] Brown | Derivative preprocessing and optimal corrections for baseline drift in multivariate calibration[END_REF], the baseline drift is characterized as a "colored" noise, with a lowfrequency dominance in the noise power spectrum. In the following, we restrict the term "baseline" to refer to the smoothest part of the trend or bias (The portion of the chromatogram recording the detector response when only the mobile phase emerges from the column, [START_REF] Mcnaught | IUPAC. Compendium of Chemical Terminology[END_REF]), while we call "noise" the more stochastic part. Peak line shapes are of possibly various nature, from Gaussian to asymmetric empirical models [17, p. 97 sq.]. Meanwhile, they can easily be described as short-width, steep-sided up-and-down bumps. They therefore also possess relatively broad frequency spectra, albeit localized and behaving differently from the drift noise disturbance. Leaving peak artifacts (fronting and tailing, co-elution, etc.) aside, their quantitative analysis (peak area, width, height quantification) is thus hindered by the possibility to accurately remove both the smooth baseline and the random noise [START_REF] Laeven | Optimal peak area determination in the presence of noise[END_REF]. Indeed, these problems are often addressed independently, in two different steps (which could, in turn, "introduce substantial levels of correlated noise" [START_REF] Brown | Derivative preprocessing and optimal corrections for baseline drift in multivariate calibration[END_REF]): a generally low-order approximation or smoothing for the baseline, and forms of filtering for the noise on the residual chromatogram with background removed.

First, although seemingly simple, the problem of baseline subtraction remains a long-standing issue, which can be traced back to [START_REF] Wilson | The elimination of errors due to baseline drift in the measurement of peak areas in gas chromatography[END_REF][START_REF] Pearson | A general baseline-recognition and baseline-flattening algorithm[END_REF]. Recent overviews are presented in [START_REF] Schulze | Investigation of selected baseline removal techniques as candidates for automated implementation[END_REF][START_REF] Fredriksson | An objective comparison of pre-processing methods for enhancement of liquid chromatography-mass spectrometry data[END_REF][START_REF] Komsta | Comparison of several methods of chromatographic baseline removal with a new approach based on quantile regression[END_REF]. Spectral information processing [START_REF] Smit | Specification and estimation of noisy analytical signals: Part I. Characterization, time invariant filtering and signal approximation[END_REF][START_REF] Smit | Specification and estimation of noisy analytical signals: Part II. Curve fitting, optimum filtering and uncertainty determination[END_REF][START_REF] Wentzell | Signal processing in analytical chemistry[END_REF] has been a major course of action. Methods based on linear and non-linear [START_REF] Moore | Median filtering for removal of low-frequency background drift[END_REF][START_REF] Kneen | Algorithm for fitting XRF, SEM and PIXE X-ray spectra backgrounds[END_REF][START_REF] Ruckstuhl | Baseline subtraction using robust local regression estimation[END_REF] filtering, or multiscale forms of filtering with wavelet transforms [START_REF] Chau | Application of wavelet transform in processing chromatographic data[END_REF][START_REF] Hu | A background elimination method based on wavelet transform for Raman spectra[END_REF][START_REF] Cappadona | Wavelet-based method for noise characterization and rejection in high-performance liquid chromatography coupled to mass spectrometry[END_REF][START_REF] Liu | Intelligent background correction using an adaptive lifting wavelet[END_REF] have been proposed. The relative overlap between the spectra of the peaks, the baseline, and the noise has led to alternative regression Time index (sample)

Figure 1: Gas chromatogram from two-dimensional gas chromatography [START_REF] Vendeuvre | Comprehensive twodimensional gas chromatography for detailed characterisation of petroleum products[END_REF], with flat baseline and low noise: (a) analytical signal x, (b) first-order difference D 1 x and (c) second-order difference D 2 x. The chromatographic signal is sparse, as are its first-and second-order derivatives.

models, based on various constraints. The low-pass part of the baseline may be modeled by regular functions, such as low-degree polynomials [START_REF] Mazet | Background removal from spectra by designing and minimising a non-quadratic cost function[END_REF][START_REF] Zhang | Baseline correction using adaptive iteratively reweighted penalized least squares[END_REF] or (cubic) spline models [START_REF] Fischer | Background estimation in experimental spectra[END_REF][START_REF] Gulam Razul | Bayesian model selection and parameter estimation of nuclear emission spectra using RJMCMC[END_REF][START_REF] De Rooi | Mixture models for baseline estimation[END_REF], in conjunction with manual identification, polynomial fitting or iterative thresholding methods [START_REF] Gan | Baseline correction by improved iterative polynomial fitting with automatic threshold[END_REF]. Related algorithms based on signal derivatives [START_REF] Brown | Derivative preprocessing and optimal corrections for baseline drift in multivariate calibration[END_REF][START_REF] Danielsson | Matched filtering with background suppression for improved quality of base peak chromatograms and mass spectra in liquid chromatography-mass spectrometry[END_REF] have been devised. In many approaches, modeling and constraints are laid on the potential features of the baseline itself: shape, smoothness, transformed domain properties. Consequently, it appears beneficial to investigate generalized penalizations [START_REF] Eilers | A perfect smoother[END_REF][START_REF] Boelens | New background correction method for liquid chromatography with diode array detection, infrared spectroscopic detection and Raman spectroscopic detection[END_REF][START_REF] Mazet | Background removal from spectra by designing and minimising a non-quadratic cost function[END_REF][START_REF] Zhang | Baseline correction using adaptive iteratively reweighted penalized least squares[END_REF], with less stringent models on either the signal, the background or the noise. This is the very motivation of this paper: a joint estimation of these three chromatogram components while avoiding overly restrictive parametric models. Specifically, in this work, the baseline is modeled as a low-pass signal, while chromatogram peaks of interest are deemed to be sparse up to second order derivatives, leaving random noise as a residual.

In the past decade, this concept of parsimony, or sparsity, has been an active and fruitful drive in signal processing and chemistry. It entails the possibility of describing a signal of interest with a restricted number of non-zero parameters or components. Sparsity trades accuracy (of the description) with concentration (of the decomposition). Many algorithms based on sparsity have been developed for reconstruction, denoising, detection, deconvolution. Most sparse modeling techniques arose from the "least absolute shrinkage and selection operator" (better known under the lasso moniker [START_REF] Tibshirani | Regression shrinkage and selection via the Lasso[END_REF][START_REF] Rasmussen | A tutorial on the lasso approach to sparse modeling[END_REF]), basis pursuit methods [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF], total variation [START_REF] Chan | The digital TV filter and nonlinear denoising[END_REF], and compound regularization [START_REF] Bioucas-Dias | An iterative algorithm for linear inverse problems with compound regularizers[END_REF]. While the latter essentially promote sparsity, different problems require, simultaneously, other constraints, like signal smoothness or residual stochasticity.

More specifically, recent works in signal [START_REF] Mazet | Background removal from spectra by designing and minimising a non-quadratic cost function[END_REF][START_REF] Selesnick | Resonance-based signal decomposition: A new sparsity-enabled signal analysis method[END_REF][START_REF] Selesnick | Polynomial smoothing of time series with additive step discontinuities[END_REF][START_REF] Ning | ECG enhancement and QRS detection based on sparse derivatives[END_REF] and image processing [START_REF] Granai | Sparse decomposition over multi-component redundant dictionaries[END_REF][START_REF] Elad | Simultaneous cartoon and texture image inpainting using morphological component analysis (MCA)[END_REF][START_REF] Fadili | Image decomposition and separation using sparse representations: An overview[END_REF][START_REF] Starck | Image decomposition via the combination of sparse representations and a variational approach[END_REF][START_REF] Briceño-Arias | Proximal algorithms for multicomponent image processing[END_REF] have promoted a framework for decomposing potentially complex measurements into "sufficiently" distinct components. Such non-linear decompositions are termed "morphological component analysis", "geometric separation" or "clustered sparsity" [START_REF] Kutyniok | Geometric separation by single-pass alternating thresholding[END_REF]. Such approaches are amenable to analytical chemistry issues, relying on morphological properties of baselines and chromatographic peaks. Figure 1(a) displays a chromatogram x obtained from two-dimensional gas chromatography [START_REF] Vendeuvre | Comprehensive twodimensional gas chromatography for detailed characterisation of petroleum products[END_REF]. It consists of abrupt peaks returning to a relatively flat baseline, and hence exhibits a form of sparsity. Moreover, as illustrated in Fig. 1(b) and (c), the second and third-order derivatives of x are also sparse; often sparser than x itself. We thus model the peaks of a chromatogram as a sparse signal, whose first several derivatives are also sparse. In addition, baselines are sometimes approximated by polynomials or splines [START_REF] Mazet | Baseline spectrum estimation using half-quadratic minimization[END_REF][START_REF] Mazet | Background removal from spectra by designing and minimising a non-quadratic cost function[END_REF][START_REF] Zhang | Baseline correction using adaptive iteratively reweighted penalized least squares[END_REF]]. However, most baseline signals in practice do not follow polynomial laws faithfully over a long range. We thus instead model slowly varying baseline drifts as low-pass signals. The more generic low-pass model for the baseline provides a convenient and flexible way to specify the behavior of the smoothing operator, in comparison with polynomial or spline approximations.

Following aforementioned works on morphological component analysis and its variations, in particular [START_REF] Selesnick | Simultaneous low-pass filtering and total variation denoising[END_REF], we formulate an approach for the decomposition of measured chromatogram data into its modeled components: sparse peaks, low-pass baseline, and a stochastic residual. It is termed BEADS, for Baseline Estimation And Denoising with Sparsity. To this end, we pose an optimization problem wherein the terms of the objective function encapsulate the foregoing model. We develop a fast-converging iterative numerical algorithm, drawing on techniques of convex optimization. Due to its formulation as a convex problem, its solution is unique and the proposed algorithm is guaranteed to converge to the unique solution regardless of its initialization. Furthermore, we express the algorithm exclusively in terms of banded matrices, such that the proposed iterative algorithm can be implemented with high computational efficiency and low memory. Accordingly, the proposed algorithm is suitable for long data series.

Preliminaries

In this paper, lower and upper case bold are used to denote vectors and matrices, respectively, e.g., vector x and matrix A. The N point signal x is denoted as x = [x 0 , x 1 , . . . , x N -1 ]

T . The n-th element of the vector x is denoted

x n or [x] n . Element (i, j) of the matrix A is denoted A i,j or [A] i,j .
The setting of this work is in the discrete-data domain, so all derivatives will be expressed by finite differences and the words 'derivative' and 'difference' are used interchangeably. The first-order difference matrix of size (N -1) × N is defined as:

D 1 =       -1 1 -1 1 . . . . . . -1 1       (1)
and similarly, the second-order difference matrix of size (N -2) × N is defined as:

D 2 =       -1 2 -1 -1 2 -1 . . . . . . . . . -1 2 -1       . (2)
Generally, the difference operator of order k, of size (N -k) × N , is denoted as D k . For convenience, when k = 0, we define D 0 as the identity matrix, i.e., D 0 = I.

The 1 and 2 norms of x are defined as the sums:

x 1 = n |x n |, x 2 2 = n |x n | 2 . ( 3
)
To minimize a challenging cost function F , the Majorization-Minimization (MM) approach [START_REF] Figueiredo | Majorization-minimization algorithms for wavelet-based image restoration[END_REF][START_REF] Lange | Optimization transfer using surrogate objective functions[END_REF] solves a sequence of simpler minimization problems,

x (k+1) = arg min x G(x, x (k) ) (4)
where k 0 denotes the iteration counter. The MM method requires that in each iteration a convex function G(x, v) be a majorizer of F (x) and that it coincides with F (x) at x = v. That is,

G(x, v) F (x) for all x (5a) G(v, v) = F (v). ( 5b
)
With initialization x (0) and under suitable assumptions, the MM update equation (4) produces a sequence x (k) converging to the minimizer of F (x).

Baseline Estimation and Denoising: Problem Formulation

The proposed approach is based on modeling an N -point noise-free chromatogram data vector as

s = x + f , s ∈ R N . (6)
The vector x, consisting of numerous peaks, is modeled as a sparse-derivative signal (i.e., x and its first several derivatives are sparse). The vector f , representing the baseline, is a low-pass signal. We further model the observed (noisy) chromatogram data as

y = s + w (7) = x + f + w, y ∈ R N ( 8
)
where w is a stationary white Gaussian process with variance σ 2 . Our goal is to estimate the baseline, f , and peaks, x, simultaneously, from observation y. We assume that if peaks are absent, then the baseline can be approximately recovered from a noisecorrupted observation by low-pass filtering, i.e., f ≈ L(f + w) where L is a suitable low-pass filter. Hence, given an estimate x of the peaks, we may obtain an estimate f of the baseline by filtering y in ( 7) with low-pass filter L,

f = L(y -x). (9)
In this case, we may obtain an estimate ŝ by adding x,

ŝ = f + x (10) = L(y -x) + x (11) = L y + H x ( 12
)
where H is the high-pass filter,

H = I -L. (13)
In order to obtain an estimate x of the chromatogram peaks from observed data y, we will formulate an inverse problem with the quadratic data fidelity term yŝ 2 2 . Note that

y -ŝ 2 2 = y -L y -H x 2 2 (14) = H(y -x) 2 2 . (15)
Hence, the data fidelity term in the proposed formulation depends on the high-pass filter H. Also note that the data fidelity term does not depend on the baseline estimate f . The optimization problem, formulated below, produces an estimate of the chromatogram peaks x. The baseline estimate is then obtained by (9).

Low-pass and High-pass Filters

We take the filters L and H to be zero-phase, non-causal, recursive filters. In other words, they filter relatively symmetric signals (such as chromatogram peaks), without introducing shifts in peak locations. A procedure for defining such filters is given in [START_REF] Selesnick | Simultaneous low-pass filtering and total variation denoising[END_REF]. A filter is specified by two parameters: its order 2d and its cutoff frequency f c . The high-pass filter H described in [START_REF] Selesnick | Simultaneous low-pass filtering and total variation denoising[END_REF] is of the form

H = A -1 B ( 16
)
where A and B are banded convolution (Toeplitz) matrices. Expressing H in terms of banded matrices leads in [START_REF] Selesnick | Simultaneous low-pass filtering and total variation denoising[END_REF] to the development of algorithms that utilize fast solvers for banded linear systems [START_REF] Press | Numerical recipes: The Art of Scientific Computing[END_REF]Sec. 2.4]. Being convolution matrices, A and B represent linear, time-invariant (LTI) systems, and H represents a cascade of LTI systems where the LTI system B is followed by the LTI system A -1 . Using the commutative property of LTI systems, in this paper we define H as

H = BA -1 . (17)
Due to A and B being finite matrices, they are not exactly commutative. However, the difference between ( 16) and ( 17) is confined to the beginning and end of the signal and hence, for long signals, is negligible. In Sections 4.1 and 4.2, we will see that (17) serves our purposes better than [START_REF] Fadili | Image decomposition and separation using sparse representations: An overview[END_REF], allowing the derivation of a computationally efficient optimization algorithm.

Compound sparse derivative modeling

According to [START_REF] Brown | Comprehensive Chemometrics: Chemical and Biochemical Data Analysis[END_REF], the first i derivatives, i = 0, . . . , M , of the estimated peaks, x, should be sparse. In sparse signal processing, sparse signal behaviour is typically achieved through the use of suitable nonquadratic regularization terms. Therefore, to obtain an estimate x, the following optimization problem is proposed:

x = arg min x F (x) = 1 2 H(y -x) 2 2 + M i=0 λ i R i (D i x) . ( 18
)
The assumption that the observed data is corrupted by additive white Gaussian noise (AWGN) is reflected in the use of a quadratic data fitting term as is classical. The quadratic data fidelity term is given by [START_REF] Elad | Simultaneous cartoon and texture image inpainting using morphological component analysis (MCA)[END_REF]. In ( 18), D i is the order-i difference operator defined in Section 2. Functions R i : R

N -i → R, are of the form R i (v) = n φ (v n ) (19)
where the function φ : R → R, referred to as a penalty function, is designed to promote sparsity. Substituting [START_REF] Fischer | Background estimation in experimental spectra[END_REF] in [START_REF] Figueiredo | Majorization-minimization algorithms for wavelet-based image restoration[END_REF], we obtain

x = arg min x F (x) = 1 2 H(y -x) 2 2 + M i=0 λ i Ni-1 n=0 φ ([D i x] n ) ( 20
)
where N i denotes the length of D i x. The constants λ i 0 are regularization parameters. Increasing λ i has the effect of making D i x more sparse. More discussion of how to specify λ i will be given in Section 5.

Compared with [START_REF] Selesnick | Simultaneous low-pass filtering and total variation denoising[END_REF], this work introduces several modification to adapt the approach therein to chromatograms. The novel features of the approach proposed here include: 1) An M -term compound regularizer to model chromatogram peaks. 2) The modeling of the positivity of chromatogram peaks by non-symmetric penalties. 3) An improved algorithm based on MM (in contrast to ADMM as in [START_REF] Selesnick | Simultaneous low-pass filtering and total variation denoising[END_REF]) for compound regularization. The improved algorithm converges faster and does not require any user-specified step-size parameter as does the earlier algorithm of [START_REF] Selesnick | Simultaneous low-pass filtering and total variation denoising[END_REF].

Symmetric penalty functions

For many applications, samples of the target signal x and its derivative D i x are positive or negative with equal probability, or this information is not available. In such cases, the penalty function should be symmetric about x = 0. One such function is the absolute value function,

φ A (x) = |x| (21)
which leads to 1 norm regularization. While widely used, one drawback of ( 21) is that it is non-differentiable at zero, which can lead to numerical issues, depending on the utilized optimization algorithm. To address this issue, we utilize a smoothed (differentiable) approximation of the 1 penalty function, for example, the hyperbolic function As the constant > 0 approaches zero, the functions φ B and φ C approach the absolute value function. When = 0, then φ B and φ C reduce to the absolute value function. Functions φ A , φ B , and φ C are illustrated in Fig. 2 for comparison. The penalty functions and their first-order derivatives are listed in Table 1.

φ B (x) = |x| 2 + (22) or φ C (x) = |x| -log (|x| + ). (23)
In order that the smoothed penalty functions maintain the effective sparsity-promoting behavior of the original non-differentiable penalty function, should be set to a sufficiently small value. On the other hand, should be large enough so as to avoid the afore-mentioned numerical issues arising in some optimization algorithms (in particular, the MM algorithm developed below). Fortunately, we have found that the numerical issues are reliably avoided with small enough so that its impact on the optimal solution is negligible. The same holds true for the asymmetric functions discussed below. We have found that = 10 -5 works well with both φ B and φ C .

Asymmetric penalty functions

For some applications, the signal x may be known to be sparse in an asymmetric manner. For example, it is known that chromatogram data have positive peaks above a relatively flat baseline. In such cases, it is preferable to use an asymmetric penalty function that penalizes positive and negative values differently, as in [START_REF] Eilers | Unimodal smoothing[END_REF][START_REF] Mazet | Background removal from spectra by designing and minimising a non-quadratic cost function[END_REF][START_REF] Mazet | Baseline spectrum estimation using half-quadratic minimization[END_REF]. We start with the function θ : R → R defined by

θ(x; r) = x, x 0 -rx, x < 0 ( 24
)
where r > 0 is a parameter. The function θ(x; r) is convex in x and penalizes the positive and negative amplitude asymmetrically. If r = 1, then θ(x; r) = φ A (x) = |x|.

The definition [START_REF] Hu | A background elimination method based on wavelet transform for Raman spectra[END_REF] has the same drawback as [START_REF] Gan | Baseline correction by improved iterative polynomial fitting with automatic threshold[END_REF], that is, it is non-differentiable at x = 0. To alleviate this, a differentiable version of ( 24) is also proposed. In contrast to the differentiable symmetric penalties φ B and φ A , the function we propose is of the form

θ (x; r) =      x, x > f (x), |x| -rx, x < -. (25)
The function θ (x; r) differs from θ(x; r) in a way that, in the small interval [-, ], a function f (x) is defined. In Section 4.2, we will set f (x) to be a fixed second-order polynomial function and with this modification, 1) θ (x; r) will remain convex in (-∞, +∞); 2) θ (x; r) will be continuously differentiable on (-∞, +∞).

Asymmetric penalty functions are also used in algorithm backcor [START_REF] Mazet | Background removal from spectra by designing and minimising a non-quadratic cost function[END_REF][START_REF] Mazet | Baseline spectrum estimation using half-quadratic minimization[END_REF] to reflect the positivity of the peaks. We note that backcor uses non-convex penalties, while the method described here uses convex penalties (although, it can be modified to use non-convex penalties).

φ(x) φ (x) φ A (x) |x| sign(x) φ B (x) |x| 2 + x √ |x| 2 + φ C (x) |x| -log (|x| + ) x |x|+

Algorithms

Symmetric penalty functions

The first form of the proposed approach, BEADS, is defined through the minimization of the objective function F in [START_REF] Fredriksson | An objective comparison of pre-processing methods for enhancement of liquid chromatography-mass spectrometry data[END_REF], where φ is a differentiable symmetric penalty function as in Table 1. In this section, we use the MM procedure (4) to derive an iterative algorithm for this optimization. Hence, we seek a majorizer G(x, v) of F (x). First, we find a majorizer g(x, v) for φ(x) : R → R such that

g(x, v) φ(x), (26a)
g(v, v) = φ(v), for all x, v ∈ R. (26b)
Since φ(x) is symmetric, we set g(x, v) to be an even second-order polynomial,

g(x, v) = mx 2 + b (27)
where m and c are to be determined so as to satisfy [START_REF] Kneen | Algorithm for fitting XRF, SEM and PIXE X-ray spectra backgrounds[END_REF]. From ( 26), we have

g(v, v) = φ(v) and g (v, v) = φ (v), (28)
that is, g(x, v) and its derivative should agree with φ(x) at x = v. Using ( 27) and ( 28), (26) becomes

mv 2 + b = φ(v) and 2mv = φ (v). ( 29
)
Solving for m and b we obtain

m = φ (v) 2v and b = φ(v) - v 2 φ (v). (30)
Substituting [START_REF] Lange | Optimization transfer using surrogate objective functions[END_REF] in [START_REF] Komsta | Comparison of several methods of chromatographic baseline removal with a new approach based on quantile regression[END_REF], we obtain

g(x, v) = φ (v) 2v x 2 + φ(v) - v 2 φ (v) (31)
which gives

n g(x n , v n ) = n φ (v n ) 2v n x 2 n + φ(v n ) - v n 2 φ (v n ) (32) = 1 2 x T [Λ(v)] x + c(v) (33)
Ni-

1 n=0 φ(x n ) ( 34
)
where Λ(v) denotes a diagonal matrix with diagonal elements

[Λ(v)] n,n = φ (v n ) v n (35)
and c(v) is the scalar,

c(v) = n φ(v n ) - v n 2 φ (v n ) . (36)
Using [START_REF] Mazet | Baseline spectrum estimation using half-quadratic minimization[END_REF] to [START_REF] Moore | Median filtering for removal of low-frequency background drift[END_REF], we can write

M i=0 λ i Ni-1 n=0 g([D i x] n , [D i v] n ) = M i=0 λ i 2 (D i x) T [Λ(D i v)] (D i x) + c i (v) M i=0 λ i Ni-1 n=0 φ ([D i x] n ) ( 37
)
where Λ(D i v) are diagonal matrices,

[Λ(D i v)] n,n = φ ([D i v] n ) [D i v] n (38)
and c i (v) are scalars,

c i (v) = n φ([D i v] n ) - [D i v] n 2 φ ([D i v] n ) . (39)
Equality holds when x = v. Equation ( 37) implies that

G(x, v) = 1 2 H(y -x) 2 2 + M i=0 λ i 2 (D i x) T [Λ(D i v)] (D i x) + c(v) (40)
is a majorizer for F in [START_REF] Fredriksson | An objective comparison of pre-processing methods for enhancement of liquid chromatography-mass spectrometry data[END_REF]. Minimizing G(x, v) with respect to x leads to the explicit solution

x = H T H + M i=0 λ i D T i [Λ(D i v)] D i -1 H T Hy. (41)
Equation ( 41) explains why we use a differentiable penalty function. Suppose that the penalty function φ were taken to be the absolute value function, i.e., φ(x) = |x|, then (38) becomes

[Λ(D i v)] n,n = 1 [D i v] n sign([D i v] n ). ( 42
)
As the iterative algorithms progresses, some values [D i v] n go to zero due to the sparsifying property of the penalty. Consequently, 'divide-by-zero errors' are encountered in the implementation of [START_REF] Schulze | Investigation of selected baseline removal techniques as candidates for automated implementation[END_REF]. This is due to the absolute value function being non-differentiable at zero. The use of smoothed (differentiable) penalty functions such as those in Table 1, avoids this issue. For example, if we let φ = φ B , then [START_REF] Pearson | A general baseline-recognition and baseline-flattening algorithm[END_REF] becomes

[Λ(D i v)] n,n = 1 |[D i v] n | 2 + ( 43
)
and if we let φ = φ C , we obtain

[Λ(D i v)] n,n = 1 |[D i v] n | + . ( 44
)
Another issue in implementing [START_REF] Ruckstuhl | Baseline subtraction using robust local regression estimation[END_REF] resides in the computational complexity of solving the linear system represented by the matrix inverse. The computational cost increases dramatically with the length of the signal y. To address this, recall from (17) that we take the high-pass filter H to have the form H = BA -1 . Hence, (41) can be written as

x = A -T B T BA -1 + M i=0 λ i D T i [Λ(D i v)] D i -1 A -T B T BA -1 y (45) = A B T B + A T M i=0 λ i D T i [Λ(D i v)] D i A -1 B T BA -1 y (46) = AQ -1 B T BA -1 y (47)
Table 2: Algorithm to minimize cost function [START_REF] Fredriksson | An objective comparison of pre-processing methods for enhancement of liquid chromatography-mass spectrometry data[END_REF].

Input: y, A, B, λ i , i = 0, . . . , M

1. b = B T BA -1 y 2. x = y (Initialization) Repeat 3. [Λ i ] n,n = φ ([D i x] n ) [D i x] n , i = 0, . . . , M, 4. M = M i=0 λ i D T i Λ i D i 5. Q = B T B + A T MA 6. x = AQ -1 b Until converged 8. f = y -x -BA -1 (y -x) Output: x, f
where

Q = B T B + A T M i=0 λ i D T i [Λ(D i v)] D i A. (48)
Note that Q is a banded matrix. (Sums and product of banded matrices are banded.) Hence, x can be obtained with high computational efficiency and low memory using fast solvers for banded systems [39, Sec. 2.4]. Note that this can not be achieved when H is written H = A -1 B as in [START_REF] Selesnick | Simultaneous low-pass filtering and total variation denoising[END_REF].

Asymmetric and symmetric penalty functions

To account for the positivity of the chromatogram peaks, we propose a second form of BEADS which uses an asymmetric penalty that penalizes negative values of x more than positive values. Hence, in this section we use the MM approach to derive an algorithm solving the problem:

x = arg min x F (x) = 1 2 H(y -x) 2 2 + λ 0 N -1 n=0 θ (x n ; r) + M i=1 λ i Ni-1 n=0 φ ([D i x] n ) . (49)
The penalty θ (x n ; r), given by ( 25), is a differentiable version of the asymmetric penalty function θ(x n ; r). The function φ is a differentiable symmetric penalty functions such as φ B or φ C in Table 1.

To obtain a majorizer for F , we first find a majorizer for θ(x; r) : R → R, defined in [START_REF] Hu | A background elimination method based on wavelet transform for Raman spectra[END_REF]. We seek a majorizer as illustrated in Fig. 3, i.e., a function of the form

g(x, v) = ax 2 + bx + c (50)
which is an upper bound of θ(x; r) that not only coincides with θ(x; r) at x = v, but also coincides with θ(x; r) at x = s on the opposite side of zero as illustrated in Fig. 3. That is,

g(v, v) = θ(v; r), g (v, v) = θ (v; r), (51)
g(s, v) = θ(s; r), g (s, v) = θ (s; r). ( 52
)
The differentiation is with respect to the first argument. Note that a, b, c, s are all functions of v. Solving for them gives The smoothed asymmetric penalty function θ ε (x; r), r = 3

a = 1 + r 4|v| , b = 1 -r 2 , c = (1 + r)|v| 4 , s = -v. (53)
(-ε, f(-ε)) (ε, f(ε))
x Substituting [START_REF]Modern practice of gas chromatography[END_REF] in [START_REF] Tibshirani | Regression shrinkage and selection via the Lasso[END_REF], we obtain a majorizer for θ(x; r),

g(x, v) = 1 + r 4|v| x 2 + 1 -r 2 x + (1 + r)|v| 4 . (54)
Again, an issue with ( 54) is that, as v approaches zero, 'divide-by-zero' numerical errors arises. To address this issue, we define θ (x; r) to be a continuously differentiable approximation to θ(x; r). In a neighborhood of x = 0, we define θ (x; r) to be the second order polynomial [START_REF] Schottky | Über spontane Stromschwankungen in verschiedenen Elektrizitätsleitern[END_REF]

with v = , i.e., θ (x; r) =          x, x > 1+r 4 x 2 + 1-r 2 x + 1+r 4 , |x| -rx, x < - ( 55
)
where > 0 is a small constant. The new function θ (x; r), illustrated in Fig. 4, behaves similarly to θ(x; r) but is continuously differentiable. The majorizer given by ( 54) is still valid for [START_REF] Campbell | The Study of Discontinuous Phenomena[END_REF] in the domain of (-∞, ) and (-, +∞). In the domain [-, ], we use θ (x; r) itself as its majorizer. Hence, a majorizer of θ (x; r) is found to be

g 0 (x, v) =    1+r 4|v| x 2 + 1-r 2 x + |v| 1+r 4 , |v| > 1+r 4 x 2 + 1-r 2 x + 1+r 4 , |v| . (56)
A proof that g 0 (x, v) is a majorizer of θ (x; r) is given in Appendix A.

Using [START_REF] Vendeuvre | Characterisation of middle-distillates by comprehensive two-dimensional gas chromatography (GC × GC): A powerful alternative for performing various standard analysis of middle-distillates[END_REF], we then have

N -1 n=0 g 0 (x n , v n ) = x T [Γ(v)]x + b T x + c(v) (57) N -1 n=0 θ (x n ; r) (58)
where Γ(v) is a diagonal matrix with diagonal elements

[Γ(v)] n,n =    1+r 4|vn| , |v n | 1+r 4 , |v n | ( 59
)
and b is a vector with elements

[b] n = 1 -r 2 ( 60
)
and c(v) is a scalar that does not depend on x. Using ( 37) and ( 57), we find that a majorizer for F in ( 49) is given by:

G(x, v) = 1 2 H(y -x) 2 2 + λ 0 x T [Γ(v)]x + λ 0 b T x + M i=1 λ i 2 (D i x) T [Λ(D i v)] (D i x) + c(v). ( 61
)
Minimizing G(x, v) with respect to x yields

x = H T H + 2λ 0 Γ(v) + M i=1 λ i D T i [Λ(D i v)] D i -1 H T Hy -λ 0 b . (62)
Using H = BA -1 as in [START_REF] Selesnick | Simultaneous low-pass filtering and total variation denoising[END_REF], we can write (62) as

x = AQ -1 B T BA -1 y -λ 0 A T b ( 63
)
where Q is the banded matrix,

Q = B T B + A T MA, ( 64
)
and M is the banded matrix,

M = 2λ 0 Γ(v) + M i=1 λ i D T i [Λ(D i v)] D i . ( 65
)
As above, the system of equations represented by Q in (63) is banded and thus a fast solver for banded systems can be used to implement the MM update equation. Using the above equations, the MM iteration takes the form:

M (k) = 2λ 0 Γ(x (k) ) + M i=1 λ i D T i Λ(D i x (k) ) D i . ( 66
)
Q (k) = B T B + A T M (k) A (67) x (k+1) = A[Q (k) ] -1 B T BA -1 y -λ 0 A T b (68)
The complete BEADS algorithm for solving cost function ( 49) is detailed in Table 3. The run-time of BEADS, with a fixed number of 30 iterations, as implemented in MATLAB, is tabulated in Table 4. Note that applying the algorithm to a signal of length 10,000 requires less than one second of computation. Run-times were measured using MATLAB version 2010b on an i7 3.2GHz PC with 8 GB of RAM. The convergence property of the MM approach is discussed in [START_REF] Hunter | A tutorial on MM algorithms[END_REF][START_REF] Lange | Optimization transfer using surrogate objective functions[END_REF]. In particular, when the objective function is strictly convex, as is the case in ( 20) and ( 49), the MM algorithm is guaranteed to converge to the unique optimum. [START_REF] Starck | Image decomposition via the combination of sparse representations and a variational approach[END_REF]. BEADS 0.040 0.124 0.677 7.266

Input: y, r 1, A, B, λ i , i = 0, . . . , M 1. [b] n = 1 -r 2 2. d = B T BA -1 y -λ 0 A T b 3. x = y (Initialization) Repeat 4. [Γ] n,n =    1+r 4|xn| , |x n | 1+r 4 , |x n | 5. [Λ i ] n,n = φ ([D i x] n ) [D i x] n , i = 0, . . . , M, 6. M = 2λ 0 Γ + M i=1 λ i D T i Λ i D i 7. Q = B T B + A T MA 8. x = AQ -1 d Until converged 9. f = y -x -BA -1 (y -x) Output: x, f

Experiments

Baseline correction of simulated chromatograms in Gaussian noise

This example illustrates the use of BEADS for the removal of baselines in simulated chromatograms. The algorithm is compared with two other methods, airPLS [START_REF] Zhang | Baseline correction using adaptive iteratively reweighted penalized least squares[END_REF] and backcor [START_REF] Mazet | Baseline spectrum estimation using half-quadratic minimization[END_REF][START_REF] Mazet | Background removal from spectra by designing and minimising a non-quadratic cost function[END_REF]. We conduct two experiments. In both experiments, the chromatogram peaks are generated as described in [START_REF] Komsta | Comparison of several methods of chromatographic baseline removal with a new approach based on quantile regression[END_REF] (similar also peak simulation of [START_REF] Zhang | Baseline correction using adaptive iteratively reweighted penalized least squares[END_REF][START_REF] Mazet | Baseline spectrum estimation using half-quadratic minimization[END_REF][START_REF] Mazet | Background removal from spectra by designing and minimising a non-quadratic cost function[END_REF]). Namely, the signal x is created as a superposition of Gaussian functions with different amplitudes, positions and widths. For Experiments 1 and 2, we simulate the baseline in two ways, respectively:

1. Type 1 simulated baseline: as in [START_REF] Komsta | Comparison of several methods of chromatographic baseline removal with a new approach based on quantile regression[END_REF][START_REF] Mazet | Background removal from spectra by designing and minimising a non-quadratic cost function[END_REF], a baseline signal is generated as the sum of an order-p polynomial and a sinusoidal signal of frequency f . Specifically, for each realization, the order p and frequency f are uniformly distributed in a prescribed range. 2. Type 2 simulated baseline: a baseline is generated as a stationary random process with a power spectrum limited to a low-pass range of [0, f c ] Hz. Specifically, such a signal is obtained by applying a low-pass filter with cut-off frequency f c to a white Gaussian process. Sample realizations are illustrated in Figs. 5 and6 for Experiments 1 and 2, respectively. White Gaussian noise is added to each realization. We generate 500 realizations and vary the variance of the signal such that the SNR ranges from -5dB to 25dB. For each realization and SNR level, we apply the BEADS, airPLS, and backcor algorithms to estimate the baseline. The accuracy of the baseline estimation is evaluated by computing the SNR of the output, i.e, the energy of the generated baseline divided by the energy of the difference between the generated and the estimated baselines, measured in decimal.

For Type 1 simulated baselines (Experiment 1), the results are shown in Fig. 7 and Table 5. Algorithm backcor and BEADS outperform airPLS. Since the morphology of Type 1 baselines is relatively simple, backcor and BEADS perform quite similarly; however, BEADS yields a slightly smaller error on average. For Type 2 simulated baselines (Experiment 2), the results are shown in Fig. 8 and Table 6. As in Experiment 1, airPLS yields the greatest error and BEADS the least error, on average; but, the improvement of BEADS in comparison with backcor and airPLS is more significant than in Experiment 1. BEADS is better able to estimate the more challenging Type 2 baseline, because it models the baseline as a low-pass signal rather than as a parametric function (e.g., polynomial).

To clarify the use of BEADS as applied in this example: we have set the low-pass filter L = I -H to be a second order filter (d = 1 in [START_REF] Selesnick | Simultaneous low-pass filtering and total variation denoising[END_REF]) with a cut-off frequency of f c = 0.0035 cycles/sample. Although the variance of the noise can affect the choice of f c , in practice, we have the effect insignificant.

We model the chromatogram peaks as having two sparse derivatives, i.e., we set M = 2 in the cost function F in [START_REF] Starck | Image decomposition via the combination of sparse representations and a variational approach[END_REF]. We use the differentiable penalty φ B with = 10 -5 . The regularization parameters λ i , i = 0, 1, 2 should match the sparsity of x and its derivatives. For each derivative order, i: the sparser D i x, the larger the corresponding λ i should be set. For an individual signal, a simple rule is to choose λ i inversely proportional to D i x 1 for i = 0, 1, 2, and to set the proportionality constant α according to the noise variance. In this experiment, we have set λ i to be inversely proportional to the empirical mean of D i x 1 , i = 0, 1, 2, where the mean is computed over the 500 realizations. We manually tune the proportionality constant α so as to minimize the average SNR. It is also conceivable that other approaches such as cross-validation or bootstrapping can be used to select the λ i parameters.

The backcor algorithm requires two user-specified parameters: a threshold value and the order of the approximation polynomial. To set these parameters, we optimized them, via a search, to minimize the SNR. The algorithm, airPLS, also requires two user-specified parameters which we likewise set so as to minimize the SNR.

Baseline correction: Poisson observation process

With the variety of chromatographic detectors [53, p. 277-337], noise distributions may more aptly be characterized by Poisson statistics, under the denomination of shot noise [START_REF] Schottky | Über spontane Stromschwankungen in verschiedenen Elektrizitätsleitern[END_REF][START_REF] Campbell | The Study of Discontinuous Phenomena[END_REF]. Gaussian fluctuations can be regarded as a limit of a Poisson process. In modern detectors the noise can often be modeled as Poisson distributed, proportional to the square root of intensity (peaks + baseline). More specifically, if at time n we denote the Poisson observation y(n), the peak signal x(n) and the baseline f (n), then y(n) may be modeled as y(n) = (P (n)/c) 2 where c is a proportionality constant and

P (n) ∼ Poisson(λ(n)) is a Poisson random variable with mean λ(n) = c x(n) + f (n).
Although the proposed algorithm is developed under the assumption that noise is additive Gaussian, we also test its performance when the observed data follows a Poisson model. A simulated signal is shown in With this Poisson noise model, two experiments (3 and 4) are conducted following the same setup as Experiments 1 and 2 respectively. The noise level is quantified by the constant c and performance is evaluated using SNR. The results are detailed in Fig. 10 and Fig. 11. Similar to Experiment 1 and 2, BEADS uniformly outperforms the two methods to which it is compared.

Baseline correction of real chromatogram data

This example illustrates baseline correction using the proposed method, BEADS, as applied to real chromatogram data. The chromatogram data, from [START_REF] Zhang | Baseline correction using adaptive iteratively reweighted penalized least squares[END_REF], is shown in gray in Fig. 12. The algorithms, BEADS, backcor, and airPLS, are applied to estimate the baseline. The parameters for backcor and airPLS were manually set so as to obtain their best possible result. Figs. 12(a)-(c) display the estimated baseline produced by each of the three methods, and (d)-(f) show the corresponding estimated peaks (obtained by subtracting the estimated baseline from the data). The three methods are able to capture the baseline trend. However, close examination shows that BEADS exhibits less distortion in some regions. In the interval 2200-2500, backcor overestimates the baseline, while airPLS slightly underestimates the baseline.

Joint baseline correction and denoising

Some methods, e.g., backcor and airPLS, have been developed specifically for baseline removal; while other methods have been devised for the reduction of random noise. However, in many operational scenarios, both baseline drift and noise are present in the measured data. In this case, BEADS can be used to jointly perform baseline correction and noise reduction.

For illustration, white Gaussian noise is added to the chromatogram signal from [START_REF] Zhang | Baseline correction using adaptive iteratively reweighted penalized least squares[END_REF] (Fig. 13(a)). The new signal exhibits both baseline drift and additive noise. The output obtained using BEADS is illustrated in Fig. 13 Note that the estimated chromatogram peaks, illustrated in Fig. 13(b), are well delineated. The baseline is well estimated as illustrated in Fig. 13(c). The residual constitutes an estimate of the noise. The cost function history of the iterative algorithm, shown in Fig. 14, indicates that the algorithm converges in about 20 iterations.

Conclusion

This paper addresses the problems of chromatogram baseline correction and noise reduction. The approach, denoted 'BEADS,' is based on formulating a convex optimization problem designed to encapsulate non-parametric models of the baseline and the chromatogram peaks. Specifically, the baseline is modeled as a low-pass signal and the series of chromatogram peaks is modeled as sparse and as having sparse derivatives. Moreover, to account for the positivity of chromatogram peaks, both asymmetric and symmetric penalty functions are utilized (symmetric ones for the derivatives). In order that the iterative optimization algorithm be computationally efficient, use minimal memory, and can be applied to very long data series, we formulate the problem in terms of banded matrices. As such, the algorithm leverages the availability of fast solvers for banded linear systems and the majority of the computational work of the algorithm resides therein. Furthermore, due to the formulation of the problem as a convex problem and the properties of the majorization-minimization (MM) approach by which the iterative algorithm is derived, the proposed algorithm is guaranteed to converge regardless of its initialization.

The baseline correction performance is evaluated on simulated chromatogram data and compared with two state-of-art methods. In particular, the proposed method, BEADS, outperforms methods based on polynomial modeling when low-order polynomials are not efficient representations of the baseline. Experiments suggest that BEADS can better estimate the baselines of real chromatograms. Finally, since BEADS jointly estimates the baseline and the chromatogram peaks, it can be used to perform both baseline correction and noise reduction simultaneously.

Figure 2 :

 2 Figure 2: Penalty functions φ A (x), φ B (x), and φ C (x).

Figure 3 :

 3 Figure 3: Asymmetric penalty function and its majorizer. (v = 0.8, r = 3)

Figure 4 :

 4 Figure 4: Continuously differentiable asymmetric penalty function θ (x; r) in (55). The function is a second-order polynomial on [-, ].

Figure 5 :Figure 6 :Figure 7 :Figure 8 :

 5678 Figure 5: Simulated chromatograms with Type 1 baseline [27].

Figure 9 :

 9 Figure 9: An example of chromatogram data (peaks + baseline) in Poisson noise.

Fig. 9 (

 9 Fig. 9(a); Fig. 9(b) and (c) show the Poisson observation and the difference between the observation and the noise-free simulated signal (peaks + baseline).With this Poisson noise model, two experiments (3 and 4) are conducted following the same setup as Experiments 1 and 2 respectively. The noise level is quantified by the constant c and performance is evaluated using SNR. The results are detailed in Fig.10and Fig.11. Similar to Experiment 1 and 2, BEADS uniformly outperforms the two methods to which it is compared.

  (b)-(d). (We have used BEADS parameters r = 6, f c = 0.006 cycles/sample, d = 1 and φ = φ C .)

Figure 10 :

 10 Figure 10: Experiment 3. Output SNR as a function of C.

Figure 11 :

 11 Figure 11: Experiment 4. Output SNR as a function of C.

Figure 14 :

 14 Figure 14: Cost function history of the BEADS algorithm.

Table 1 :

 1 Symmetric penalty functions and their derivatives.

Table 3 :

 3 BEADS algorithm to minimize cost function

Table 4 :

 4 Run-time (in sec.) of BEADS for N -point data.

	N	10 2	10 3	10 4	10 5

Table 5 :

 5 Experiment 1. The mean and standard deviation (std) of SNR. The table shows the result when input SNR is 0dB, 10dB and 20dB.

		0 dB		10 dB		20 dB
		mean	std	mean	std	mean	std
	BEADS	28.1 8.52	32.64	8.02	38.33 6.74
	backcor 24.91 9.75	31.27	8.33	36.47 6.53
	airLPS 20.26 9.65	22.54 10.15	26.71 7.76

Table 6 :

 6 Experiment 2. The mean and standard deviation (std) of SNR. The table shows the result when input SNR is 0dB, 10dB and 20dB.

	0 dB		10 dB		20 dB
	mean	std	mean	std	mean	std
	BEADS 18.75 3.71	19.99 3.17	20.89 3.32
	backcor 17.20 4.57	18.93 3.74	19.54 3.18
	airLPS 16.71 4.80	17.52 5.54	17.98 4.82

$ This research is supported by the NSF under grant CCF-1018020.

The proposed method also provides a more general framework for different types of signals. By adjusting regularization parameters λ i and penalty functions, the model can be customized for signals of different kinds. For example, an electrocardiography (ECG) signal is only sparse in its first several derivatives (the signal itself is not sparse) [START_REF] Ning | ECG enhancement and QRS detection based on sparse derivatives[END_REF], therefore, by setting λ 0 = 0 and choosing proper λ i , i = 1, 2, 3 and using symmetric penalty functions, BEADS can be used for ECG baseline estimation and denoising.

BEADS uses symmetric/asymmetric convex penalty function to promote positivity of the estimated peaks, however, other penalty functions, such as non-convex penalties, may provide further improvements. In addition, due to the growing interest on hyphenated techniques in analytical chemistry [START_REF] Vendeuvre | Comparison of conventional gas chromatography and comprehensive two-dimensional gas chromatography for the detailed analysis of petrochemical samples[END_REF][START_REF] Vendeuvre | Characterisation of middle-distillates by comprehensive two-dimensional gas chromatography (GC × GC): A powerful alternative for performing various standard analysis of middle-distillates[END_REF], with an increasing need laid on repeatability, extensions to two-dimensional chromatography are advisable. These will all be considered in future work.

Appendix A.

We prove that g 0 (x, v) in ( 56) is a majorizer of θ (x; r). Define

For θ (x; r) to be a majorizer of g 0 (x, v), we need to show that If v > , then

which can be written as

where we have used v > and |x| . If v < -, then

which can be written as

where we have used v < -and |x| .


Amelia Carolina Sparavigna

Recurrence plots of sunspots and solar radio flux

Keywords: Recurrence plots, Solar physics, Sunspots, Solar radio flux

The paper shows the recurrence plots of two time series, concerning data of the solar activity. Data ranging from 1947 to 2015 are about the sunspot number and the values of solar radio flux, which are known as highly correlated. Recurrence plots clearly display the oscillating behaviour of these quantities. The corresponding plots are also compared, evidencing some difference in the behaviour of the two time series.

Introduction Several physical processes and dynamical systems have recurrent behaviours that can be analysed using methods based on recurrence plots. It was in 1987 that Eckmann, Kamphorst and Ruelle proposed these plots, showing that they were able of visualising the recurrence of states of dynamical systems with arbitrary number of dimensions [START_REF] Eckmann | Recurrence plots of dynamical systems[END_REF][START_REF] Eckmann | Ergodic theory of chaos and strange attractors[END_REF]. In fact, a recurrence plot visualises the trajectory of a given system in its phase space through a twodimensional representation of its recurrences. At first, such recurrence of a state at time i at a different time j, was depicted within a two-dimensional squared matrix (i,j) with black and white dots, where black dots marked a recurrence. The matrix had a diagonal black line, called the line of identity, simply meaning that each point is recurrent with itself. Besides black and white plots, today software exists which gives coloured images [START_REF] Kononov | Visual Recurrence Analysis[END_REF], software which is easy to apply to any time series, as proposed in [START_REF] Sparavigna | Carbon dioxide concentration and emissions in atmosphere: Trends and recurrence plots[END_REF][START_REF] Sparavigna | Recurrence plots of geolocated time series from satellite maps of NOAA STAR Vegetation Health Index[END_REF][START_REF] Sparavigna | Recurrence plots from altimetry data of some lakes in Africa[END_REF]. Cross recurrence plots can also be used for comparing time series [START_REF] Marwan | Cross recurrence plots and their applications[END_REF][START_REF] Zbilut | Embeddings and delays as derived from quantification of recurrence plots[END_REF]. In a previous paper [START_REF] Sparavigna | Recurrence plots of sunspots, solar flux and irradiance[END_REF], we have discussed data concerning sunspots, solar flux and irradiance, between 1955 and 2007 using the recurrence plots with black and white images; here we will see the plots from 1947 to 2015, obtained from a software giving coloured images. Number of sunspots are from Ref.10 and data of radio flux are from the Ref.11. By a simple image processing we can compare the behaviours of sunspot and solar flux. The comparison allows evidencing a difference in the behaviour of the two time series.

Sunspots and solar radio data

As told in [START_REF] Vv | [END_REF], humans have observed sunspots for a very long time. The first written record of sunspots was made by Chinese astronomers around 800 B.C. In 1128, an English monk, John of Worcester, made the first drawing of sunspots. However, it was only in 1843 that Samuel Schwabe discovered the cycle of the average number of sunspots. The average duration of the sunspot cycle is 11.1 years. In this paper we will use data of sunspots from [START_REF]NASA's Marshall Space Flight Center[END_REF], but we can find them from other sources too. For instance, in [START_REF] Sparavigna | Recurrence plots of sunspots, solar flux and irradiance[END_REF], we used data on sunspots from the National Geophysical Data Center (NGDC), NOAA Satellite and Information Service [13]. Let us remember that a sunspot is a region on the photosphere with a lower temperature than its surroundings and having an intense magnetic activity. Data show that the number of sunspots oscillates; other data showing an oscillation are the solar radio emission and the irradiance data. The measurement of radio waves is used to monitor the structure of the solar corona, the outer most regions of the Sun's atmosphere. Solar radio emissions at different frequencies allow observing the radiation from different heights in the solar atmosphere. A lower frequency corresponds to a higher height of origin, because the frequency decreases uniformly outwards: 245 MHz originates high in the corona, while 15,400 MHz originates in the low corona. The 5 MHz emission corresponds to about 10 solar radii height [START_REF] Mclean | Solar radiophysics: studies of emission from the sun at metre wavelengths[END_REF]. The microwave wavelength 2800 MHz daily radio flux highly correlates with the daily sunspot number [START_REF]The numerical data for graphs and selected bibliography are given[END_REF]. The two databases are usually considered as interchangeably. Moreover, in 1995, Schmahl and Kundu found that the solar radio fluxes in the spectral range 1000-9400 MHz correlate well with the total solar irradiance [START_REF] Schmahl | Microwave proxies for sunspot blocking and total irradiance[END_REF]. Solar Radio 10.7 cm Daily Flux (2800 MHz) noon values had been obtained from the National Research Council of Canada since 1947. Until 1991 the observations were made at the Algonquin Radio Observatory, near Ottawa [START_REF]The data of the solar radio flux[END_REF]. Over 1990-1991 the research program was transferred to the Dominion Radio Astrophysical Observatory, near Penticton, British Columbia. The characteristics of the observations are surveyed in [START_REF] Covington | Solar Radio Emission at 10.7 cm[END_REF]. The accuracy of flux was discussed in Ref. [START_REF] Tapping | Limits to the Accuracy of the 10.7cm Flux[END_REF]. The sources and emission mechanisms contributing to the 10.7 cm flux are proposed in Ref. [START_REF] Tapping | Recent Solar Radio Astronomy at Centimeter Wavelengths: The Temporal Variability of the 10.7-cm Flux[END_REF]. In [START_REF] Sparavigna | Recurrence plots of sunspots, solar flux and irradiance[END_REF], because we had the data on sunspot number reported with monthly averaged values, we had to evaluate the same average for the 10.7 cm flux data. In this paper, we simply use data from Ref.11, which are already monthly averaged, to have the recurrence plot. So we will see the plots from data [START_REF]NASA's Marshall Space Flight Center[END_REF] and [START_REF]Government of Canada, Solar radio flux, Archive of measurements[END_REF] and compare them. Moreover, in [START_REF] Sparavigna | Recurrence plots of sunspots, solar flux and irradiance[END_REF], we simply used black and white recurrence plots RP for time-discrete variables, as given by: ( )

j i j , i x x RP - - ε Θ = (1)
where i,j are the discrete time variables, corresponding to months. i x is the variable (the sunspots number or the 10.7 cm flux). The Θ-function depends on the difference between the distance of the two points and on a predefined cut-off distance ε. We had two-dimensional squared matrices with black and white dots, where black dots marked the recurrence. The Figure 1 is giving the results from Ref.9, from 1955 to 2007. On the left, it is shown the recurrence plot (RP) for the sunspots, and on the right the RP for the radio flux. Image processing can be useful too. In the larger panel of Fig. 1, we show an image obtained combining the single RPs of sunspots and radio flux. In green, they are shown the points where the two RPs are coincident and in red, the points where they are different. Two marked regions, a and b, corresponds to periods where we find a high solar activity. In a (cycle 20), the recurrences are almost equal, but in b, (cycle 23) the red pixels evidence the difference in time when the maximum of activity is reached. This means that the 10.7 cm radio signal and the sunspot number are highly correlated but not interchangeable, especially near solar maximum, as reported in Ref. [START_REF] Wilson | 10.7 cm solar radio flux and magnetic complexity of active regions[END_REF]. The speculative possibility proposed in Ref. [START_REF] Wilson | 10.7 cm solar radio flux and magnetic complexity of active regions[END_REF], is that the radio-flux lag discriminates between long-period and short-period cycles, so it is an indicator of the presence of two solar cycle modes. Further analyses of sunspots and radio flux had been given in [START_REF] Johnson | Power law relating 10.7 cm flux to sunspot number[END_REF][START_REF] Ghosh | Multi-technique Analysis of the Solar 10.7 cm Radio Flux Time-Series in Relation to Predictability[END_REF][START_REF] Fragkou | Spatiotemporal time series analysis methods for the study of turbulent magnetohydrodynamic channel flows[END_REF].

Data from 1947 to 2014

As previously told, in the present paper we will use data ranging from 1947 (February) to 2015 (May). So we have a number of 819 discrete times. To have the recurrence plots we use software developed by E. Kononov.

In the Figure 2 we can see the images reproducing two layouts of the recurrence plot of the sunspot numbers. On the left, we have a two-dimensional square matrix with coloured pixels, on the right a greyscale images. Both axes of the images are the months from February 1947 to May 2015. In the colour image, the red pixels are marking a recurrence, whereas, in the greyscale image, are the white pixels to mark the recurrence. The Figure 3 is showing the same for the solar radio flux. To compare the recurrence plots of Figures 2 and3, we can use the following approach (used in [START_REF] Sparavigna | Cloud tracking in METEOSAT images using GIMP program[END_REF][START_REF] Sparavigna | The GNU Image Manipulation Program applied to study the sand dunes[END_REF]), the result of which is given in the Figure 4. Let us consider the greyscale image on the right of Figure 2. We can easily invert the grey tones by using GIMP software, for instance. Then, this image with inverted tones is added to the greyscale image on the right of Figure 3. The result is given in the Figure 4. In it, we can see that darker and brighter lines exist, which are the lines where we have a great difference of the two plots. The darkest line corresponds to year 2001 and then corresponds to region b in the Figure 1 [9]. We can do the same with the coloured images having the result in the Figure 5, which is in agreement with that of Figure 4. As concluded in [START_REF] Sparavigna | Recurrence plots of sunspots, solar flux and irradiance[END_REF], the data of sunspots and radio flux show clearly an oscillating behaviour, both in graphs and in recurrence plots. However, the recurrence plots are able to highlight some features, which are not easily seen in the graphs, because of their strong dispersion of values, dispersion which is present also in monthly averaged data. In [START_REF] Sparavigna | Recurrence plots of sunspots, solar flux and irradiance[END_REF], from the recurrence plots, we have seen that, for cycle number 23, there is a difference in time when the maximum of activity is reached by the sunspots and by the solar radio flux. This is in agreement with the observation of Ref.21: radio signal and sunspot number are highly correlated but not interchangeable, especially near solar maximum. Here again, in the Figures 4 and5 we can see this feature. Let us conclude remarking the speculative possibility, proposed in Ref. [START_REF] Wilson | 10.7 cm solar radio flux and magnetic complexity of active regions[END_REF], of a presence of two solar cycle modes. 2 and3, we can use the following approach. Let us consider the greyscale image on the right in Figure 2. We can invert the grey tones by using GIMP software for instance. Then, the image with inverted grey tones is added to the grey tone image on the right of Figure 3. The resulting image is given here on the left. In it, we can see some dark and bright lines: these lines mean that there, we have a great difference of the two plots. The darkest line corresponds to year 2001 (cycle number 23) and then corresponds to region b in the Figure 1. On the right we have the same image with adjusted contrast and brightness.

Figure 1 :

 1 Figure 1: The image shows the two-dimensional square matrix with black and white dots, where black dots mark a recurrence for data of sunspot numbers (on the left) and solar flux (on the right). Both axes are time axes. The recurrence is within a threshold of 10 percent. Note the almost equal behaviour. Data are from the first of January 1955 till December 31, 2007, monthly averaged values of sunspot number and 10.7 cm flux data. In the larger panel, we show an image obtained combining the two single recurrence plots: in green, the points where the RPs are coincident and in red, the points where they are different. The regions a (cycle 20) and b (cycle 23) correspond to two maxima of the solar activity. In a, the recurrences are almost equal, but in b, the red pixels evidence a difference in time lag. Note that b is centred about 2001.

Figure 2 :

 2 Figure 2:The image shows two layouts of the recurrence plot of SUNSPOT numbers we can obtain using the software by E. Kononov. We have two-dimensional squared matrices with coloured and grey-tone pixels. Both axes are the numbers of months from February 1947 to May 2015, which are 819. On the left, the red pixels are marking a recurrence, whereas, on the right, are the white pixels to mark the recurrence.

Figure 3 :

 3 Figure 3:The image shows two layouts of the recurrence plot of SOLAR FLUX we can obtain using the software by E. Kononov. We have two-dimensional squared matrices with coloured and grey-tone pixels. As in the Figure2, time series are ranging from February 1947 to May 2015. On the left, the red pixels are marking a recurrence, whereas, on the right, are the white pixels to mark the recurrence.

Figure 4 :

 4 Figure 4: To compare the recurrence plots of Figures2 and 3, we can use the following approach. Let us consider the greyscale image on the right in Figure2. We can invert the grey tones by using GIMP software for instance. Then, the image with inverted grey tones is added to the grey tone image on the right of Figure3. The resulting image is given here on the left. In it, we can see some dark and bright lines: these lines mean that there, we have a great difference of the two plots. The darkest line corresponds to year 2001 (cycle number 23) and then corresponds to region b in the Figure1. On the right we have the same image with adjusted contrast and brightness.

Figure 5 :

 5 Figure 5: We can do the same as in the Figure4for the coloured images of Figures2 and 3. Let us consider the coloured image on the left in Figure2. We can invert the colour tones. Then, the image with inverted colour tones is added to the colour image on the right of Figure3. The resulting image is given here on the left. On the right we have the same image with adjusted contrast and brightness. The result is in agreement with that of Figure4.

Acknowledgements

Many thanks for data collections of solar activity. We acknowledge the receipt of the dataset (composite_d25_07_0310a.dat) from PMOD/WRC, Davos, Switzerland, unpublished data from the VIRGO Experiment on the cooperative ESA/NASA Mission SOHO. Many thanks to the site Natural Resources Canada, of the Government of Canada. Many thanks to the Solar Physics Group at NASA's Marshall Space Flight Center