TFAW: wavelet-based signal reconstruction to reduce photometric noise in time-domain surveys

TL;DR

TFAW is a wavelet-enhanced version of the Trend Filtering Algorithm that significantly improves the detection and characterization of astrophysical signals in noisy time-domain survey data.

Contribution

TFAW introduces wavelet-based noise filtering into TFA, enhancing signal detection and parameter estimation without altering intrinsic signal characteristics.

Findings

Increases signal detection efficiency up to 2.5 times for low SNR data.

Improves transit detection rate by 2-5 times in low-SNR regimes.

Reduces noise standard deviation by approximately 40 times compared to TFA.

Abstract

There have been many efforts to correct systematic effects in astronomical light curves to improve the detection and characterization of planetary transits and astrophysical variability. Algorithms like the Trend Filtering Algorithm (TFA) use simultaneously-observed stars to remove systematic effects, and binning is used to reduce high-frequency random noise. We present TFAW, a wavelet-based modified version of TFA. TFAW aims to increase the periodic signal detection and to return a detrended and denoised signal without modifying its intrinsic characteristics. We modify TFA's frequency analysis step adding a Stationary Wavelet Transform filter to perform an initial noise and outlier removal and increase the detection of variable signals. A wavelet filter is added to TFA's signal reconstruction to perform an adaptive characterization of the noise- and trend-free signal and the noise contribution at each iteration while preserving astrophysical signals. We carried out tests over simulated sinusoidal and transit-like signals to assess the effectiveness of the method and applied TFAW to real light curves from TFRM. We also studied TFAW's application to simulated multiperiodic signals, improving their characterization. TFAW improves the signal detection rate by increasing the signal detection efficiency (SDE) up to a factor ~2.5x for low SNR light curves. For simulated transits, the transit detection rate improves by a factor ~2-5x in the low-SNR regime compared to TFA. TFAW signal approximation performs up to a factor ~2x better than bin averaging for planetary transits. The standard deviations of simulated and real TFAW light curves are ~40x better than TFA. TFAW yields better MCMC posterior distributions and returns lower uncertainties, less biased transit parameters and narrower (~10x) credibility intervals for simulated transits. We present a newly-discovered variable star from TFRM.

Tables

Table 1. Table 1 : Planet parameters used for TFAW simulations

Planet	$R_P\left(R_J\right)$	$M_P\left(M_J\right)$	$M_{\ast }\left(M_{\odot }\right)$	$R_{\ast }\left(R_{\odot }\right)$	P (d)
1	1.98	1.40	1.36	1.23	0.8468
2	1.23	1.10	1.36	1.23	0.4842

Table 2. Table 2 : Mutually exclusive detections and mean SDE values for simulated planetary transits 1 and 2 as shown in Table  1 . N TFA subscript N TFA \mathrm{N}_{\mathrm{TFA}} : not detected using TFAW light curves, but detected using TFA data. N TFAW subscript N TFAW \mathrm{N}_{\mathrm{TFAW}} : detected using TFAW light curves, but not detected using TFA data. N mut subscript N mut \mathrm{N}_{\mathrm{mut}} : simultaneous detections with TFA and TFAW. SDE TFA subscript SDE TFA \mathrm{SDE}_{\mathrm{TFA}} : mean TFA SDE. SDE TFAW subscript SDE TFAW \mathrm{SDE}_{\mathrm{TFAW}} : mean TFAW SDE. Percentage values in parenthesis are with respect to the 5000 tested transits.

Planet	$N_{\mathrm{TFA}}$	$N_{\mathrm{TFAW}}$	$N_{\mathrm{mut}}$	${\mathrm{SDE}}_{\mathrm{TFA}}$	${\mathrm{SDE}}_{\mathrm{TFAW}}$
1	265	784	561	16.99	17.86
	(5.3$\%$)	(15.7$\%$)	(11.2$\%$)
2	79	503	612	14.61	15.92
	(1.58$\%$)	(10.1$\%$)	(12.2$\%$)

Table 3. Table 3 : Detection distributions for simulated planetary transits 1 and 2 as shown in Table  1 for three bins of total noise contributions.

$\sigma$				$0.01-0.04$			$>0.04$
	$N_{\mathrm{TFA}}$	$N_{\mathrm{TFAW}}$	$N_{\mathrm{mut}}$	$N_{\mathrm{TFA}}$	$N_{\mathrm{TFAW}}$	$N_{\mathrm{mut}}$	$N_{\mathrm{TFA}}$	$N_{\mathrm{TFAW}}$	$N_{\mathrm{mut}}$
Planet 1	-	-	75	76	131	318	189	653	167
	$N_{\mathrm{TFA}}$	$N_{\mathrm{TFAW}}$	$N_{\mathrm{mut}}$	$N_{\mathrm{TFA}}$	$N_{\mathrm{TFAW}}$	$N_{\mathrm{mut}}$	$N_{\mathrm{TFA}}$	$N_{\mathrm{TFAW}}$	$N_{\mathrm{mut}}$
Planet 2	10	41	107	67	393	491	2	69	14

Table 4. Table 4 : Noise parameters for the simulated light curves

Signal type	SNR	${\sigma }_w$ (mag)	${\sigma }_r$ (mag)
Sinusoidal	High	0.01	0.005
Sinusoidal	Low	0.1	0.005
Transit	High	0.01	0.005
Transit	Low	0.04	0.005

Table 5. Table 5 : Top table: Actual parameters values used to simulate the transit of Planet 1. Posterior transit parameters values and their uncertainties (with the 25% and 75% quantile as the upper and lower errors) for TFA and TFAW MCMC fits. Bottom table: 95% confidence intervals of the highest probability density for Planet 1 transit parameters TFA and TFAW MCMC fits.

Parameters	a ($R_{\ast }$)	q	i (∘)	per (days)	p	l
Simul. values	3.392862	0.214	88	0.8468	0.16542286	0.312
TFA MCMC	3.53363${}_{-0.08075}{}^{}{}_{}{}^{+0.06460}$	0.37952${}_{-0.09784}{}^{}{}_{}{}^{+0.10736}$	87.7895${}_{-0.28489}{}^{}{}_{}{}^{+0.19976}$	0.84628${}_{-0.00008}{}^{}{}_{}{}^{+0.0009}$	0.19564${}_{-0.0032}{}^{}{}_{}{}^{+0.00347}$	0.34398${}_{-0.08413}{}^{}{}_{}{}^{+0.09770}$
TFAW MCMC	3.39438${}_{-0.00228}{}^{}{}_{}{}^{+0.01072}$	0.21317${}_{-0.00591}{}^{}{}_{}{}^{+0.00206}$	88.0075${}_{-0.00212}{}^{}{}_{}{}^{+0.00687}$	0.84687${}_{-0.0001}{}^{}{}_{}{}^{+0.00008}$	0.16656${}_{-0.00102}{}^{}{}_{}{}^{+0.00312}$	0.31259${}_{-0.00189}{}^{}{}_{}{}^{+0.00603}$
	95% confidence intervals of the highest posterior density
TFA MCMC	3.368 - 3.672	0.186 - 0.577	87.31 - 88.21	0.8461 - 0.8465	0.186 - 0.204	0.161 - 0.564
TFAW MCMC	3.384 - 3.432	0.198 - 0.219	87.99 - 88.02	0.8467 - 0.8470	0.164 - 0.172	0.305 - 0.335

Equations

$$F(i)=\sum _{j=1}^M{c}_j{X}_j(i),$$

$$D=\sum _{i=1}^N[Y(i)-A(i)-F(i){]}^2,$$

$$\hat{Y}(i)=Y(i)-\sum _{j=1}^M{c}_j{X}_j(i),$$

$$P(k)=\sum _t\mid W_k(t){\mid }^2$$

$$maxscale=\lfloor \frac{\log (\frac{n}{l-1})}{\log 2}\rfloor ,$$

$$\widehat{Y^{\prime }}(i)=\hat{Y}(i)-ISWT(\widehat{Y(i)},\text{noise level}),$$

Full text

*Aims. *We present TFAW, a wavelet-based modified version of TFA. TFAW aims first, to increase the periodic signal detection and second, to return a detrended and de-noised signal without modifying its intrinsic characteristics.

*Methods. *We modify TFA’s frequency analysis step adding a Stationary Wavelet Transform filter to perform an initial noise and outlier removal and increase the detection of variable signals. A wavelet-based filter is added to TFA’s signal reconstruction to perform an adaptive characterization of the noise- and trend-free signal and the underlying noise contribution at each iteration while preserving astrophysical signals. We carried out tests over simulated sinusoidal and transit-like signals to assess the effectiveness of the method and applied TFAW to real light curves from TFRM. We also studied TFAW’s application to simulated multiperiodic signals.

*Results. *TFAW improves the signal detection rate by increasing the signal detection efficiency (SDE) up to a factor $\sim$2.5$\times$ for low SNR light curves. For simulated transits, the transit detection rate improves by a factor $\sim$2-5$\times$ in the low-SNR regime compared to TFA. TFAW signal approximation performs up to a factor $\sim$2$\times$ better than bin averaging for planetary transits. The standard deviations of simulated and real TFAW light curves are $\sim$40$\%$ better compared to TFA. TFAW yields better MCMC posterior distributions and returns lower uncertainties, less biased transit parameters and narrower ( $\sim$10$\times$) credibility intervals for simulated transits. TFAW is also able to improve the characterization of multiperiodic signals. We present a newly-discovered variable star from TFRM.

1 Introduction

The photometric precision and accuracy achieved by an astronomical survey is a key factor in transit signal detection. The photometric precision that can be attained depends on several factors such as the atmosphere conditions, the instrument, the number of non-distorted stellar PSFs that can be imaged on the CCD, the reduction pipeline and the photometric errors of the measurements. Detectors with wide field of view (FoV) often suffer from additional issues such as undersampled PSFs, atmospheric extinction, source blending and optical distortion at the edges of the detector. These and other effects increase the noise of the light curves and also introduce systematic variations that decrease the detection probability of any periodic signal in the data.

Many of the systematic variations in a given light curve are shared by light curves of other stars in the same data set. In order to remove those systematics, one can identify the objects in the field that suffer from the same kind of variations as the target (correlated noise) and then build and apply a filter based on the light curves of these comparison stars. Trend Filtering Algorithm (hereafter TFA) (Kovács et al. 2005) and SysRem (Tamuz et al. 2005) are often applied to remove systematic variations in time-domain surveys, in particular wide FoV and/or multi-telescope transiting exoplanetary surveys, such as SuperWASP (Pollacco et al. 2006), MEarth (Irwin et al. 2009), HAT-South (Bakos et al. 2008), TFRM-PSES (Fors et al. 2013) or NGTS (Wheatley et al. 2013).

Wavelets have been applied in several astronomy areas like astronomical signal processing (Starck & Murtagh 1994; Núñez & Otazu 1996; Starck et al. 1998); redshift spectra study (Machado et al. 2013); cosmic microwave background (Moudden et al. 2005), baryon acoustic oscillation analysis (Arnalte-Mur et al. 2012); solar activity (Aschwanden et al. 1998; Giménez de Castro et al. 2001); stellar activity, pulsation and rotation analysis (de Freitas et al. 2010; Bravo et al. 2014); signal detection (Szatmary et al. 1994; Fors et al. 2008; Otazu et al. 2002; Régulo et al. 2007); galaxy distribution morphology (Antoja et al. 2012) and light curve noise analysis (Cubillos et al. 2017), and filtering (Carter & Winn 2009; Grziwa et al. 2014, 2016; Grziwa & Pätzold 2016; Waldmann 2014).

TFA implements correlated noise removal using multiple stars in the field. The first part of the algorithm consists of removing trends from trend- and noise-dominated time series in order to increase the probability of detecting weak signals (called frequency analysis in the original article). Secondly, if a significant periodic signal is present in the target light curve, TFA reconstructs the shape of that signal using an iterative process (called signal reconstruction). One key aspect of the signal reconstruction step is that it should minimize the chances of suppressing the signal or modifying its intrinsic astrophysical characteristics. Several efforts have extended TFA to improve the signal recovery, such as using a meticulous selection of the template stars (for example, Kim et al. (2009) through hierarchical methods or Principal Component Analysis, (Petigura & Marcy 2012), or taking into account stellar variability or multiperiodic signals (Kovacs & Bakos 2008)). TFA is further described in Sect 2.1.

The TFA modification we present here (hereafter TFAW) introduces a wavelet filter in order to 1) remove outliers using a wavelet-inferred estimation of the signal, 2) search for periods in this outlier-free and de-noised signal during the frequency analysis step, 3) estimate the shape of the trend- and noise-free phase folded signal and 4) iteratively de-noise the trend-free light curve during the signal reconstruction process. TFAW combines TFA detrending and systematics removal capabilities with the wavelet transform’s signal decoupling and de-noising potential. The filter is built using the Stationary Wavelet Transform algorithm (hereafter SWT) or à trous algorithm (Holschneider et al. 1989). TFAW differs from other wavelet-based noise-filtering algorithms in that it does not require any parametric model fitting as in Carter & Winn (2009) nor any extra computational method (Waldmann 2014). Also, the noise contribution of the signal is estimated directly from its SWT at each iteration step and the de-noising is done through the subtraction of this contribution from the signal. This allows TFAW to de-noise the signal without modifying any of its intrinsic properties in contrast to wavelet coefficient thresholding (Grziwa et al. 2014, 2016) that can lead to distortions of the signal and introduce artificial oscillations (ripples) around discontinuities (Mallat 2008). The aim of this article is to compare the signal detection efficiency, the de-noising and signal reconstruction capabilities of TFAW and TFA, by assessing its performance on a set of simulated and real light curves. Both planet detection and characterization conducted in current and near-future surveys, such as CoRoT (Auvergne et al. 2009) and Kepler (Borucki et al. 2003), could be direct candidates to benefit from TFAW improvements in signal detection and signal reconstruction without distortion shown in this work.

In Sects 2.2 and 2.3 the SWT and the Wavelet Power Spectrum are defined. In Sects  2.4 and 2.5, the TFAW and a SWT-based outliers removal algorithm are described, respectively. In Sects  2.6 and  2.7 the selection criteria for the signal and noise levels and the mother wavelet are introduced.

In Sect 3, we apply TFA and TFAW to sinusoidal and transit-like simulated light curves affected by combinations of Gaussian and correlated noises. To assess the performance of both algorithms we compare them in a number of tests. First, in Sect 3.1, the TFA and TFAW planetary transit detection efficiency is evaluated. Second, the intrinsic signal shape and depth invariance after TFAW detrending and de-noising is validated in Sect 3.2. Next, the application of the algorithm to multiperiodic superimposed signals is studied in Sect. 3.3. In Sect. 3.4, the better signal approximation of TFAW with respect of TFA bin averaging is shown. The performance of TFA and TFAW is assessed in terms of the bias of the fitted transit parameters values and their uncertainties in Sect. 3.5.

We apply TFAW to light curves from a wide FoV ground-based data set: the Telescope Fabra-ROA at Montsec (TFRM) in Sect 4.1, to confirm that TFAW outperforms TFA. One of the light curves detrended with TFAW from TFRM is found to be newly-discovered as variable star. In Sect. 4.2, a TFAW vs. TFA quantitative assessment is performed for a set of real light curves from TFRM.

2 The TFAW algorithm

2.1 TFA basics and formulation

The key point behind TFA is the assumption that many stars within a photometric data set suffer from the same or similar systematic effects. Taking a representative set of template stars, one can build a filter function that the removal of correlated systematics from a target light curve. Given a target N-point light curve and a set of M zero-averaged template stars {$X_{j}(i)$;$i=1,2,...,N$;$j=1,2,...,M$}, the filter for that target is defined as:

[TABLE]

The set of coefficients {$c_{j}$} is determined by minimizing the following expression:

[TABLE]

Where $\{Y(i)\}$ represents the target light curve and $\{A(i)\}$ is either constant in the case of the frequency analysis step of the algorithm or, the best representation of the trend- and noise-free signal to be found in the iterative signal reconstruction step. Once the filter has been computed, the corrected light curve is defined by:

[TABLE]

Although TFA is able to remove trends and systematics from photometric time series, this algorithm does not completely decouple stochastic and correlated noises from the underlying signals.

2.2 The Discrete Wavelet Transform and the Stationary Wavelet Transform

Unlike the Fourier transform, the wavelet transform decomposes a given time series into its wavelets, i.e., highly localized impulses obtained from scaling and shifting the mother wavelet function. These scaling and shifting operations allow us to calculate the wavelet coefficients, which represent the correlation between the wavelet and a localized section of the signal. The wavelet coefficients are calculated for each wavelet segment, giving a time-series function measuring the wavelets’ correlation to the signal. In comparison to the sine wave used in the Fourier transform, which is smooth and of infinite length, the wavelet is irregular in shape and compactly supported. These properties make wavelets an ideal tool for analyzing signals of a non-stationary nature. Their irregular shape allows us to analyze signals with discontinuities or sharp changes, while their compactly supported nature allows temporal localization of the signal’s features.

In the Discrete Wavelet Transform (DWT) the signal is decomposed into dyadic blocks (i.e. shifting and scaling is based on a power 2). The DWT of a given signal is calculated by applying a series of filters to the signal. In the first level of decomposition the signal passes through a low pass filter resulting in a convolution of the original signal and the filter. Simultaneously, the signal is also decomposed using a high-pass filter. The outputs give us the detail coefficients (related to the high-pass filter) and the approximation coefficients (related to the low-pass filter). The filter outputs are then sub-sampled by 2. In an iterative process, the original signal is then successively decomposed into components of lower scale, while the high frequency components are retained but not further decomposed. The decomposition is represented as a binary tree with the nodes representing a sub-space with a different time-frequency localization. The tree is known as a filter bank while the wavelet multi-scale decomposition and reconstruction is known as multi-resolution analysis (MRA) (Mallat 1989). The maximum number of decomposition levels that can be performed is dependent on the input size of the data: a signal of length $2^{k}$ can be decomposed into k discrete levels or scales.

The DWT is not a time-invariant transform and is therefore very sensitive to the alignment of the signal in time (the DWT of a translated version of a signal might not be the translated version of the DWT of the signal). The Stationary Wavelet Transform (SWT) is a wavelet transform algorithm designed to overcome the lack of translation-invariance of the DWT. As with the DWT, high and low pass filters are applied to the signal at each level. However, in order to achieve translation-invariance, the downsamplers and upsamplers are removed and the filter coefficients are upsampled by a factor of $2^{(k-1)}$ in the kth level of the algorithm. By padding the filters at each level with zeroes, the two new sequences at each level have the same length as the original time series.

An example of a MRA using the SWT is shown in Fig. 1. Our test signal is a combination of two sinusoidal signals, one of high frequency and another of low frequency. They both have the same amplitude and are affected by Gaussian random noise. As can be seen in Fig. 1, SWT is able to decouple each of these signal contributions at different decomposition levels. This is shown by computing the inverse SWT (hereafter ISWT) at each level. High-frequency noise has been separated from the sinusoidal signals, appearing as an extra contribution to the lower decomposition levels. In this example, noise can be fairly represented by the combination of the first two decomposition levels. The high frequency sinusoidal signal can be recovered by summing the third and fourth decomposition levels while the low frequency one would be recovered through the combination of the remaining levels. Crucially, the wavelet decomposition allows the signal characteristics to change through the time series, while retrieving both time and frequency information from the original signal. Other techniques such as the discrete Fourier Transform convert data from time into frequency domain but in doing so, time information is lost.

2.3 Wavelet Power Spectrum

The wavelet transform converts signal $f$ into a series of wavelet coefficients $W_{k}(t)$. Each of these coefficients represents the amplitude of the wavelet function at a particular location within the signal at a particular wavelet scale. A useful way to determine the distribution of energy within the signal is to compute the wavelet power spectrum (WPS). The local WPS at a particular decomposition level $k$ is calculated by summing up the squares of wavelet coefficients for that level (Torrence & Compo 1998):

[TABLE]

In the case of the Continuous Wavelet Transform (CWT), one can easily compute a relationship between the scale and the corresponding period. Following the method of Meyers et al. (1993), the relationship between the equivalent Fourier period and the wavelet scale can be derived analytically for a particular wavelet function, $\psi$, by substituting a cosine wave of a known frequency into the convolution of the mother wavelet with the Discrete Fourier Transform (DFT) (Torrence & Compo 1998). The use of the CWT power spectrum as a way to detect periodicities and study other astrophysical properties such as pulsation and rotation has already been studied by de Freitas et al. (2010) and Bravo et al. (2014).

When using the SWT to compute the WPS, the relationship between the scale and the frequency cannot be so easily computed as the signal energy can be distributed among different scales. However, the fact that the SWT power spectrum also presents peaks, will be of use for the outlier removal and signal detection methods described in Sect. 2.5 and Sect. 2.6.

2.4 TFAW: application of the SWT

In the original version of TFA, once a significant periodic signal has been found in the target time series during the frequency analysis step, the phase-folded light curve is used to iteratively estimate $\{A(i)\}$ (the noise-free signal). They employ the simple bin average method, using a fixed number of bins (100 for $\sim$3000 data points) to ensure statistical stability and to obtain a reasonable noise averaging. They point out, however, that more accurate methods can be used to approximate $\{A(i)\}$. For example, in the case of planetary transits, one could use a Mandel & Agol (2002) model fit to get a more precise estimation of the noise- and trend-free signal.

In TFAW we follow a similar approach to the original TFA. First, an initial filter is computed using $\{A(i)\}$ = $\langle{Y}\rangle$. Then, the SWT and WPS of the filtered data $\hat{Y}(i)$ = $Y(i)-\{A(i)\}$ are computed. Using the method explained in Sect. 2.5 outliers are removed and an estimated signal is obtained. The latter is then used to run Box Least Square (BLS) (Kovács et al. 2002) (or Lomb-Scargle (LS) Scargle (1982)) to search for any periodic signal in our target light curve.

If any significant period is found, the light curve is phase folded to that period. Instead of using a bin average, $\{A(i)\}$ is estimated using the SWT of the phase folded signal. We decompose the phase folded light curve up to a SWT decomposition level given by the following equation:

[TABLE]

where n is the signal length and $l$ is the wavelet filter length.

At each decomposition level, we set all the detail coefficients and approximation coefficients to zero except the ones corresponding to that level. Then, the ISWT is computed for each level separately (as seen in Fig. 1). In order to estimate the shape of the noise- and trend-free signal we set a threshold decomposition level or signal level (see details in Section 2.6). $\{A(i)\}$ is then computed as the sum of the ISWTs of those levels from signal level to the last decomposition level. In this way we ensure that our signal’s estimate is separated from the high frequency noise while preserving the high-frequency shape of the signal (unlike normal binning; see Sect. 3.4 for details in SWT signal approximation versus bin averaging): the high frequency noise components are better characterized by the lower decomposition levels as they have higher frequency resolution.

This new signal estimate, $\{A(i)\}$, can then be used to compute a new set of {$c_{j}$} coefficients by means of equation 2. The new filter obtained from this set of coefficients gives us the new corrected light curve at each step of iteration through equation 3. We set an additional filter at each iteration given by the SWT of the phase folded light curve. In this case, instead of using those decomposition levels above signal level to reconstruct the noiseless signal, we set a noise level so that the noise is characterized by the sum of the ISWTs of those levels below it. Thus equation 3 at each iteration step becomes:

[TABLE]

The iteration of the algorithm stops when the relative difference between the standard deviations of the residuals in the successive iterations falls under a certain limit. As with the original TFA algorithm, we set the limit at $10^{-3}$.

To summarize, the main steps of TFAW are the following:

An initial filter is computed using $\{A(i)\}$ = $\langle{Y}\rangle$ as with the original TFA to remove trends and other systematics. 2. 2.

The SWT and WPS are computed from the filtered light curve obtained in 1. 3. 3.

Outliers are removed and an outlier-free and de-noised estimated signal is obtained. 4. 4.

The estimated signal is used to search for periodicities. 5. 5.

If a significant period is found, the signal is phased folded to run the iterative signal reconstruction. Otherwise, no reconstruction is performed. 6. 6.

Using the SWT a new signal estimation $\{A(i)\}$ is computed by means of the signal level. The noise contribution of the light curve is estimated using the noise level. 7. 7.

The new $\{A(i)\}$ and noise contribution are used to compute the new filtered signal $\hat{Y^{\prime}}(i)$. 8. 8.

Iteration continues until the convergence criterion is met.

2.5 Outlier Removal Using the SWT

An outlier can be defined as an observation point that is distant from the norm of the sample. Their origin can be diverse, but the presence of outliers in a data sample leads to effects in the mean which can display a bias toward the outlier value. Also, depending on their nature, they may impact the time series analysis with respect to modeling, estimation or forecasting. There are several methods for outlier detection: model-based  (Peirce 1852; Grubbs 1950); interquartile range (Tukey 1949); distance-based (Knorr et al. 2000) and wavelet-based (Mallat & Hwang 1992; Bilen & Huzurbazar 2002; Grané & Veiga 2010).

One way of detecting outliers is through wavelet thresholding (Bilen & Huzurbazar 2002). First, the wavelet transform of the target signal is computed. Then, a threshold value is computed at each decomposition level. Any wavelet coefficient that exceeds the threshold corresponding to that scale can be considered a potential outlier. Once the outliers have been detected, their associated coefficients are set to zero and the inverse wavelet transform is computed. However, this approach can introduce artificial oscillations (ripples) around outliers and discontinuities (Mallat 2008; Grziwa et al. 2014, 2016) and also, if the threshold is not correctly computed, modifications to the signal’s shape can result.

The outlier detection and removal approach developed for this work is a distance-based method that relies on an estimated signal. As stated in Sect. 2.4, a noise-free signal can be obtained using the sum of the ISWT for those levels from the signal level to the last decomposition level. Instead of using a pre-set signal level to obtain our estimated signal, we can define it using the WPS defined in Section 2.3. We know that the WPS determines the distribution of energy within the signal. In the case of the SWT, it allows to determine the decomposition levels at which the correlation between the signal and the wavelet transform at that scale was higher. Thus, the scale of the highest peak in the WPS defines our signal level. However, given that the nature of our light curves could be multi-periodic resulting in peaks at different scales, the peak at the lowest scale (i.e. highest frequency) is used to define the signal level. In addition, given the discrete nature of the SWT and, also depending on the number of decomposition levels, the energy of our signal can be distributed within more than one scale. Thus, the signal level is always set to also consider the previous scale in order to avoid losing as less signal contributions as possible. With this signal level we can then reconstruct the estimated signal, as explained before, by summing the corresponding ISWTs. A threshold can then be built in such a way that, given a point in the time series, if its distance to the estimated signal is above the threshold, then, it is considered an outlier and removed (or replaced by any given value). In this work, the threshold value is computed using the Universal Threshold (Donoho & Johnstone 1994) and, every point exceeding 5 times this value, are considered as outliers. Once the target light curve is free of outliers, the estimated signal is recomputed. Except in the case of a very high frequency signal (for which the WPS peak at the lowest scale could coincide with the scale characterizing the noise), this final estimated signal represents an outlier-free and de-noised approximation to our target signal that, as stated before, will be used to search for any periodicities within our light curve.

Fig. 2 shows an example of the outlier removal capabilities of the method. Left plot corresponds to a simulated phase folded light curve of Planet 1 in Table 1. Using the signal estimation explained above, once the threshold value has been estimated, the algorithm is able to effectively remove any outlier as seen in Fig. 2 right.

2.6 Signal and noise level choice criteria

The selection of the signal level and noise level parameters must be carefully performed. As described above, high-frequency noise appears as an extra contribution to the lower decomposition levels of the SWT. However, high frequency signals hidden within the data (i.e. second order sinusoidal modulations, etc…) could be present in some of the levels in which the noise also appears. One way of solving this issue would be to increase the number of data points to a higher power of two. The greater the number of levels we decompose our signal, the easier it is to separate the noise contribution from the high frequency signals. Nonetheless, even if part of the signal is mixed with the noise at some level, using TFAW signal reconstruction capabilities, one could be able to recover the signal with almost any distortion given that the correct period has been found except in those cases in which the signal has exactly the same (or very close) frequency as the noise.

One way to diminish the signal contribution lost at each iteration step and avoid as much signal distortion is to select a noise level value that uses just the lowest SWT decomposition levels (usually the first one), and accept the noise contribution from the lower decomposition levels. We recommend setting a noise level equal to 1 (lowest SWT decomposition level) for time series shorter than 2048 points (i.e. 10 decomposition levels) and increase it with care for each extra power of two of the time series length.

Regarding the signal level for the phase folded light curve, given that it should provide the best signal approximation, we recommend using the method explained in Sect. 2.5. First computing the WPS of the phase folded signal (which should be already outlier-free) and then estimating the signal level from the WPS’ peaks.

2.7 Mother wavelet selection criteria

There are several mother wavelets with different analytic properties that can be used for signal decomposition. For our application of the SWT-modified version of TFA we use the non-orthogonal, symmetric base of the biorthogonal (bior) wavelet family, first proposed by Cohen et al. (1992). More specifically, we used the pre-computed values of the bior 3.9 wavelet included in the PyWavelets111https://pywavelets.readthedocs.io module. We selected this wavelet family for TFAW as it allows the construction of symmetrical wavelet functions and because the shape of the reconstruction scaling function is very similar to the characteristic shape of a planetary transit (searches for other types of astrophysical variability may benefit from other wavelet shapes).

3 TFAW performance

In Kovács et al. (2005) several tests are presented demonstrating the signal detection and reconstruction capabilities of TFA. In this section we explore how the inclusion of the SWT filter in TFAW can improve the detection of variable signals, assess the impact that this filter has on the frequency power spectra, study the effects on the SNR of the light curves. Also we compare the use of the SWT to compute the signal approximation against bin averaging, test the likeness of the final TFAW light curve to the inserted variable signal and compare the results with the ones obtained using the original TFA. For all the simulated light curves hereafter and the real ones in Section 4.1 we have used 2048 data points and set the noise and signal levels to the lower decomposition level and the one given by the WPS peak as explained in Sect. 2.5 and 2.6.

3.1 TFAW vs TFA transit detection efficiency

In order to assess the transit detection capabilities of BLS during TFAW and TFA frequency analysis step, we generate a set of 5000 time series with different random combinations of Gaussian and correlated noises. We also introduce linear and exponential trends (i.e. to simulate changes in airmass, changes in the CCD position of the object, etc), as well as gaps and jumps in the data. We use a set of 250 simulated template stars suffering from the same jumps, gaps and trends as the target time series, affected by random distributions of $\sigma_{\mathrm{w}}$ and $\sigma_{\mathrm{r}}$ combinations. To each of these light curves, we add the transit signal given by the parameters of planets 1 and 2 in Table 1 separately. We use batman222https://www.cfa.harvard.edu/~lkreidberg/batman Kreidberg (2015) package to simulate a Mandel & Agol (2002) planetary transit model. The total noise contribution for each light curve ranges from 0.005 to 0.2 mag for planet 1 and from 0.005 to 0.1 mag for planet 2.

We want to compare the detection rates obtained with TFA and TFAW light curves. We use the same two criteria as Kovács et al. (2005) and define a detection as the following:

• •

The highest peak in the BLS power spectrum must have a frequency between $[\mathrm{f}_{\mathrm{p}}-0.001,\mathrm{f}_{\mathrm{p}}+0.001]$, where $\mathrm{f}_{\mathrm{p}}$ is the frequency of the simulated planetary transit.

• •

The signal detection efficiency (SDE) (Kovács et al. 2002) of the highest peak in the power spectrum must be greater than 6.

We run BLS using 100 bins and 99000 frequency steps to ensure statistical stability  (Kovács et al. 2002). The transit search is done in the (0.01, 12) days range. We count the number of times the signal is detected in the TFAW light curves but not detected in the TFA ones. By performing this and the opposite test using a template sample of 250 stars, we obtain the results shown in Table 2.

The example shown in Table 2 shows that the number of non-simultaneous detections (both in absolute and percentage) is always higher (a factor $\sim$3$\times$ for Planet 1 and $\sim$6$\times$ for Planet 2) in the case of TFAW light curves for both planetary transits. Also the mean SDE values are higher for TFAW than for TFA. As can be seen in Fig. 3 the SDE improvement for a given light curve can be up to a factor $\sim$2.5$\times$ for low SNR signals. Table 3 shows the distribution of the TFA and TFAW detections in three total noise contribution bins. In the high-SNR regime both TFA and TFAW behave in a similar way. TFAW performs better than TFA in detecting the transit signal in the mid-SNR regime ($\sigma$=0.01-0.04) by a factor $\sim$1.16$\times$ for Planet 1 and a factor $\sim$1.6$\times$ for Planet 2. In the low-SNR case ($\sigma>$0.04), TFAW detects the transit in $\sim$2.3$\times$ more light curves than TFA for Planet 1 and in $\sim$5.2$\times$ for Planet 2. Taking into account the whole noise range, TFAW improves the detection rate by a factor $\sim$1.6$\times$ for both planetary transits. Figure 4 shows two examples of the increase in the SDE of the BLS power spectrum peaks for Planet 1 in Table 1 in the low SNR regime. Top two panels correspond to a simultaneous detection for which TFAW yields an increased SDE for the true period of 11.2 versus 6.2 for TFA. The bottom two panels present the power spectrum for a TFAW mutually exclusive detection. In this case, the peak corresponding to the true period of the transit is hidden within the noise peaks of the TFA power spectrum with a SDE of $\sim$3.9 (below the detection threshold defined above). For TFAW, the peak can easily be identified and has a SDE of $\sim$8.3 (a factor $\sim$2.1$\times$ higher than for TFA). In addition, TFAW is able to detect the planetary signals, with a SDE above 6, up to a total noise of $\sim$0.19 mag for planet 1 and up to $\sim$0.07 mag for planet 2 whereas TFA is able to detect them up to $\sim$0.1 and $\sim$0.045 mag respectively.

3.2 TFAW vs TFA signal reconstruction

To illustrate the noise and trend filtering efficiency of the TFAW method during the iterative signal reconstruction step, we generate a sinusoidal signal affected by a combination of different levels of simulated Gaussian, $\sigma_{\mathrm{w}}$, and correlated, $\sigma_{\mathrm{r}}$ noises. As in Sect. 3.1 we also introduce linear and exponential trends, as well as gaps and jumps in the data. The period of the simulated signal is 0.63 days and 0.03 mag amplitude. Again, we use a template of 250 simulated template stars suffering from the same trends and systematics as the target star, each of them affected by a random distribution of $\sigma_{\mathrm{w}}$ and $\sigma_{\mathrm{r}}$ combinations (see Table 4 for details).

Fig. 5 shows the noise filtering capabilities of TFAW compared to TFA for a simulated target star with a sinusoidal modulation. While the noise dispersion is clearly diminished, there is no modification to either the amplitude or the phase of the signal, and the sampling timescale of the signal remains unchanged (unlike in simple binning). We also show the LS power spectrum obtained in the frequency analysis step. As can be seen the height and amplitude of the peaks for the TFAW case have a higher power value compared to the TFA ones.

In Fig. 6, the same sinusoidal signal is simulated but with lower SNR. As before, the noise dispersion is also significantly diminished while the period and amplitude of the signal remains unchanged. In this case, the LS power spectrum is improved by a greater factor than the one for the high SNR case. However, in this case, although the amplitude of the power spectrum peaks is greater in the case of TFAW, the SDE values for both TFA and TFAW are very similar (10.55 and 10.6, respectively) due to the low LS continuum.

Fig. 7 simulates a planetary transit. We use batman, with the parameters of Planet 1 in Table 1 to simulate a Mandel & Agol (2002) planetary transit model, and the noise parameters in Table 4. The high-frequency noise is filtered by TFAW during the iterative reconstruction step, resulting in a better defined transit without modifying its shape and depth. This result overcomes the artificial ripples and transit depth modification introduced in Fig. 6 of Grziwa & Pätzold (2016) due to wavelet coefficient thresholding.

The TFAW BLS power spectrum for this high SNR planetary transit remains almost unchanged with respect to the TFA version (with a lower noise floor), because the BLS power spectrum is not primarily limited by high-frequency random noise. The SDE values for this case are also very similar; 23.6 for TFA and 24.5 for TFAW.

Fig. 8 shows the same transit but with a much larger noise component (see Table 4). As in the case of the simulated sinusoidal with low SNR (see Fig. 6), TFAW improves the characterization of the signal’s shape while diminishing the noise, recovers the correct period and transit depth, and improves the detectability of the transit signal in the BLS periodogram (with and SDE of 17.8 for TFAW and 13.3 for TFA).

3.3 Application to multiperiodic signals

Kovacs & Bakos (2008) presents a variation of TFA to extend its application to multiperiodic signals using their Fourier representation (for the cases in which such representation is adequate). However, they state that if the signal has additional components such as transients or planetary transits, a more complicated model should be used to approximate each of those extra signals present in the light curve. TFAW can be used to separate the different signal contributions directly from the signal after TFA’s frequency analysis step (provided that the periods of each of the signals have been previously and correctly determined using BLS or another method). As a first example, we simulate the same high SNR-type sinusoidal signal as in Fig. 5 and Fig. 6 modulated with a 0.1438-day sinusoid with an amplitude of 0.006 mag, and with noise parameters according the first row of Table 4. As can be seen in Fig. 9, if we are able to find the correct frequency of the secondary signal (notice the peak around 6.9 d*-1* in the corresponding LS periodogram), we can fully recover it by applying TFAW directly (i.e. with no need of subtracting the primary signal) to the raw data without modifying the amplitude and shape of any of the underlying signals.

The second example in Fig. 10 shows the results of applying TFAW to a high SNR-type light curve affected by two planetary transits simulated with batman using the parameters of planets 1 and 2 in Table 1, and noise parameters according the second-to-last row of Table 4. Again, if we are able to find the period of the secondary transit after the frequency analysis step, we can separate the secondary transit from the primary simply by phase folding the raw light curve and applying TFAW. As in the other examples, we are also able to improve the SNR of both the primary and secondary signals.

3.4 Wavelet versus bin average signal approximation

In order to test the likeness of the signal approximation, $\{A(i)\}$, to the inserted signal TFAW compared to the bin averaging used by TFA, we use Planet 1 in Table 1. We inject it in 500 light curves, using batman and a Mandel & Agol (2002) planetary transit model, for increasing values of noise contribution (i.e. lower transit depth). For each simulation, we compute the signal estimation given by the SWT and bin averaging with 100 bins and obtain their deviation with respect to the simulated input signal. Fig. 11 shows that, in general, the estimation of the transit shape given by the sum of the ISWTs defined by signal level provides a better representation than the one given by the TFA bin average method. This is especially true for low $\sigma_{\mathrm{signal}}$, with up to a factor of $\sim$2$\times$ of improvement.

In the case of planetary transits, as can be seen in Fig. 12, the improvement in the signal approximation is due to the fact that the wavelet reconstruction of the signal better fits the ingress and egress profiles of the transit even in the cases of low signal-to-noise ratios compared to the bin average method.

3.5 Comparison of TFA and TFAW transit parameters fit values and uncertainties

We want to compare TFA vs. TFAW performance in terms of assessing the bias of the fitted transit parameters values and their uncertainties. To quantify those, a low SNR-type planetary transit of the Planet 1 according Table 1 and last row of Table 4 was considered. This was modeled with batman, following a Mandel & Agol (2002) analytic transit model. We used the Markov Chain Monte Carlo (MCMC) sampler provided by emcee333http://dfm.io/emcee (Foreman-Mackey et al. 2013) to sample the posterior distribution of the 6 transit parameters (a, q, i, per, p, l). Keeping the eccentricity fixed, we consider a uniform distribution for the priors and run the sampler with 200 chains and 5000 iterations with a burn-in phase of 1000 iterations.

Fig. 13 shows the 1-D and 2-D projections of the posterior probability distributions of the 6 MCMC fitted parameters for the TFA and TFAW detrended light curves. Similarly, at top panel of Table 5 we compare the injected transit parameters values with the ones obtained through MCMC for TFA and TFAW posterior probability distributions. MCMC parameter values correspond to the 50% quantile while the uncertainties are computed from the 25% and 75% quantiles as the upper and lower errors. In the case of TFA, while some of the distributions are fairly behaved, for parameters a, q and p, they present wider features either characterized by larger uncertainties (a, q) and/or larger biases (a, p). On the other hand, for TFAW light curves, the biases with respect to the initial parameters values are strongly diminished compared to TFA ones. As for the uncertainties, they are, in general, smaller for TFAW than for TFA. In some parameters, such as i, l and a, they are largely decreased. At the bottom panel of Table 5, the 95% confidence highest probability density credibility intervals for both TFA and TFAW MCMC cases, are shown. The width of the credibility intervals is $\sim$10$\times$ narrower (except for the period which is better determined) in the case of TFAW compared to the TFA ones.

4 Application to Real Survey Light Curves

In this section we test the TFAW algorithm on real data from the wide FoV ground-based TFRM-PSES survey.

4.1 The Telescope Fabra-ROA at Montsec

The Telescope Fabra-ROA at Montsec (TFRM) is a refurbished f/0.96 0.50m Baker-Nunn Camera (Henize 1957) for robotic CCD wide-FoV surveying purposes (Fors et al. 2013). The project is a joint collaboration between the Reial Acadèmia de Ciències i Arts de Barcelona - Observatori Fabra and the Real Instituto y Observatorio de la Armada. The telescope is installed at the summit of the Montsec d’Ares (Spain) as part of the Observatori Astronòmic del Montsec.

The TFRM is conducting a targeted M-dwarfs transiting exoplanetary survey (TFRM-PSES) since its deployment at the end of 2011. TFRM-PSES benefits from the wide FoV of 19.4 sq. deg. for a 16 MPix camera, to deliver at 25 sec-cadence few-mmag precision-level light curves for every object brighter than R=14 mag (del Ser et al. 2015; Fors et al. 2013).

4.1.1 Data description

The data from the TFRM-PSES survey used for the TFAW performance assessment comprises 2048 data points from 30 nights observed during 2013, 2014 and 2015 for a field centered at $(\alpha,\leavevmode\nobreak\ \delta)$ = (10:14:44, +48:30:00). All light curves in the data set have been observed with a 475 nm cutoff frequency glass filter (Schott GG475), 17 sec exposure time and an altitude above 50*∘*. The light curves were generated with the APEX-2 reduction pipeline (Devyatkin et al. 2010). The pipeline includes aperture photometry, atmospheric extinction correction, differential photometry (using the BESTRED algorithm  (Voss 2006)) and astrometric calibration.

4.1.2 Results

In Fig. 14 and Fig. 15 we present two examples obtained from the TFRM-PSES data. Fig. 14, 2MASS J10144313+5018191, is a newly-discovered variable star, detected both with TFA and TFAW, probably a $\delta$-Scuti star with a 0.1592-day period. As can be seen, after TFAW is applied, the noise contribution decreases significantly and outliers are efficiently removed. The LS power spectrum during the frequency analysis step greatly improves with respect the TFA one presenting higher power of the peaks and an increase in the SDE (10.1 for TFAW in front of 9.8 for TFA). The one in Fig. 15, a W Ursae Majoris-type variable star with a 0.371018-day period (the maximum peak in the LS power spectrum corresponds half the cataloged period) is cataloged as NSVS 4921994 (Woźniak et al. 2004). As in the previous case, the SNR of the signal is improved and the correct period is recovered from the TFAW LS power spectrum. Also, the power spectrum itself presents higher peaks and an increased SDE (10.31 for TFAW in front of 10.1 for TFA).

The small increase in the SDEs of both example light curves for TFAW compared to TFA is to be expected from Figure 3 as they are in the high SNR regime.

4.2 TFAW real light curve quantitative performance

Tests made on both simulated and real light curves show that after TFAW is applied, signals affected by systematics and noise are recovered with no shape and period distortion. Fig. 16 top panel shows the TFRM-PSES photometric precision (standard deviation $\sigma$) vs. R magnitude for both TFA and TFAW light curves. The set of the 6485 displayed light curves correspond to the same field, observed over three years. TFAW photometric performance is not only consistently better than TFA over all the R magnitude range, but also shows a closer fit with the expected blue-lined trend-free stochastic noises performance, especially the faint end. Fig. 16 bottom panel shows, similarly to Fig.5 in Kovács et al. (2005), the decrease in the standard deviations of the same sample of 6485 TFAW-filtered TFRM-PSES light curves compared to the original TFA ones. It can be seen, that for almost all the TFA standard deviation range (conversely brighter to fainter magnitudes), the standard deviations of TFAW light curves is $\sim$40$\%$ better than that for TFA light curves. In the case of low standard deviations (i.e. bright stars), the improvement is smaller because the noise contribution is also smaller. The presence of real variable stars in the sample could also explain a fraction of such TFA and TFAW with similar standard deviations, as their intrinsic variabilities dominate the noise contribution (see Fig. 15 as an example).

The results presented in Sect. 4.1.2 in the LS power spectra demonstrate that TFAW does not introduce false periodicities or eliminate any of the signal peaks during the frequency analysis step. For the examples presented before as well as for the rest of TFRM-PSES light curves analyzed, TFAW outlier removal and frequency analysis step results in an overall improvement of the power spectrum and of the SDEs.

Regarding the number of template stars, tests run both with simulated and real light curves show that using 250 template stars gives the best results. TFAW can also be run with less stars (¡10 stars). However, the noise filtering is less efficient and some trends and systematics could still be present in the filtered data. As with the original TFA, it is not necessary to compute the template matrix for each target star, only in the cases in which the latter accidentally coincides with one of the template ones.

The algorithm has been fully implemented and parallelized in Python. The CPU time needed for a 2048 data points, single light curve TFAW analysis using 250 template stars is typically around 60 seconds and about 100 iterations. A full run of the TFAW on 6485 TFRM-PSES light curves with 2048 points per time series, 250 template stars and 90,000 BLS frequency steps takes 7.5 hours of 24 $\times$ CPUs @ 2.00 GHz. As a comparison, TFA runs the same set of light curves with the same parameters and hardware capabilities in 0.5 hours, i.e., 15$\times$ less than TFAW. The main reason for this performance loss is the extra time in the computation of the SWT and the ISWT at each iteration step. The current PyWavelets implementation we are using for these transforms is coded in Cython. Future effort is planned in order to significantly speed up TFAW.

5 Conclusions

TFAW is a totally generic, Python-based, parallelized algorithm useful for any kind of survey which seeks to improve the performance of signal detection, reconstruction and characterization, provided that a set of comparison light curves sharing the same systematics and trends as the target time series is available. TFAW makes use of an improved frequency analysis step to determine if a significant period is present in the data. This improvement consist on the introduction of a outlier-removal tool as well as the search of periods using a de-noised light curve obtained through the SWT of the original one. Once this period is found, the use of the SWT allows a more precise estimation of the signal’s shape than bin averaging. The algorithm then decouples a large component of the noise contribution from the raw signal, making use of a modified version of the original TFA iterative signal reconstruction approach based on the SWT estimation of the noise at each iteration step, leading to an overall SNR improvement without alteration of the signal’s time sampling or astrophysical characteristics.

Tests conducted on both simulated and real TFAW-filtered light curves show an improvement of $\sim$40$\%$ (although it can be higher) in their standard deviations with respect to the ones detrended with TFA, leading to a better characterization of the signal, without modifying its features (i.e. amplitude, phase, shape or depth). Results obtained for the frequency analysis step show that TFAW does not introduce any false periodicity, and that it can improve the overall power spectra of non-multiperiodic signals and the SDE of the peaks. In the case of simulated transits, TFAW frequency analysis step improves the transit detection rate a factor $\sim$2-5$\times$ for the low SNR signals with respect TFA and increases the SDEs up to a factor $\sim$2.5$\times$. In addition, TFAW is able to detect the transits for signals with $\sim$2$\times$ higher standard deviation. We also show that the SWT signal approximation provides a closer representation of the underlying signal with respect to bin averaging. In the case of planetary transits, this improvement is due to a better fit of the ingress and egress profiles. Also, we demonstrate that the TFAW-filtered light curve yields better MCMC posterior distributions, diminishes the bias in the fitted transit parameters and their uncertainties and narrows the credibility intervals up to a factor $\sim$10$\times$ for simulated transits.

We have also shown that TFAW is able to isolate the different underlying signals within a light curve with multiple periodic signals, such as multi-transit signals, transients, modulations or other kinds of stellar variabilities. As for single-period light curves, the SNR of each signal contribution improves due to the noise filtering capabilities of TFAW.

It is important to note that TFAW needs an even number of data points to work that will depend on the number of levels in which the signal wants to be decomposed (ideally, the length of the time series should be a power of two). Also, as described in Section 2.4, better results are obtained if the number of SWT levels is above ten, ensuring that the noise and signal contributions are as separated as possible in the SWT decomposition of the signal. Compared to the original TFA, TFAW requires fewer reference stars to create the template used for the signal reconstruction and the noise filtering.  Kovács et al. (2005) suggests the use of $\sim$600 template stars, while TFAW operates stably with 250 template stars (in the data presented in this paper).

In this paper we have applied TFAW to thousands of real light curves from TFRM-PSES survey.

We are planning to present a dedicated TFAW test run on large numbers of The Evryscope (Law et al. 2015, 2016) light curves in a forthcoming paper.

TFAW algorithm also has the potential to improve the planet detection and characterization capabilities of space-based data sets, such as from CoRoT and Kepler; this study will be detailed in a forthcoming paper.

Acknowledgements.

DdS acknowledges financial support from RACAB and Universitat de Barcelona. OF acknowledges funding support by the grant MDM-2014-0369 of the ICCUB (Unidad de Excelencia ‘María de Maeztu’). The TFRM project is supported, in part, by the Departament d’Empresa i Coneixement de la Generalitat de Catalunya. We also acknowledge support by the Spanish Ministerio de Economía y Competitividad (MINECO/FEDER, UE) under grants AYA2013-47447-C3-1-P, AYA2016-76012-C3-1-P, MDM-2014-0369 of ICCUB (Unidad de Excelencia ’María de Maeztu’), and the Catalan DEC grant 2014 SGR 86. The authors would like to acknowledge the comments and suggestions of the referee that have helped to improve this work. Also, we would like to acknowledge Ignasi Ribas, Jaime Boloix, Albert Rosich, and Juan Carlos Morales for their contribution to the TFRM-PSES survey, the support team of the TFRM for the commissioning and maintenance of the telescope, as well as to José M. Solanes for contributing with his computing resources.