Frequency (Hz)

An Efficient Forecasting Approach to Reduce Boundary Effects in Real-Time Time-Frequency Analysis

Adrien Meynard, Hau-Tieng Wu

Abstract-Time-frequency (TF) representations of time series are intrinsically subject to the boundary effects. As a result, the structures of signals that are highlighted by the representations are garbled when approaching the boundaries of the TF domain. In this paper, for the purpose of realtime TF information acquisition of nonstationary oscillatory time series, we propose a numerically efficient approach for the reduction of such boundary effects. The solution relies on an extension of the analyzed signal obtained by a forecasting technique. In the case of the study of a class of locally oscillating signals, we provide a theoretical guarantee of the performance of our approach. Following a numerical verification of the algorithmic performance of our approach, we validate it by implementing it on biomedical signals.

Index Terms-Boundary effects, time-frequency, forecasting, nonstationarity

I. INTRODUCTION

I N any digital acquisition system, the study and the interpretation of the measured signals generally require an analysis tool, which enables researchers to point out the useful characteristics of the signal. The need for signal analysis arises from various signals, ranging from audio [1], [2], mechanical [START_REF] Peng | Vibration signal analysis and feature extraction based on reassigned wavelet scalogram[END_REF], or biomedical signals [START_REF] Akay | Detection and Estimation Methods for Biomedical Signals[END_REF]. For instance, biomedical signals, such as photoplethysmogram (PPG), contain several characteristics, including respiratory rate or blood pressure, that cannot be interpreted from its run-sequence plot in the time domain. An analysis tool would make possible the extraction of these useful characteristics.

Usually, the measured signals exhibit nonstationary behavior, and the observed quantities might be interfered with by transient phenomena that can vary rapidly and irregularly. In this paper, we focus on oscillatory time series. These signals might, for example, oscillate fast with large amplitudes at one moment, and then oscillate slowly with small amplitudes at the next moment. In order to adapt the analysis to nonstationarities, local spectral analysis is generally performed [START_REF] Stoica | Spectral analysis of signals[END_REF], [START_REF] Matz | Generalized evolutionary spectral analysis and the Weyl spectrum of nonstationary random processes[END_REF] A. Meynard is the corresponding author (e-mail: adrien.meynard@duke.edu).

short-time Fourier transform [START_REF] Gröchenig | Foundations of time-frequency analysis, ser. Applied and Numerical Harmonic Analysis[END_REF] (STFT), a typical tool built for this purpose, enables the determination of the local frequency content of a nonstationary signal.

Windowing is a common method for performing local analysis. Among many others, STFT [START_REF] Flandrin | Time-frequency/time-scale analysis[END_REF], continuous wavelet transform (CWT) [START_REF] Daubechies | Ten lectures on wavelets, ser. CBMS-NSF Regional Conference Series in Applied Mathematics[END_REF], synchrosqueezing transform (SST) [START_REF] Daubechies | Synchrosqueezed wavelet transforms: An empirical mode decomposition-like tool[END_REF], and reassignment [START_REF] Auger | Time-frequency reassignment and synchrosqueezing: An overview[END_REF] (RS) are representations that fall back on the use of an analysis window. Let x : I → R denote the observed signal, where I denotes the finite interval where the signal is measured. Let g s : R → R denote the analysis window, where s is a shape parameter. The support of g s is localized around the origin and is small with respect to |I|. The translation operator is T τ defined as:

T τ f = f (t -τ) , ∀ f : R → R .
Then, the local analysis of x around the instant τ ∈ I relies on the evaluation of the following inner product: V x (s, τ) = x, T τ g s I .

(

A major shortcoming of this technique occurs when analyzing the signal x near the boundaries of the interval I. Clearly, at these points, half of the information is missing. Consequently, the results of the inner product (1) are distorted. This phenomenon is usually understood as the boundary effect. The result of the SST of a PPG is shown in the lower-left corner of Fig. 1 (see Section IV-D3 for a comprehensive description). The distortion resulting from boundary effects is clearly visible on the right side of this representation-this area is highlighted by the red dashed rectangle. Indeed, while in the major left part of the image, clear lines stand out, they become blurred as they approach the right boundary of the image. Therefore, estimations of signal characteristics, like instantaneous frequencies [START_REF] Delprat | Asymptotic wavelet and Gabor analysis: extraction of instantaneous frequencies[END_REF] or amplitudes, from this TF representation appear to be imprecisely determinable (or even likely to fail) in the vicinity of the boundaries. Moreover, this boundary effect would unavoidably limit applying these window-based TF analysis tools for realtime analysis purposes. It is thus desirable to have a solution to eliminating the boundary effects. Attempts to minimize the boundary effects generally consists in softening the discontinuity on signal edges. Usually, the approaches can be classified into two main classes, choosing a proper analysis window and extending/forecasting the signal. In the first class, judicious choice of analysis window whose support does not interfere with the boundary points can minimize the occurrence of aberrant patterns near the boundaries of the TF plane, for instance, [START_REF] Chui | Numerical Methods in Approximation Theory[END_REF], [START_REF] Depczynski | The fast wavelet transform on compact intervals as a tool in chemometrics: II. Boundary effects, denoising and compression[END_REF]. Due to the specific relationship between the chosen analysis windows and the TF analysis tool, these techniques do not make it possible to reduce the boundary effects of all TF representations. Another natural idea, the second class, consists of carrying out a preliminary step of extending, or forecasting, the signal beyond its boundaries. Due to its flexibility, various forecasting schemes have been proposed. For example, there exist simple extension schemes that do not take into account the dynamical behavior of the signal, such as zero-padding, periodic extension, symmetric extension [START_REF] Kharitonenko | A wavelet transform with point-symmetric extension at tile boundaries[END_REF], [START_REF] Chen | Symmetric extension methods for m-channel linear-phase perfect-reconstruction filter banks[END_REF], or polynomial extrapolation [START_REF] Williams | A discrete wavelet transform without edge effects using wavelet extrapolation[END_REF]. There exist extension schemes based on physically relevant dynamical models, such as the Extended Dynamic Mode Decomposition [START_REF] Williams | A Data-Driven Approximation of the Koopman Operator: Extending Dynamic Mode Decomposition[END_REF] (EDMD) and the Gaussian process regression [START_REF] Rasmussen | Gaussian Processes for Machine Learning, ser. Adaptive Computation and Machine Learning series[END_REF], [START_REF] Roberts | Gaussian processes for time-series modelling[END_REF] (GPR). In speech processing, dynamic mode predictors have also been proposed [START_REF] Vargas | Speech analysis and synthesis based on dynamic modes[END_REF] to forecast signals falling into the socalled source-filter model. This gradient-based technique relies on the shadowing approach proposed in [START_REF] Grebogi | Shadowing of physical trajectories in chaotic dynamics: Containment and refinement[END_REF]. There also exist extension schemes based on stochastic models, such as the Trigonometric, Box-Cox transformation, ARMA errors, Trend and Seasonal components (TBATS) algorithm [START_REF] Livera | Forecasting time series with complex seasonal patterns using exponential smoothing[END_REF], the dynamic linear models [START_REF] West | Bayesian forecasting and dynamic models[END_REF]. In the physically relevant dynamical models, the oscillation or trend, are usually modeled as the mean of a stochastic process, while in the stochastic models, the oscillation or trend are modeled by the covariance structure of the stochastic process. It is also possible to consider polynomial regression [START_REF] Fan | Local polynomial modelling and its applications: monographs on statistics and applied probability 66[END_REF] or modeling the mean by splines [START_REF] Hall | Theory for penalised spline regression[END_REF], kernel functions [START_REF] Chang | Training and testing low-degree polynomial data mappings via linear svm[END_REF], and wavelets [START_REF] Marron | Exact risk analysis of wavelet regression[END_REF],

or nonparametric regression [START_REF] Fan | Local polynomial modelling and its applications: monographs on statistics and applied probability 66[END_REF], to estimate the mean function before forecasting the signal. Neural network models, such as the long short-term memory model [START_REF] Vlachas | Data-driven forecasting of high-dimensional chaotic systems with long short-term memory networks[END_REF], is another approach. The above list is far from exhaustive, and we refer readers with interest to [START_REF] Hyndman | Forecasting: principles and practice[END_REF] for a friendly monograph on the general forecasting topic. While the model-based extension scheme gives betterextended signals than those simple extension without dynamical models, they generally have a great computational cost.

In this paper, we propose a fast extension algorithm based on a simple dynamical model that optimizes the trade-off between the extension quality and the computational cost, so that the real-time analysis can be achieved. The algorithm is composed of two steps. Assuming that the support of the analysis window is smaller than 2∆, the local analysis near the boundary of I is now possible without lack of information thanks to knowledge brought by the extension.

Thus, assuming that the quality of the extension step is sufficient, the analysis results obtained that way will be less sensitive to the boundary effects than the result of the analysis tool applied directly to the nonextended signal. We claim and prove that forecasting oscillatory signals based on the simple dynamical model combined with the simple least square approach is sufficient for reducing the boundary effects for the TF analysis, or other kernel-based analysis. To our knowledge, such theoretical analysis does not exist in the literature. See the bottom of Fig. 1-in particular, the regions indicated by the dashed boxes-for a snapshot of the result. The main benefit of this simple approach is a numerically efficient solution with a theoretical guarantee for realtime analysis purposes.

The paper is organized in the following way. In section II, we provide an extension method based on a linear dynamic model. We derive the corresponding algorithm for boundary effects reduction. In Section III, we show that the dynamic model we consider is sufficient to extend signals taking the form of sums of sine waves. An evaluation of the theoretical performance of our algorithm on this class of signals is given in Section III. In Section IV, we compare our extension method with more sophisticated methods such as EDMD, GPR, or TBATS. We show that our algorithm gives fast results of reasonable quality. Finally, we evaluate the performance of our boundary effects reduction algorithm on biomedical signals, such as respiratory signals, and compare it to the theoretical results.

II. ALGORITHM

As explained above, the algorithm for the reduction of boundary effects on TF representations relies on extending the signal by forecasting it before applying the TF analysis.

We start with the notation. Let x : R → R denote a continuous-time signal. In this work, we consider a finite-length discretization of that one. Thus, the sampled signal x, whose length is denoted by N, is such that

x[n] = x n f s , ∀n ∈ {0, . . . , N -1} ,
where f s denotes the sampling frequency. Let M and K be two positive integers such that M < K and K + M ≤ N. Then, for all k ∈ {0, . . . , K -1}, we extract from x ∈ R N the subsignal x k ∈ R M given by:

x k =    x[N -K -M + k] . . . x[N -K + k -1]    . ( 2
)
These subsignals are gathered into the matrix X ∈ R M×K such that:

X = x 0 • • • x K-1 .
Notice that these subsignals are overlapping each other. Indeed, x k+1 is a shifting of x k from one sample. We also consider the matrix Y ∈ R M×K given by:

Y = x 1 • • • x K .
The boundary effect reduction algorithm is based on manipulating X and Y.

The pseudo-code of the proposed real-time algorithm to reduce boundary effects on windowing-based TF representations is shown in Algorithm 1. We coined the algorithm BoundEffRed. Below, we detail the algorithm, particularly the signal extension SigExt in Algorithm 2.

Algorithm 1

Tackling boundary effects of a TF representation in real-time.

F x = BoundEffRed(x, M, K, L, F ) Inputs: x N , M, K, N 0 , L, F
while N increases do Real-time input: x N Forecasting step.

• Signal extension: xN+L = SigExt(x N ). Representation estimation step.

• Extended representation evaluation: F ( xN+L ).

• Restriction of F ( xN+L ) to the current time interval (see [START_REF] Daubechies | Ten lectures on wavelets, ser. CBMS-NSF Regional Conference Series in Applied Mathematics[END_REF]) to obtain F

(N) x = F ext (x N ).
Y = f (X) . (3)
In a general framework, forecasting means estimating the function f from the observed values taken by the signal, in order to predict its future values. For instance, the dynamic mode decomposition [START_REF] Schmid | Dynamic mode decomposition of numerical and experimental data[END_REF], [START_REF] Williams | A Data-Driven Approximation of the Koopman Operator: Extending Dynamic Mode Decomposition[END_REF] or other more complicated models [START_REF] Roberts | Gaussian processes for time-series modelling[END_REF], [START_REF] Livera | Forecasting time series with complex seasonal patterns using exponential smoothing[END_REF], [START_REF] West | Bayesian forecasting and dynamic models[END_REF], [START_REF] Vlachas | Data-driven forecasting of high-dimensional chaotic systems with long short-term memory networks[END_REF], allow this by setting additional constraints on the behavior of f . We will see, in Section III, that considering such a complex dynamic model is not necessary for the study of the oscillatory signals of interest to us. That is why we consider here a naive dynamical model, assuming that we have the following relation:

Y = AX , (4)
where A ∈ R M×M . In other words, we adopt a classical strategy in the study of dynamical systems, the linearization of a nonlinear system, when the system is sufficiently regular. Notice that this linearized dynamical model can be written equivalently according to the subsignals x k , as:

x k+1 = Ax k , , ∀k ∈ {0, . . . , K -1} . (5)
b) Forecasting: Determining the forecasting amounts to estimating the unknown matrix A. Indeed, let à denotes the estimate of A. We then obtain the forecasting of x k+ by recursively applying the linear relation (4):

xK+ = Ã Ã • • • Ã times x K = Ã x K .
Let e M denote the unit vector of length M given by e M = 0 • • • 0 1

T . Consequently, given definition (2), the forecasting of the signal at time N-1+ f s is provided by

x[N -1 + ] = e T M xK+ = α ( ) x K , (6)
where α ( ) denotes the last row of à , i.e., α ( ) = e T M à . c) Model Estimation: To estimate the matrix A, we consider the simple but numerically efficient least square estimator. That is, we solve the following problem:

à = arg min α L(A) , (7)
where the loss function L is given by:

L(A) = Y -AX 2 = K-1 ∑ k=0 x k+1 -Ax k 2 .
Therefore, solving the problem [START_REF] Gröchenig | Foundations of time-frequency analysis, ser. Applied and Numerical Harmonic Analysis[END_REF], i.e., ∇L( Ã) = 0, gives the following estimate à of the dynamical model matrix A:

à = YX T (XX T ) -1 . ( 8
)
Remark 1. This expression clearly shows that the matrix à takes the following form:

à =         0 1 0 • • • 0 . . . . . . . . . . . . . . . . . . . . . . . . 0 0 • • • • • • 0 1 α 1 • • • • • • • • • α M         .
Then, except for the row vector α = (α 1 • • • α M ), the matrix A is fully determined by the dynamical model. d) Signal Extension: Since we are building a realtime algorithm, we consider that only the right "side" of the signal has to be extended, the left "side" being fixed since it only concerns the past values of the signal. We therefore construct the extended signal x ∈ R N+L concatenating the observed signal x, and the forward forecasting xfw ∈ R L . We summarize the extension step in Algorithm 2.

Algorithm 2 Signal extension. x = SigExt(x, M, K, L)

Inputs: x N , M, K, L Forward forecasting.

• Estimation of the matrix à via equation [START_REF] Flandrin | Time-frequency/time-scale analysis[END_REF].

• Forecasting xfw ∈ R L obtained applying equation [START_REF] Matz | Generalized evolutionary spectral analysis and the Weyl spectrum of nonstationary random processes[END_REF] with ∈ {1, . . . , L}.

Output:

Extended signal xN+L = x N xfw T .

B. Step 2: Extended Time-Frequency Representation

Let F : R N → C F×N generically denotes the TF representation of interest to us, which could be, for instance, STFT, CWT, SST, or RS. Here, F typically denotes the size of the discretization along the frequency axis. Due to the boundary effects, the representation F (x N ) shows undesired patterns near its edges. To alleviate the boundary effects, we apply the representation to the estimated extended signal x. This strategy moves the boundary effects out of the time interval I = [0, N-1 f s ]. Finally, the boundary-effects insensitive representation

F ext : R N → C F×N of x N is given for all ν ∈ {0, • • • , F - 1}, n ∈ {0, • • • , N -1} by: F ext (x N )[ν, n] = F ( xN+L )[ν, n] . ( 9
)
This amounts to restricting the representation F ( xN+L ) to the current measurement interval of x N . For the sake of simplicity, we denote the restriction operator by R, where R : C F×(N+L) → C F×N . Consequently, we have:

F ext (x N ) = R (F ( xN+L )) . ( 10
)
We call F ext the boundary-free TF representation.

F (N) x
∈ C F×N is the estimation of the boundary-free TF representation at iteration N in Algorithm 1. For numerical purposes, and to make the real-time implementation achievable, we do not perform a full re-estimation of F ext at iteration N + 1. Instead, the additional knowledge provided by x[N + 1] only influences the values of the last L win columns of F ext (x N+1 ), where L win denotes the window half-length used by the TF representation. Thus, F (N+1) x is obtained by the concatenation of the first N -L win + 1 columns of F (N) x with the last L win columns of F ext (x N+1 ).

III. THEORETICAL PERFORMANCE

A. Signal Model

We model the meaningful part of the observed signal as a deterministic multicomponent harmonic signal; that is, a sum of sine waves:

z[n] = J ∑ j=1 Ω j cos 2π f j n f s + ϕ j , (11)
where J denotes the number of components, Ω j > 0 the amplitude of the jth component, f j its frequency, and ϕ j ∈ [0, 2π) its initial phase. For the sake of simplicity, we make an additional assumption on the frequencies of each component. We assume that for all j ∈ {1, . . . , J}:

∃ p j , p j ∈ N * : f j = p j M f s = p j K f s . (12)
In addition, the observed signal is assumed to be corrupted by an additive Gaussian white noise. Therefore, the measured discrete signal x is written as:

x = z + σw , ( 13
)
where z follows model [START_REF] Auger | Time-frequency reassignment and synchrosqueezing: An overview[END_REF], w is a Gaussian white noise, whose variance is normalized to one, and σ > 0. Clearly, σ 2 denotes the variance of the additive noise σw.

B. Forecasting Error

On the forecasting interval, we decompose the estimated signal x as follows:

x[n] = z[n] + [n] , (14)
where is the forecasting error. When n ∈ I = {0, . . . , N -1}, this error contains only the measurement noise, that is [n] = σw [n]. Outside the interval I, the importance of the forecasting error is also affected by the loss of information resulting from the linearization of the dynamical model we consider in (4). To evaluate the actual behavior of the forward forecasting error [n] when n ≥ N, we determine its first two moments.

1) The mean, or estimation bias, is such that:

µ[n] ∆ = E{ [n]} = E{ x[n]} -z[n] . ( 15
)
Given the forecasting strategy, we have µ[n] = 0 when n ∈ I and

µ[n] = E{α (n-N+1) }z K + σE{α (n-N+1) w K } -z[n] (16) when n ≥ N.
2) The covariance is given by:

γ[n, n ] ∆ = E{( [n] -µ[n]) [n ] -µ[n ] } = E{ x[n] x[n ]} -z[n]z[n ] -µ[n]z[n ] -µ[n ]z[n] -µ[n]µ[n ] .
Thus by definition of the noise, we have γ[n, n ] = σ 2 δ n,n when (n, n ) ∈ I 2 . When n ≥ N, let us denote = n -N + 1. Then, we have two cases. (i) If n ∈ I:

γ[n, n ] = σE{w[n ]α ( ) }z K + σ 2 E{w[n ]α ( ) w K } . ( 17
) (ii) If n = N -1 + λ ≥ N: γ[n, n ] = z T K E α ( ) T α (λ) z K + σE{α ( ) w K α (λ) }z K + σE{α (λ) w K α ( ) }z K + σ 2 E{α ( ) w K α (λ) w K } -z[n]z[n ] -z[n]µ[n ] -z[n]µ[n ] -µ[n]µ[n ] . (18)
Besides, we recall that γ[n, n ] = γ[n , n]. In Theorem 1, we specify the asymptotic behavior of the forecasting bias and covariance when the dataset size K is great. Theorem 1. Let x ∈ R N be a discrete-time random signal following model [START_REF] Chui | Numerical Methods in Approximation Theory[END_REF]. Let x denotes its forecasting, obtained using the extension Algorithm 2. Let n ≥ N be a sample index. Then, the first-order moment of the forecasting error

[n] in (14) defined in [START_REF] Kharitonenko | A wavelet transform with point-symmetric extension at tile boundaries[END_REF] satisfies

|µ[n]| ≤ a (n) 0 σ 2 + 1 K a (n) 1 σ 2 + a (n) 2 + o 1 K (19)
as K → ∞, where a

(n) 0 , a (n)
1 , and a

(n)

2 are positive quantities, independent of K and σ. Its second-order moment γ[n, n ] satisfies the following approximation equations:

(i) if n ∈ I = {0, . . . , N -1}:

γ[n, n ] ≤ b (n,n ) 0 σ 2 + 1 K b (n,n ) 1 + b (n,n ) 2 σ 2 + o 1 K (20) as K → ∞, where b (n,n ) 0 , b (n,n ) 1
, and b

(n,n ) 2 are positive quantities, independent of K and σ; (ii) if n ≥ N: γ[n, n ] ≤ c (n,n ) 0 σ 2 + 1 K c (n,n ) 1 σ 2 + c (n,n ) 2 + c (n,n ) 3 σ 2 + o 1 K (21)
as K → ∞, where c

(n,n ) 0 , c (n,n ) 1 , c (n,n ) 2
and c

(n,n ) 3 are positive quantities, independent of K and σ.

Proof. See the Supplementary Materials. The proof is mainly based on Taylor expansions of the forecasting error. Nonoptimal values of a

(n) 0 , a (n) 1 ,. . . , b (n,n ) 0 ,. . . , c (n,n ) 3
are also provided in the proof.

The forecasting bias and covariance asymptotically depend linearly on the ratio 1 K . This shows the need to use a sufficiently large dataset to obtain an accurate forecast. Ideally when K → ∞, the forecasting error would behave like the measurement noise σw, i.e., a zero-mean white noise whose variance is of the order of σ 2 . However, Theorem 1 shows that the estimator may remain asymptotically biased (the first term in the bound of equation ( 19) being independent of K). Concerning the covariance of the forecasting error, the expected asymptotic behavior is verified. Indeed, when K → ∞, the variance of the forecasting estimator increases at most linearly with the noise variance σ 2 , since is bounded by b

(n,n ) 0 σ 2 .
Assume now that K is sufficiently great and fixed. The noise influences the quality of the forecasting estimator in two ways. On the one hand, when the noise variance σ 2 is great, the bias and the variance of the estimator potentially increase linearly with σ 2 . This reflects the classical behavior of an estimator of a signal degraded by additive noise. On the other hand, due to terms in 1/σ 2 , when the noise variance σ 2 is small, the bias and the variance of the forecasting estimator are no longer controlled via equations ( 19) and ( 21). This illustrates the necessary presence of noise to obtain an accurate prediction. Indeed, if σ is too small, the matrix XX T is ill-conditioned and its inversion in ( 8) is sensitive. The forecasting is then of poor quality. Noise therefore appears as a regularization term to ensure the invertibility of XX T . This seemingly counterintuitive fact indicates the challenge we encounter to do forecasting with the dynamic model (3), the degeneracy. Note that from the lag map perspective, X and Y consist of points located on a low dimensional manifold. Thus, locally the ranks of X and Y are low, which leads to the degeneracy. However, when noise exists, it naturally fills in deficient ranks, and stabilizes the algorithm.

The dependencies of the forecasting quality on the subsignals lengths M and the forecasting index = n -N + 1 are hidden in the expression of the parameters a

(n) 0 , a (n) 1 ,. . . , b (n,n ) 0 ,. . . , c (n,n ) 3
. The numerical results, discussed in Section IV-C1, allow us to better evaluate these dependencies.

C. Adaptive Harmonic Model

One can extend the previous result to the case where the instantaneous frequencies and amplitudes of the components of the deterministic part of the observed signal are slowly varying. It is an AM-FM model that, in its is continuous-time version, takes the following form

z(t) = J ∑ j=1 a j (t) cos(2πφ j (t)) , ( 22
)
where a j (t) and φ j (t) describe how large and fast the signal oscillates at time t. Clearly, ( 11) is a special case satisfying the AM-FM model when the a j and φ j are both constants. In many practical signals, the amplitude and frequency do not change fast. It is thus reasonable to further restrict the regularity and variations of the instantaneous amplitudes and frequencies of these AM-FM functions. When the speeds of variation of the instantaneous amplitudes a j and frequencies φ j are small, the signal can be "locally" well approximated by a harmonic function in [START_REF] Auger | Time-frequency reassignment and synchrosqueezing: An overview[END_REF]; that is, locally at time t 0 , z(t) can be well approximated by z 0 (t) := ∑ J j=1 a j (t 0 ) cos(2π(φ j (t 0 )t 0 φ j (t 0 ) + φ j (t 0 )t)) by the Taylor expansion. An AM-FM function satisfying the slow variation of the instantaneous amplitudes a j and frequencies φ j is referred to the adaptive harmonic model (AHM) (see [START_REF] Chen | Non-parametric and adaptive modelling of dynamic periodicity and trend with heteroscedastic and dependent errors[END_REF], [START_REF] Daubechies | ConceFT: concentration of frequency and time via a multitapered synchrosqueezed transform[END_REF] for mathematical details). It is thus clear that the forecasting error caused by the SigExt algorithm additionally depends on the speed of variation of the instantaneous amplitudes a j and frequencies φ j . For the forecasting purpose, it is thus clear that when the speed of variation of a j and φ j are slow K can be chosen large while the signal can be well approximated by [START_REF] Auger | Time-frequency reassignment and synchrosqueezing: An overview[END_REF]. Hence, Theorem 1 can still be applied to approximate the error. However, how to determine the optimal K based on the signal depends on the application and is out of the scope of this paper. It will be explored in future work for specific application problems.

The models described in Sections III-A and III-C consider the meaningful part of the signal to be a deterministic component. These models are purposely adapted to signals showing local line spectra. The physiological signals we are interested in, such as respiratory or cardiac signals, have this characteristic (see Section IV-D). Signals with wider local spectra, such as electroencephalography signals, do not fall into this category, and are more faithfully modeled as random signals. Thus, the theoretical justifications proposed above are no longer applicable to guarantee the forecasting quality of SigExt.

D. Performance of the Boundary Effects Reduction

While we do not provide a generic proof of the boundary effect reduction on any TF analysis tools, we discuss a particular case of SST. Since SST is designed to analyze signals satisfying the AHM, as is discussed in Section III-C, Theorem 1 ensures that the forecasting error defined in ( 14) is controlled and bounded in terms of mean and covariance. Recall Theorem 3 in [START_REF] Chen | Non-parametric and adaptive modelling of dynamic periodicity and trend with heteroscedastic and dependent errors[END_REF], which states that when the additive noise is stationary and may be modulated by a smooth function with a small fluctuation, the SST of the observed signal remains close to the SST of the clean signal, throughout the TF plane. We refer the reader to [START_REF] Chen | Non-parametric and adaptive modelling of dynamic periodicity and trend with heteroscedastic and dependent errors[END_REF] for a precise quantification of the error made in the TF plane, which depends notably on the covariance of the additive noise and on the speed of variation of the amplitudes and instantaneous frequencies composing the signal. In our case, while the noise of the historical data is stationary, the forecasting error depends on the historical noise, and hence the overall "noise" is nonstationary. However, the dependence only appears in the extended part of the signal. We claim that the proof of Theorem 3 in [START_REF] Chen | Non-parametric and adaptive modelling of dynamic periodicity and trend with heteroscedastic and dependent errors[END_REF] can be generalized to explain the robustness of SST to the noise, and a systematic proof will be reported in our future theoretical work. We verify experimentally that Algorithm 1 is efficient for a large number of representations in the following Section IV. Therefore, in our case, this means that the boundary effect is strongly reduced since the impact of the forecasting error does SST is limited. An immediate application is that the instantaneous frequencies can now be well estimated continuously up to the edges of the TF plane, and real-time acquisition of the instantaneous frequency and amplitude information is possible.

IV. NUMERICAL RESULTS

For the reproducibility purpose, the MATLAB code and datasets used to produce the numerical results in this section are available online at https://github.com/ AdMeynard/BoundaryEffectsReduction.

A. Extension Methods for Comparison

To highlight the fact that the linear dynamical model is sufficient to catch most of the dynamical behavior of signals following the AHM, we compare the performance of Algorithm 1 with a simple extension obtained by pointwise symmetrization [START_REF] Kharitonenko | A wavelet transform with point-symmetric extension at tile boundaries[END_REF]. We also evaluate the performance of reference forecasting algorithms that could be used for extending such signals. These methods are:

• The EDMD has been developed by Williams et al. [START_REF] Williams | A Data-Driven Approximation of the Koopman Operator: Extending Dynamic Mode Decomposition[END_REF]. The proposed algorithm is a way to obtain an approximation of the so-called Koopman operator of the observed system, which theoretically allows catching dynamic of nonlinear systems [START_REF] Korda | Linear predictors for nonlinear dynamical systems: Koopman operator meets model predictive control[END_REF].

• The GPR [START_REF] Rasmussen | Gaussian Processes for Machine Learning, ser. Adaptive Computation and Machine Learning series[END_REF] is a method relying on a probabilistic dynamical model. That one is based on the Gaussian process structure, and therefore offers more flexibility in the type of dynamic that could be modeled than the linear model (4).

• The TBATS method [START_REF] Livera | Forecasting time series with complex seasonal patterns using exponential smoothing[END_REF] is based on a classical decomposition of times series into a trend, seasonal and ARMA components, with a specific dynamic for the seasonal component. This model demands the estimation of numerous parameters and, by implication, may require a long computation time.

B. Evaluation Metrics

The first evaluation metric is quantifying the forecasting quality (i.e. not depending on ) in the time domain. We evaluate the Experimental Mean Square Error MSE xp ( x) of the forward forecast extended signals, namely,

MSE xp ( x) = 1 L x -x ext 2 (23) = 1 L L ∑ =1 µ xp [N-1+ ] 2 +γ xp [N-1+ , N-1+ ].
where x ext is the ground-truth extended signal; that is,

x ext := x[-L] • • • x[N -1 + L] .
Then, as long as the bias µ[N -1 + ] and the variance γ[N -1 + , N -1 + ] of the forecasting estimator remain small for all , the MSE takes small values either.

The second evaluation metric is evaluating the quality of the boundary effect reduction directly on the TF representation. We compare the obtained TF representation to the optimal TF representation F opt N (x), defined as the restriction of the representation of the ground-truth extended signal x ext ; that is, we have

F opt (x N ) = R F (x ext ) , ( 24
)
where R is defined in [START_REF] Daubechies | Synchrosqueezed wavelet transforms: An empirical mode decomposition-like tool[END_REF]. To compare different techniques, we use a criterion proposed in [START_REF] Daubechies | ConceFT: concentration of frequency and time via a multitapered synchrosqueezed transform[END_REF] that quantifies the optimal transport (OT) distance between a given TF representation and the optimal one. In general, the OT distance is a quantity that measures how different two probability density functions are. Denote a TF representation as Q. For t fixed, we consider the probability density function defined by

p t Q (ξ) = |Q(ξ, t)| 2 R |Q(ν, t)| 2 dν
; that is, normalizing the spectral information at time t. At each time t, the OT distance d t between two TF representations Q and F opt at time t is defined by

d t (Q, F opt ) = R Pt Q (ξ) -P t F opt (ξ) dξ , where P t Q (ξ) = ξ -∞ p t Q (ν)dν and Pt F opt (ξ) = ξ -∞ pt F opt (ν)dν.
In other words, the OT distance between Q and F opt at time t is given by the L 1 norm of the associated cumulative density function. Note that d t (Q, F opt ) quantifies the proximity between the estimated and actual TF representations at time t. With the OT distance, we define the performance index D(Q, F ) for the reduction of boundary effects of a given TF representation Q over another TF representation F by

D(Q, F ) = I d t Q, F opt dt I d t F , F opt dt . ( 25
)
In other words, D(Q, F ) is the ratio between its averaged OTD to the optimal TF representation F opt (see [START_REF] West | Bayesian forecasting and dynamic models[END_REF]) and the averaged OTD of the original TF representation F to F opt . Thus, D(Q) < 1 means a reduction of the boundary effects. Let us evaluate the quality of the boundary effects reduction on biomedical signals.

C. Simulated Signals

We first evaluate the quality of the forecasting step and compare it to the theoretical results provided by Theorem 1. The level of the forecasting error depends on at least two parameters:

• The noise variance σ 2 .

• The size of the training dataset K. In Sections IV-C1 and IV-C2, we study the influence of these parameters. A comparison with the theoretical results of Section III is also available. 1) Sum of sine waves: We proved that the linear dynamic model is sufficient to catch the dynamical behavior of signals taking the form [START_REF] Auger | Time-frequency reassignment and synchrosqueezing: An overview[END_REF]. In order to validate this theoretical result, we apply the forecasting Algorithm 2 to a large number of realizations of the random vector x of size N = 10 4 , following model ( 13), and such that the deterministic component z takes the form:

z[n] = cos 2π p 1 n M +A cos 2π p 2 n M , ∀n ∈ {1, . . . , N},
with M = 150, p 1 = 10, p 2 = 33 and A = 1.4. Besides, the additive noise is chosen to be Gaussian: w ∼ N (0, I). a) Influence of the Noise Variance σ 2 : Here, the size of the training dataset is set to K = 450. Then, the forecasting algorithm is run on 500 realizations on the discrete signal x for 200 different values of σ, logarithmically equispaced from 10 -7 to 10 -1 . For each of these values, we determine the experimental bias, denoted as µ xp [N -1 + ], and experimental variance, denoted as 2 shows the experimental forecasting variance as a function of the noise variance for three different values the forecasting sample index: = 1, = 10, and = 100. First, when the noise is weak, the variance of the forecasting estimator is dramatically high. For instance, in this specific example, the estimator x[N -1 + ] for = 10 is inaccurate as soon as the noise variance σ 2 is lower than 4 × 10 -14 (see the red line in Fig. 2). This validates the result of Theorem 1, highlighting the necessary presence of noise to ensure the numerical stability of the estimator. Furthermore, as soon as σ 2 is higher than this critical threshold, the estimator variance increases linearly with σ 2 , as predicted by the bound [START_REF] Vargas | Speech analysis and synthesis based on dynamic modes[END_REF] provided in Theorem 1. Note that the parameter has little influence on the quality of the estimator. This is due to the stationarity of the studied signal, allowing longrange forecasts without deterioration of the estimator.

γ xp [N -1 + , N -1 + ]. Fig.
b) Influence of the Training Dataset Size K: Here, the noise variance σ is set to σ = 10 -2 . Then, SigExt is run on 500 realizations of the discrete signal x for 200 dif- As anticipated by Theorem 1, the variance of the estimator decreases linearly with 1 K as long as the dataset size is large enough. Indeed, in the log-log plot displayed in Fig. 3, as soon as K > 700, each of the three curves representing the forecasting variance as a function of K tends towards an asymptote of slope -1. In practice, the asymptotic behavior of SigExt, illustrated here on a synthetic signal, can be generalized to the study of reallife signals, as long as the duration of the recorded signal largely exceeds that of the desired extension, which allows the user to set K L. This study overlooks the analysis of the influence of the parameter M, whose influence on the value of the experimental variance is numerically not significant as long as M

2K and < M. The choice of this parameter is especially crucial when the deterministic component of the signal is no longer stationary. The AHM, discussed below, is an example.

2) Adaptive Harmonic Model: Consider a signal satisfying the AHM so that the instantaneous frequencies and amplitudes of its components vary over time. The deterministic component z of the random vector x, constructed following the model [START_REF] Chui | Numerical Methods in Approximation Theory[END_REF], takes the following form, for all n ∈ {1, . . . , N}:

x[n] = cos (2πφ 1 [n]) + R[n] cos (2πφ 2 [n]) ,
where the instantaneous amplitude R is given by:

R[n] = 1.4 + 0.2 cos 4π n N ,
and the instantaneous phases are such that: Besides, the noise is chosen to be Gaussian, w ∼ N (0, I), and we take N = 10 4 , P = 750, p 1 = 10, p 2 = 23. a) Influence of the subsignals length M: To evaluate the forecasting quality, SigExt is applied to 1000 realizations of the above-described signal. We forecast signal extensions of 0.1 s, i.e., L = 700. Table I shows the resulting mean and standard deviation (SD) of the experimental MSE, according to three different values of the subsignals length M. Note that choosing too small an M does not provide enough information to forecast the signal satisfactorily. Furthermore, too large an M makes the algorithm insensitive to local nonstationarities produced by frequency variations. The user must then make a compromise on the choice of M in order to optimize the performance of the algorithm. A choice of M such that the subsignals contain about a dozen oscillations is generally appropriate.

φ 1 [n] = p 1 P n + 0.01 2π cos 2π n N φ 2 [n] = p 2 n P + 20 2N f s n 2 .
b) Comparison of the extensions methods: The same set of experiments is run on the extension methods outlined in Section IV-A. The results, displayed in Table I, show that the naive extension we propose gives satisfying results. Indeed, when the subsignal length M is optimally chosen, SigExt performance is of the same order as that of ore sophisticated methods, like GPR or TBATS. Besides, even though these methods could be slightly more robust to perturbations, they are substantially limited by the computing time they require, which prevents them from being used to exploit realtime data. Thus, SigExt is the extension method that optimizes the trade-off between forecasting quality and computing time. Besides, the last column of Table I gives the performance index D of the boundary-free STFT associated with each extension method. This illustrates the ability of our algorithm to reduce boundary effects to a level comparable with what is possible with the GPR or TBATS extensions. The performance of BoundEffRed on other TF representations is more thoroughly analyzed on real-life signals in the following.

D. Real Physiological Signals 1) Respiratory Signal:

We first consider a thoracic movement signal recorded by the piezo-electrical sensor of length 6 hours and 20 minutes. The signal is sampled at f s = 100 Hz. A small portion of this signal is displayed in Fig. 4. From that long signal, we build a dataset of 379 nonoverlapping signals of 60 seconds, i.e., N = 6000. For each segment, we apply the forecasting methods introduced in Section IV-C2, including the SigExt method detailed in Algorithm 2. We forecast 7-second-long extensions on each segment of the signal, corresponding to L = 700. Thus, in order to catch the dynamical behavior, the size of the training signal M is chosen so that M = 1.5L . As a result of Section IV-C1, we take: K = 2.5M . The extensions obtained on one of these subsignals are shown in Fig. 4. The resulting MSEs of different methods are given in Table II. Note that the results are given in "mean ± standard deviation" format. The TBATS extension method is not implemented, as its excessive computing time makes it impossible to implement it on 6000 segments within a reasonable period of time. Note that the MSE of SigExt is on average higher than the MSE of EDMD or the symmetric extension, with a huge standard deviation. This variability is caused by the presence of few segments contaminated by artifacts so that it is unpredictable via a too simple dynamical model like (4). The left of Fig. 5 illustrates one of these outliers, where SigExt fails to catch the fast varying dynamic of the instantaneous amplitude to satisfactorily forecast the signal. The EDMD and symmetric extensions are more robust to those situations, as shown in this example, and in Table II. Nevertheless, SigExt provides a sufficiently relevant extension to give TF representations sparingly affected by boundary effects. On the right of Fig. 5, we display a comparison between the right boundary of SST of the same segment of signal (top-right), and its boundary-free SST obtained after the SigExt forecasting (bottom-right). The extension of the instantaneous frequency visible on the right side of the image illustrates the reduction of boundary effects produced despite an inaccurate signal forecasting.

We then apply BoundEffRed for diverse TF representations: STFT, SST, RS, as well as Concentration of Frequency and Time (ConceFT), a generalized multitaper SST-based representation introduced in [START_REF] Daubechies | ConceFT: concentration of frequency and time via a multitapered synchrosqueezed transform[END_REF]. We mention that unlike other representations, the standard RS is not causal and rigorously requires knowing the STFT on the whole time-frequency plane before making any reassignment. To enable real-time implementation of RS, we use here a real-time version of RS proposed in [START_REF] Lin | ConceFT for Time-Varying Heart Rate Variability Analysis as a Measure of Noxious Stimulation During General Anesthesia[END_REF]. We set the extension length L accordingly to the window length used by the TF analysis tool. For instance, the window length used to evaluate the STFT is of 1400 samples. To prevent the STFT from being sensitive to the boundaries, we set L = 700. In this way, the evaluation of the spectral content of the signal near its boundary is not limited by a lack of information in the interval delimited by the window support. From now on, all results are given for L equal to the half of the width of the window used in the TF transform. In Table II, we give the averaged performance index (25) by evaluating the whole TF representations (including boundaries). Even though SigExt performs moderately well in the sense of MSE, the boundary effects are dramatically reduced on the TF representations, and the averaged performance is in the same order of that given by EDMD or GPR. A t-test is performed to compare the performance index of SigExt with those of other methods. Under the null hypothesis that the means are equal, with a 5% significance level, the t-tests show no statistically significant difference between SigExt and EDMD or GPR, regardless of the representation considered.

2) Two-Component Respiratory Signal: Then, we consider a second respiratory signal that contains not only the respiratory cycle, but also a cardiac component known as the cardiogenic artifact [START_REF] Smith | Recognition of cardiogenic artifact in pediatric capnograms[END_REF]. The top of Fig. 6 shows an excerpt of the measured thoracic impedancebased respiratory signal, recorded during a bronchoscopic sedation procedure, along with the SigExt extension. Below, the simultaneously recorded ECG signal is depicted. The magenta dotted lines emphasize the coincidences between the cardiogenic artifact of the respiratory signal and the QRS complexes of the ECG signal. Such cardiac information can be harvested and utilized when other first-line heart rate information resource is not available or broken [START_REF] Lu | Recycling cardiogenic artifacts in impedance pneumography[END_REF].

The top of Fig. 7 displays the ordinary SST and the boundary-free SST of this respiratory signal, sampled at 64 Hz. This TF analysis brings out both components-the respiratory component, indicated by the green arrows, is located around the fundamental frequency 0.3 Hz and its multiples, while the cardiac component, indicated by the blue arrows, is located around 1.5 Hz. Note that, here, boundary effects are only reduced on the right-side boundary of the TF domain, delimited by the red dashed boxes. The middle of Fig. 7 show a zoom near the rightboundary of the ordinary SST, while the bottom of Fig. 7 respectively show the respective zoom on the boundaryfree SST. In this area, the boundary-free SST makes it possible to disentangle the components contained in the signal. Moreover, the performance index of this representation takes the value D = 0.705. This means that BoundEffRed has reduced the right-side boundary effects by about 30% with respect to the ordinary SST. This shows the ability of our algorithm to work on signals containing several nonstationary components.

3) Photoplethysmogram: We consider a 640-second photoplethysmogram (PPG) signal extracted from the Physionet dataset [START_REF] Pimentel | Toward a robust estimation of respiratory rate from pulse oximeters[END_REF], [START_REF] Goldberger | Physiobank, physiotoolkit, and physionet: Components of a new research resource for complex physiologic signals[END_REF], sampled at f s = 125 Hz. A portion of this signal is displayed in Fig. 8. The estimated 5-second extensions of this segment obtained by SigExt, EDMD, and GPR forecastings are superimposed to the ground-truth extension in Fig. 8.

We divide the signal into 32-second-long segments and apply BoundEffRed to each of them. Table III shows the performance index D of the boundary-free TF representations, averaged over the signals. For all the TF representations considered, the results show that SigExt reduces boundary effects about as efficiently as the extensions given by EDMD or GPR. As for the respiratory signal, t-tests show no statistically significant difference between SigExt and EDMD or GPR, regardless of the representation considered. For visual inspection, the TF representation of SST resulting from the BoundEffRed strategy is shown in the bottom-right panel of Fig. 1, where SST is applied to the portion of PPG displayed in Fig. 8. It produces a significant improvement in the quality of the SST near boundaries. Indeed, the blurring visible when zooming on the right boundary of the SST has almost vanished. Real-time tracking of the instantaneous frequencies contained in the measured signal is therefore greatly facilitated.

To further evaluate the influence of the noise level on the performance of BoundEffRed, we artificially add Gaussian noise to the measured PPG signal. It is thus an additional noise to the measurement noise actually contained in the signal. Fig. 9 shows the average performance index of BoundEffRed for different values of the Signal-to-Noise Ratio (SNR). Whatever the representation considered, BoundEffRed is relatively robust to noise, as long as the SNR remains above 10 dB. Below this level, the reduction of boundary effects gradually deteriorates.

E. Real-Time Implementation

We simulate the real-time implementation of the SST from BoundEffRed applied to the PPG described in Section IV-D3, subsampled by a factor of 2. The test is performed on a 2-Core Intel Core i5 CPU running at 1.7 GHz and 7.7 GB of RAM.

For this signal, sampled at f s = 65.5 Hz, the suitable window is 8 seconds long. Therefore, the extension must be over 4 seconds, i.e., L = 250. The forecasting step in each iteration of BoundEffRed then takes no more than t forecast = 46 ms. Denote by H ≥ 1 the hop size in samples, that is the number of samples between successive columns of the SST. Besides, since the support of the analysis window covers 2L samples, the update of the SST requires the calculation of 2L/H new columns of the SST. In this example, the frequency dimension of the SST being 512, the computational time of one column is t SST = 2.08 ms, on average. A general rule to determine the acceptable values of H for real-time implementation is verifying that the computational time to update the boundary-free TF representation is smaller than the lag between H samples; that is,

t forecast + L H t SST < H f s .
In our example, taking H ≥ 8 samples is sufficient to ensure the feasibility of real-time implementation. This value is reasonable because it allows a maximum overlap of 98.4% of the window length. SST variations are therefore perceived gradually. Combined with the low latency of the algorithm, this allows for smooth monitoring of the real-time SST. A video clip illustrating this is available on GitHub, at https://github.com/AdMeynard/ BoundaryEffectsReduction/tree/master/Animations.

V. CONCLUSION

In this paper, we propose an algorithm, named Bound-EffRed, for the real-time reduction of boundary effects in TF representations. This method is based on an extension of the signal obtained by a simple-minded and numerically efficient forecasting. We have shown theoretically that the chosen dynamic model is sufficient to extend signals formed by a sum of sine waves. Moreover, the low computational time allows us to switch to a real-time implementation of BoundEffRed, unlike other existing forecasting methods. The numerical results also confirmed the robustness to noise of BoundEffRed, as well as its ability to be applied to many TF representations. Additional applications to ECG and PPG are included in the Supplementary Materials.

Various improvements can be considered to make the algorithm more robust. In particular, we have noticed (see Fig. 5) that when the regular oscillations of the observed signal break, the forecasting step is no longer relevant, and only slows down the calculation of the TF representation. A preliminary step should then be added to the algorithm to detect signal activity and disable the forecasting step when possible. More fundamentally, one can also consider accelerating the computational time by optimizing the forecasting step to improve the realtime performance of BoundEffRed. Indeed, each new forecast requires the inversion of the matrix XX T of size M × M. However, the matrix X used to forecast x N+H at the current iteration differs from the one used at the previous iteration to forecast x N by only H columns. It then seems natural to take inspiration from the work of Strobach [START_REF] Strobach | Square hankel svd subspace tracking algorithms[END_REF], generalized by Badeau et al. [START_REF] Badeau | Sliding window adaptive svd algorithms[END_REF], which proposes a fast algorithm for the singular value decomposition of successive data matrices taking the same form as X. Their algorithm is limited to the case where H = 1. Developing an extension of this algorithm to variations of H > 1 columns would then allow us to efficiently update the matrix XX T -1 . These avenues will be explored in our future work.

Supplementary Materials for "An Efficient Forecasting Approach to Reduce Boundary Effects in Real-Time Time-Frequency Analysis" Adrien Meynard, Hau-Tieng Wu I. PROOF OF THEOREM 1 A. Preliminaries 1) Notations: Consider z ∈ R N the deterministic signal defined in [START_REF] Auger | Time-frequency reassignment and synchrosqueezing: An overview[END_REF], and denote by

z k ∈ R M the subsignal such that z k [m] = z[N -K -M + k + m] , ∀m ∈ {0, . . . , M -1} , k ∈ {0, . . . , K} .
Define Z ∈ R M×K and Z ∈ R M×K , the matrices such that

Z ∆ = z 0 • • • z K-1 , Z ∆ = z 1 • • • z K .
Let D ∈ R M×M be the matrix defined by D[m, m ] = δ m+1,m . Recall the model (13). Based on the definition of matrices X and Y, we have:

1 K XX T = 1 K ZZ T + σ 2 I ∆ =S (0) +E (0) (26) 1
K YX T = 1 K Z Z T + σ 2 D ∆ =S (1)
+E (1) , [START_REF] Chang | Training and testing low-degree polynomial data mappings via linear svm[END_REF] where

E (a) ∆ = σE (a) 1 + σ 2 E (a)
2 ,

E (a) 1 [m, m ] = 1 K K-1 ∑ k=0 z[N 0 + m + a + k]w[N 0 + m + k] + w[N 0 + m + a + k]z[N 0 + m + k] ,
and

E (a) 2 [m, m ] = 1 K K-1 ∑ k=0 w[N 0 + m + a + k]w[N 0 + m + k] -δ (m+a)m ,
with a ∈ {0, 1}. We call E (0) and E (1) error matrices because:

E{E (0) } = E{E (0) 1 } = E{E (0) 2 } = 0 E{E (1) } = E{E (1) 1 } = E{E (1)
2 } = 0 . Thus, the random matrix Ã, defined in equation [START_REF] Flandrin | Time-frequency/time-scale analysis[END_REF], is expressed in function of the above-defined matrices as: Ã = S (1) + E (1) S (0) + E (0) -1 .

Define A 0 the deterministic matrix such that A. Meynard is the corresponding author (e-mail: adrien.meynard@duke.edu).

A 0 ∆ = S (1) S (0) -1 . (28)
We denote by α ( ) 0 the last row of A 0 . As a result, for ∈ N * , the error vector h ( ) defined by h

( ) ∆ = α ( ) -α ( ) 0 satisfy the equation h ( ) = e T M Ã -A 0 = e T M
(S (1) + E (1) )(S (0) + E (0) )

-1 -A 0 . (29)
The randomness of h ( ) completely comes from the error matrices. Besides, notice that the first M -1 rows in E (1) equal to the last M -1 rows of E (0) . We gather all sources of randomness into a vector g ∈ R M(M+1) , defined as (1) , [START_REF] Hyndman | Forecasting: principles and practice[END_REF] where "vec" denotes the vectorization operator, that concatenates the columns of a given matrix on top of one another. Then, we can write h ( ) as h ( ) = f ( ) (g) where f ( ) is a deterministic function such that:

g = vec E (0) e T M E
f ( ) : R M(M+1) → R M g → h ( ) . (31)
In the following paragraph, we provide some useful lemmas to prove Theorem 1.

2) Lemmas:

Lemma 1 (Expressions of A 0 and S (0) -1 ). Let S (0) be the M × M matrix defined in [START_REF] Hall | Theory for penalised spline regression[END_REF]. Let A 0 the M × M matrix defined in [START_REF] Marron | Exact risk analysis of wavelet regression[END_REF]. Assume the deterministic signal z takes the form (11), and the observed noisy signal takes the form [START_REF] Chui | Numerical Methods in Approximation Theory[END_REF]. Then, the inverse of the matrix S (0) is given by

S (0) -1 [m, m ] = 1 σ 2 δ m,m - 2 Mσ 2 J ∑ j=1 1 1 + 4σ 2 MΩ 2 j cos 2π p j m -m M , (32)
and the matrix A 0 is given by

A 0 [m, m ] = δ m+1,m + 2 M J ∑ j=1 1 1 + 4σ 2 MΩ 2 j cos 2π p j m M δ m+1,M . (33)
Let • max denote the maximum norm of a matrix, i.e., M max = max n,n |M[n, n ]|. Then,

S (0) -1 max ≤ 1 σ 2 1 + 2J M , (34)
A 0 max ≤ max 1, 2J M . (35)
Proof. It follows from the signal model ( 11) that the matrices S (0) and S (1) take the following form:

S (a) [m, m ] = σ 2 δ (m+a) ,m + J ∑ j,j =1 Ω j Ω j K K-1 ∑ k=0 cos 2π f j f s (N 0 + m + a + k) + ϕ j cos 2π f j f s (N 0 + m + k) + ϕ j = σ 2 δ (m+a) ,m + J ∑ j=1 Ω 2 j 2K K-1 ∑ k=0 cos 2π f j f s (m + a -m ) + cos 2π f j f s (2k + m + a + m + 2N 0 ) = σ 2 δ (m+a) ,m + J ∑ j=1 Ω 2 j 2 cos 2π f j f s (m + a -m ) + Ω 2 j 2K K-1 ∑ k=0 cos 2π f j f s (2k + m + a + m + 2N 0 ) =0 because f j fs = p j K = σ 2 δ (m+a) ,m + J ∑ j=1 Ω 2 j 2 cos 2π f j f s (m + a -m ) . (36)
Thus, S (0) is a circulant matrix, and is therefore diagonalizable in the Fourier basis:

S (0) = UΛ (0) U * ,
where

U[m, m ] = 1 √ M e -2iπmm /M and Λ (0) = diag(λ (0) 0 , . . . , λ (0)
M-1 ) with λ (0) m = σ 2 + J ∑ j=1 Ω 2 j 2 M-1 ∑ q=0 cos 2π f j f s q e -2iπqm/M = σ 2 + M 4 J ∑ j=1 Ω 2 j (δ m,p j + δ m,M-p j ) . Therefore, S (0)
-1

= UΛ (0) -1 U * ,
which leads to

S (0) -1 [m, m ] = 1 σ 2 δ m,m - 2 Mσ 2 J ∑ j=1 1 1 + 4σ 2 MΩ 2 j cos 2π p j m -m M . (37)
Directly, we have:

S (0) -1 [m, m ] = 1 σ 2 δ m,m - 2 M J ∑ j=1 1 1 + 4σ 2 MΩ 2 j cos 2π p j m -m M ≤ 1 σ 2   1 + 2 M J ∑ j=1 1 1 + 4σ 2 MΩ 2 j    ≤ 1 σ 2 1 + 2J M .
Thus,

S (0) -1 max ≤ 1 σ 2 1 + 2J M .
Furthermore, combining equations ( 36) and ( 37), we have

A 0 [m, m ] = M-1 ∑ q=0 S (1) [m, q]S (0) -1 [q, m ] = δ m+1,m + 2 M J ∑ j=1 1 1 + 4σ 2 MΩ 2 j cos 2π p j m M δ m+1,M .
Directly, we have:

A 0 [m, m ] ≤      1 if m < M -1 , 2 M ∑ J j=1 1 1 + 4σ 2 MΩ 2 j if m = M -1 .
Thus,

A 0 max ≤ max 1, 2J M .
Lemma 2 (Moments of g). Let g ∈ R M(M+1) be the random vector defined in [START_REF] Hyndman | Forecasting: principles and practice[END_REF]. Assume the deterministic signal z takes the form (11), and the observed noisy signal takes the form (13). Then, as K → ∞, the second-order moments of g are bounded as follows:

E{g[r]g[r ]} ≤ 1 K C 2 z σ 2 + 2σ 4 + o 1 K , ∀(r, r ) ∈ {0, . . . , M(M + 1) -1} 2 , ( 38
)
where C z ∆ = 2 ∑ J j=1 Ω j . Besides, higher-order moments of g behave as o 1 K .

Proof. In the following, for all r ∈ {0, . . . , M(M + 1) -1}, we denote g

[r] = σg 1 [r] + σ 2 g 2 [r],
where

g 1 [r] = E (a r ) 1 [m r , m r ] = 1 K K-1 ∑ k=0 z k [m r + a r ]w k [m r ] + w k [m r + a r ]z k [m r ] , g 2 [r] = E (a r ) 2 [m r , m r ] = 1 K K-1 ∑ k=0 w k [m r + a r ]w k [m r ] -δ m r +a r , m r ,
and m r , m r , a r are the corresponding coordinates of the matrices associated with r through the vectorization operation [START_REF] Hyndman | Forecasting: principles and practice[END_REF]. Thus, order-two moments of this random vector is split as follows:

E{g[r]g[r ]} = σ 2 E{g 1 [r]g 1 [r ]} + σ 3 E{g 1 [r]g 2 [r ]} + σ 3 E{g 2 [r]g 1 [r ]} + σ 4 E{g 1 [r]g 1 [r ]} . ( 39
)
By definition of the signal z (see equation ( 11)), we have |z[n]| ≤ ∑ J j=1 Ω j , for all n ∈ N. Thus, by a direct bound, we have

E{g 1 [r]g 1 [r ]} = 1 K 2 E K-1 ∑ k,k =0 z k [m r + a r ]z k [m r + a r ]w k [m r ]w k [m r ] + z k [m r ]z k [m r + a r ]w k [m r + a r ]w k [m r ] + z k [m r + a r ]z k [m r ]w k [m r ]w k [m r + a r ] + z k [m r ]z k [m r ]w k [m r + a r ]w k [m r + a r ] ≤ 1 K 2 J ∑ j=1 Ω j 2 K-1 ∑ k,k =0 E w k [m r ]w k [m r ] + E w k [m r + a r ]w k [m r ] + E w k [m r ]w k [m r + a r ] + |E {w k [m r + a r ]w k [m r + a r ]}| (40)
Besides, since w is a white noise,

1 K 2 K-1 ∑ k,k =0 E w k [m r ]w k [m r ] = 1 K 2 K-1 ∑ k,k =0 δ k+m r , k +m r = 1 K 2 K -m r -m r .
Moreover, 0 ≤ m r -m r ≤ M -1. Thus,

1 K - M -1 K 2 ≤ 1 K 2 K-1 ∑ k,k =0 E w k [m r ]w k [m r ] ≤ 1 K . Therefore, 1 K 2 K-1 ∑ k,k =0 E w k [m r ]w k [m r ] = 1 K + o 1 K .
Similar calculations lead to the same results for the other three terms making up the sum [START_REF] Strobach | Square hankel svd subspace tracking algorithms[END_REF]. Therefore, we have:

E{g 1 [r]g 1 [r ]} ≤ J ∑ j=1 Ω j 2 4 K + o 1 K ≤ C 2 z K + o 1 K . ( 41
)
Besides, since odd-order moments of a zero-mean multivariate Gaussian random vector are zero, we have:

E{g 1 [r]g 2 [r ]} = 1 K 2 K-1 ∑ k,k =0 z k [m r + a r ]E w k [m r ]w k [m r + a r ]w k [m r ] + z k [m r ]E w k [m r + a r ]w k [m r + a r ]w k [m r ] -δ m r +a r ,m r E{g 1 [r]} = 0 . (42) Similarly, E{g 2 [r]g 1 [r ]} = 0 . ( 43
)
Besides, by a direct calculation, we have:

E{g 2 [r]g 2 [r ]} = 1 K 2 K-1 ∑ k,k =0 E w k [m r + a r ]w k [m r ]w k [m r + a r ]w k [m r ] -δ m r +a r ,m r 1 K K-1 ∑ k =0 E w k [m r + a r ]w k [m r ] -δ m r +a r ,m r 1 K K-1 ∑ k=0 E w k [m r ]w k [m r ] + δ m r +a r ,m r δ m r +a r ,m r = 1 K 2 K-1 ∑ k,k =0 E w k [m r + a r ]w k [m r ]w k [m r + a r ]w k [m r ] -δ m r +a r ,m r δ m r +a r ,m r (44)
Moreover, using the results of the Isserlis' theorem [1] to express fourth-order moments of a Gaussian random vector as a function of its second-order moments, we expand sum (44) as follows:

E{g 2 [r]g 2 [r ]} = 1 K 2 K-1 ∑ k,k =0 E w k [m r + a r ]w k [m r ] E w k [m r + a r ]w k [m r ] -δ m r +a r ,m r δ m r +a r ,m r + 1 K 2 K-1 ∑ k,k =0 E {w k [m r + a r ]w k [m r + a r ]} E w k [m r ]w k [m r ] + 1 K 2 K-1 ∑ k,k =0 E w k [m r + a r ]w k [m r ] E w k [m r ]w k [m r + a r ] = 1 K 2 K-1 ∑ k,k =0 δ k+m r +a r ,k +m r +a r δ k+m r,k +m r + δ k+m r +a r ,k +m r δ k+m r ,k +m r +a r = 1 K 2 δ m r -m r ,m r +a r -m r -a r (K -|m r -m r |) + δ m r +a r -m r ,m r -m r -a r (K -|m r + a r -m r |) = 1 K δ m r -m r ,m r +a r -m r -a r + δ m r +m r ,m r +a r +m r +a r + o 1 K .
Therefore,

E{g 2 [r]g 2 [r ]} ≤ 2 K + o 1 K . ( 45
)
Thus, combining results (41), (42), ( 43) and (45) into expression [START_REF] Goldberger | Physiobank, physiotoolkit, and physionet: Components of a new research resource for complex physiologic signals[END_REF] gives the following bound:

E{g[r]g[r ]} ≤ 1 K C 2 z σ 2 + 2σ 4 + o 1 K .
Concerning higher-order moments, let T ≥ 3 denote the order of the moment defined by E{∏ T θ=1 g[r θ ]}. Thus,

E T ∏ θ=1 g[r θ ] = E T ∏ θ=1 σg 1 [r θ ] + σ 2 g 2 [r θ ] = 1 K T E T ∏ θ=1 K-1 ∑ k=0 ρ θ,k w[k + m r θ ], w[k + m r θ + a r θ ] = 1 K T K-1 ∑ k 1 =0 • • • K-1 ∑ k T =0 E T ∏ θ=1 ρ θ,k θ w[k θ + m r θ ], w[k θ + m r θ + a r θ ] , (46)
where

ρ θ,k (u, v) = σz k [m r θ + a r θ ]u + σz k [m r θ ]v + σ 2 uv -δ m r θ +a r θ , m r θ .
Thus,

E T ∏ θ=1 g[r θ ] ≤ ν K K T C T ,
where

C T = max (k 1 ,...,k T ) E T ∏ θ=1 ρ θ,k θ (w[k θ + m r θ ], w[k θ + m r θ + a r θ ]) ,
and ν K is the number of nonzero terms in the sum (46). Note that C T is independent of K (but depends on z and σ). The behavior of the order T moment in function of K is therefore only determined by the ratio

ν K K T .
Let us bound ν K . Fix k 1 . For each of the other indexes of summation k θ (θ ∈ {2, . . . , T}), there are four values that make ρ

θ,k θ (w[k θ + m r θ ], w[k θ + m r θ + a r θ ]) not independent of ρ 1,k 1 (w[k 1 + m r 1 ], w[k 1 + m r 1 + a r 1 ]
). Indeed, since w is a white noise these quantities are independent except when k θ is such that:

k θ = k 1 + m r 1 -m r θ k θ = k 1 + m r 1 + a r 1 -m r θ k θ = k 1 + m r 1 -m r θ -a r θ k θ = k 1 + m r 1 + a r 1 -m r θ -a r θ .
Consequently, for each value of k 1 there exist at least (K -4) T-1 combinations of (k 2 , . . . , k T ) where we have

E T ∏ θ=1 ρ θ,k θ (w[k θ + m r θ ], w[k θ + m r θ + a r θ ]) = E ρ 1,k 1 (w[k 1 + m r 1 ], w[k 1 + m r 1 + a r 1 ]) × E T ∏ θ=2 ρ θ,k θ (w[k θ + m r θ ], w[k θ + m r θ + a r θ ]) = 0 , because E ρ 1,k 1 (w[k 1 + m r 1 ], w[k 1 + m r 1 + a r 1 ]) = 0. Therefore, at least K(K -4)
T-1 of the sum (46) are zero.

Because T ≥ 3, we develop similar arguments on k 2 and k 3 to determine other cases where this correlation term vanishes. We subtract these cases to K T , the total number of combinations of (k 1 , . . . , k T ) to obtain the following maximum bound on the number of nonzero terms in the sum (46):

ν K ≤ K T -3K(K -4) T-1 + 3K(K -4)(K -8) T-2 -K(K -4)(K -8)(K -12) T-3 .
Thus,

E T ∏ θ=1 g[r θ ] ≤ C T K T -3K(K -4) T-1 + 3K(K -4)(K -8) T-2 -K(K -4)(K -8)(K -12) T-3 K T ≤ C T 1 -3 1 - 4 K T-1 + 3 1 - 4 K 1 - 8 K T-2 -1 - 4 K 1 - 8 K 1 - 12 K T-3 ≤ C T 1 -3 + 12(T -1) K + 3 - 12 K - 24(T -2) K -1 + 4 K + 8 K + 12(T -3) K + o 1 K ≤ C T o 1 K .
Therefore,

E T ∏ θ=1 g[r θ ] = o 1 K , ∀T ≥ 3 .
Lemma 3 (Bounds on the derivatives of f ( ) at the origin). Let f ( ) : R M(M+1) → R M denote the multivariate function defined in [START_REF] Schmid | Dynamic mode decomposition of numerical and experimental data[END_REF]. Assume the deterministic signal z takes the form [START_REF] Auger | Time-frequency reassignment and synchrosqueezing: An overview[END_REF], and the observed noisy signal takes the form [START_REF] Chui | Numerical Methods in Approximation Theory[END_REF]. Then, the first-order derivatives of f ( ) at the origin are bounded as follows:

∂ f ( ) m ∂g[r] g=0 ≤ d 1,z,M, σ 2 , ∀m ∈ {0, . . . , M -1} , r ∈ {0, . . . , M(M + 1) -1} , ( 47
)
where

d 1,z,M, ∆ = (2 + ( -2)M) M -1 max 1, 2J M -1 1 + max 1, 2J M 1 + 2J M .
Besides, the second-order derivatives of f ( ) at the origin are bounded as follows:

∂ 2 f ( ) m ∂g[r]∂g[r ] g=0 ≤ d 2,z,M, σ 4 , ∀m ∈ {0, . . . , M -1} , (r, r ) ∈ {0, . . . , M(M + 1) -1} 2 , ( 48
)
where

d 2,z,M, ∆ = 1 + 2J M 2 max 1, 2J M -2 1 + max 1, 2J M d 2,M, + (d 2,M, + 2d 2,M, ) max 1, 2J M ,
and d 2,M, and d 2,M, are only depending on M and .

Proof. Concerning the first-order derivative, from (29), we have:

∂ f ( ) ∂g[r] = e T M ∂ Ã ∂g[r] = -1 ∑ λ=0 e T M Ãλ ∂ Ã ∂g[r] Ã -1-λ . Thus, ∂ f ( ) ∂g[r] g=0 = -1 ∑ λ=0 e T M A λ 0 ∂ Ã ∂g[r] g=0 A -1-λ 0 . ( 49
)
Furthermore, (1) ∂g[r] S (0) + E (0) -1 + S (1) + E (1) ∂

∂ Ã ∂g[r] = ∂E
S (0) + E (0) -1 ∂g[r] =      J m r ,m r S (0) + E (0) -1 if a r = 1,
(1δ m r ,0 )J m r -1,m r + S (1) + E (1)

S (0) + E (0) -1 J m r ,m r S (0) + E (0) -1 else,
where

J m r ,m r ∈ R M×M is the matrix such that J m r ,m r [m, m ] ∆ = δ m,m r δ m ,m r . Thus, ∂ Ã ∂g[r] g=0 =    J m r ,m r S (0) -1 if a r = 1, ( 1
-δ m r ,0 )J m r -1,m r + A 0 J m r ,m r S (0) -1 else. Then, ∂ Ã ∂g[r] g=0 max ≤ (1 + A 0 max ) S (0) -1 max .
Given bounds [START_REF] Lin | ConceFT for Time-Varying Heart Rate Variability Analysis as a Measure of Noxious Stimulation During General Anesthesia[END_REF] and [START_REF] Korda | Linear predictors for nonlinear dynamical systems: Koopman operator meets model predictive control[END_REF], we have that:

∂ Ã ∂g[r] g=0 max ≤ 1 + max 1, 2J M 1 + 2J M 1 σ 2 . ( 50
)
Besides given expression (49), for all r ∈ {0, . . . , M(M + 1) -1}, we have:

∂ f ( ) m ∂g[r] g=0 ≤ 2M A -1 0 max ∂ Ã ∂g[r] g=0 max + M 2 -2 ∑ λ=1 A λ 0 max ∂ Ã ∂g[r] g=0 max A -1-λ 0 max ≤ (2 + ( -2)M) M -1 A 0 -1 max ∂ Ã ∂g[r] g=0 max .
Therefore, given bounds ( 35) and (50), we have:

∂ f ( ) m ∂g[r] g=0 ≤ d 1,z,M, σ 2 .
Concerning the second-order derivative, we have:

∂ 2 f ( ) ∂g[r]∂g[r ] = -1 ∑ λ=0 e T M ∂ Ãλ ∂g[r ] ∂ Ã ∂g[r] Ã -1-λ + e T M Ãλ ∂ 2 Ã ∂g[r]∂g[r ] Ã -1-λ + e T M Ãλ ∂ Ã ∂g[r] ∂ Ã -1-λ ∂g[r ] = -1 ∑ λ=1 λ-1 ∑ p=0 e T M Ãp ∂ Ã ∂g[r ] Ãλ-1-p ∂ Ã ∂g[r] Ã -1-λ + -1 ∑ λ=0 e T M Ãl ∂ 2 Ã ∂g[r]∂g[r ] Ã -1-λ + -2 ∑ λ=0 -λ-2 ∑ p=0 e T M Ãλ ∂ Ã ∂g[r] Ãp ∂ Ã ∂g[r ] Ã -λ-2-p .
Thus,

∂ 2 f ( ) ∂g[r]∂g[r ] g=0 = -1 ∑ λ=1 λ-1 ∑ p=0 e T M A p 0 ∂ Ã ∂g[r ] g=0 A λ-1-p 0 ∂ Ã ∂g[r] g=0 A -1-λ 0 + -1 ∑ λ=0 e T M A λ 0 ∂ 2 Ã ∂g[r]∂g[r ] g=0 A -1-λ 0 + -2 ∑ λ=0 -2-λ ∑ p=0 e T M A λ 0 ∂ Ã ∂g[r] g=0 A p 0 ∂ Ã ∂g[r ] g=0 A -λ-2-p 0 . (51)
Besides, the second-order derivative of the matrix à is given by

∂ 2 Ã ∂g[r]∂g[r ] =                                      0 if a r = 1 and a r = 1, J m r ,m r S (0) + E (0) -1 J m r ,m r S (0) + E (0) -1 if a r = 1 and a r = 0, J m r ,m r S (0) + E (0) -1 J m r ,m r S (0) + E (0) -1 if a r = 0 and a r = 1, ( 1
-δ m r ,0 )J m r -1,m r + S (1) + E (1) S (0) + E (0) -1 J m r ,m r × S (0) + E (0) -1 J m r ,m r S (0) + E (0) -1 + (1 -δ m r ,0 )J m r -1,m r + S (1) + E (1) S (0) + E (0) -1 J m r ,m r × S (0) + E (0) -1 J m r ,m r S (0) + E (0) -1 else.
Thus, the second-order derivative of the matrix à at the origin is such that

∂ 2 Ã ∂g[r]∂g[r ] g=0 =                    0 if a r = 1 and a r = 1, J m r ,m r S (0) -1 J m r ,m r S (0) -1 if a r = 1 and a r = 0, J m r ,m r S (0) -1 J m r ,m r S (0) -1 if a r = 0 and a r = 1, ( 1
-δ m r ,0 )J m r -1,m r + A 0 J m r ,m r S (0) -1 J m r ,m r S (0) -1 + (1 -δ m r ,0 )J m r -1,m r + A 0 J m r ,m r S (0) -1 J m r ,m r S (0) -1 else.
Then,

∂ 2 Ã ∂g[r]∂g[r ] g=0 max ≤      S (0) -1 2 max if a r = 1 or a r = 1 , 2 (1 + A 0 max ) S (0) -1 2 max else ≤ 2 1 + max 1, 2J M 1 + 2J M 2 1 σ 4 . (52)
Returning to equation ( 51), for all r, r ∈ {0, . . . , M(M + 1) -1} and m ∈ {0, . . . , M -1}, we have:

∂ 2 f ( ) m ∂g[r]∂g[r ] g=0 ≤ d 2,M, A 0 -2 max ∂ Ã ∂g[r] max ∂ Ã ∂g[r ] max + d 2,M, A 0 -1 max ∂ 2 Ã ∂g[r]g[r ] max ,
where d 2,M, and d 2,M, are only depending on M and . Besides, given results ( 28), ( 50) and (52), we have:

∂ 2 f ( ) m ∂g[r]∂g[r ] g=0 ≤ d 2,z,M, σ 4 .

B. Expression of the Bias µ.

By definition of the measurement noise, µ[n] = 0 when n ∈ I. Outside the measurement interval I, denote by the index such that n = N -1 + . Then, given that h ( ) = α ( )α ( ) 0 , we deduce from expression [START_REF] Chen | Symmetric extension methods for m-channel linear-phase perfect-reconstruction filter banks[END_REF] that

µ[n] = E{α ( ) }z K + σE{α ( ) w K } -z[n] = α ( ) 0 z K + E{h ( ) }z K + σE{h ( ) w K } -z[N -1 + ] ∆ = 1 [ ] + 2 [ ] + 3 [ ] , (53)
where

1 [ ] = α ( ) 0 z K -z[N -1 + ] , (54) 2
[ ] = E{h ( ) }z K , (55) 3
[ ] = σE{h ( ) w K } . (56)
Let us first determine an upper bound on

| 1 [n]|. Since α (1)
0 is the last row of A 0 , we deduce from the expression [START_REF] Daubechies | ConceFT: concentration of frequency and time via a multitapered synchrosqueezed transform[END_REF]

of A 0 that α (1) 0 [m] = 2 M J ∑ j=1 1 1 + 4σ 2 MΩ 2 j cos 2π p j m M . (57)
Besides, from equation ( 33), we also have

A 0 z K =      z[N -M + 1] . . . z[N -1] α (1) 0 z K      .
The upward-shift property is thus successively inducted when increases; that is,

A 0 z K =            z[N -M + ] . . . z[N -1] α (1) 0 z K . . . α ( ) 0 z K            .
Then, α ( ) 0 , the last row of A 0 follows the following recurrence relation:

α ( ) 0 z K = α (1) 0 Ã -1 0 z K = M- ∑ m=0 α (1) 0 [m]z[N -M + + m -1] + M-1 ∑ m=M-+1 α (1) 0 [m]α (m-M+ ) 0 z K . Hence, 1 [ ] = α ( ) 0 z K -z[N -1 + ] = M ∑ m=0 α (1) 0 [m]z[N -M + + m -1] -z[N -1 + ] + M-1 ∑ m=M-+1 α (1) 0 [m] α (m-M+ ) 0 z K -z[N -M + + m -1] = M ∑ m=0 α (1) 0 [m]z[N -M + + m -1] -z[N -1 + ] + M-1 ∑ m=M-+1 α (1) 0 [m] 1 [m -M + ] . Besides, equation (57) gives M ∑ m=0 α (1) 0 [m]z[N -M + + m -1] = J ∑ j,j =1 Ω j 2 M 1 1 + 4σ 2 MΩ 2 j M-1 ∑ m=0 cos 2π p j m M cos 2π p j N + m M + ϕ j =δ j,j M 2 cos(2π p j N M +ϕ j ) = J ∑ j=1 Ω j 1 1 + 4σ 2 MΩ 2 j cos 2π p j N M + ϕ j . Thus: | 1 [ ]| = J ∑ j=1 Ω j    1 1 + 4σ 2 MΩ 2 j -1    cos 2π p j N M + ϕ j + M-1 ∑ m=M-+1 α (1) 0 [m] 1 [m -M + ] ≤ J ∑ j=1 Ω j 1 1 + 4σ 2 MΩ 2 j -1 + 2J M -1 ∑ λ=1 | 1 [λ]| ≤ 4σ 2 M J ∑ j=1 1 Ω j + 2J M -1 ∑ λ=1 | 1 [λ]| . (58)
Then, by induction from the inequality (58), we have

| 1 [ ]| ≤ 4σ 2 M 1 + 2J M -1 J ∑ j=1 1 Ω j ∆ = c (1) σ 2 , ( 59
)
where

c (1) = 4 M 1 + 2J M -1 ∑ J j=1
1 Ω j . Note that c (1) is not depending on σ or K.

Let us now determine an upper bound on | 2 [ ]|. Since E{g[r]} = 0 and moments of order 3 and higher behave as o 1 K , a second-order Taylor expansion of h ( ) gives

E{h ( ) [m]} = 1 2 M(M+1)-1 ∑ r,r =0 ∂ 2 f ( ) m ∂g[r]∂g[r ] g=0 E{g[r]g[r ]} + o 1 K . ( 60
| 2 [ ]| ≤ C z M 3 (M + 1) 2 4 ∂ 2 f ( ) ∂g[r]∂g[r ] g=0 max 1 K C 2 z σ 2 + 2σ 4 + o 1 K ≤ 1 K c (2) 1 + c (2) 2 σ 2 + o 1 K , (61)
where

c (2) 1 ∆ = M 3 (M + 1) 2 2 C z d 2,z,M, , c (2) 2 ∆
= M 3 (M + 1) 2 4 C 3 z d 2,z,M, .
Let us now determine an upper bound on | 3 [ ]|. A second-order Taylor expansion of h ( ) gives

E{h ( ) w K } = M-1 ∑ m=0 M(M+1)-1 ∑ r=0 ∂ f ( ) m ∂g[r] g=0 E{g[r]w K [m]} + o 1 K (62)
Indeed, the third-order moments E{g

[r]g[r ]w K [m]}, behave as o 1 K . Besides, E{g[r]w K [m]} = σE{g 1 [r]w K [m]} + σ 2 E{g 2 [r]w K [m]} . (63)
Then,

|E{g 1 [r]w K [m]}| = 1 K K-1 ∑ k=0 z k [m r + a r ]E w k [m r ]w K [m] + z k [m r ]E {w k [m r + a r ]w K [m]} ≤ 2 K J ∑ j=1 Ω j = C z K . ( 64
) E{g 2 [r]w K [m]} = 1 K K-1 ∑ k,=0 E w k [m r + a r ]w k [m r ]w K [m] -δ m r +a r ,m r E {w K [m]} = 0 . ( 65
0 z K 2 + 2 α ( ) 0 z K E h ( ) z K + E h ( ) z K 2 + 2σα ( ) 0 E{w K h ( ) }z K + 2σα ( ) 0 z K E{h ( ) w K } + 2σE{h ( ) w K h ( ) }z K + σ 2 α ( ) 0 2 + σ 2 E h ( ) w K 2 + 2σ 2 α ( ) 0 E{w K h ( ) w K } -z[n] 2 -2z[n]µ[n] -µ[n] 2 = σ 2 α ( ) 0 2 -z[n] -α ( ) 0 z K + z[n] 2 + E h ( ) z K 2 + 2σα ( ) 0 E{w K h ( ) }z K + 2σE{h ( ) w K h ( ) }z K + σ 2 E h ( ) w K 2 + 2σ 2 α ( ) 0 E{w K h ( ) w K } ∆ = η 1 [n] + η 2 [n] + η 3 [n] + η 4 [n] + η 5 [n] + η 6 [n] + η 7 [n] ,
where

η 1 [n] = σ 2 α ( ) 0 2 , η 2 [n] = -z[n] -α ( ) 0 z K + z[n] 2 , η 3 [n] = E h ( ) z K 2 , η 4 [n] = 2σα ( ) 0 E{w K h ( ) }z K , η 5 [n] = σE{h ( ) w K h ( ) }z K , η 6 [n] = σ 2 E h ( ) w K 2 , η 7 [n] = 2σ 2 α ( )
0 E{w K h ( ) w K } . Let us now determine an upper bound on each of these terms.

First, since α ( ) 0 2 ≤ M A 0 2 max , we have:

η 1 [n] ≤ σ 2 M A 0 2 max ≤ σ 2 M A 0 2 max ≤ σ 2 M max 1, 2J M 2 . ( 69
)
Second, by definition of 2 and 3 (see expressions (55) and (56)), η 2 takes the following form:

η 2 [n] = ( 2 [n] + 3 [n]) 2 = o 1 K . ( 70
)
Third, second-order Taylor expansions of h ( ) give: ≤ M(M + 1)

d 1,z,M, σ 2 C z σ K + o 1 K (72)
Thus,

|η 4 [n]| ≤ M (M + 1)C 2 z d 1,z,M, max 1, 2J M 1 K + o 1 K . ( 73
)
Fifth and sixth, second-order Taylor expansions of h ( ) give

η 5 [n] = o 1 K , ( 74
)
η 6 [n] = o 1 K . ( 75
)
Seventh,

|η 7 [n]| ≤ 2σ 2 A 0 max M-1 ∑ m,m =0 E{h ( ) [m]w K [m]w K [m ]} ≤ 2σ 2 A 0 max M-1 ∑ m,m =0 M(M+1)-1 ∑ r=0 ∂ f ( ) m ∂g[r] g=0 E{g ( ) [r]w K [m]w K [m ]} + o 1 K . But, E{g ( ) [r]w K [m]w K [m ]} = σE{g ( ) 1 [r]w K [m]w K [m ]} + σ 2 E{g ( ) 2 [r]w K [m]w K [m ]} .
As before, since E{g ( )

1 [r]w K [m]w K [m ]
} is a third-order moment of a multivariate zero-mean Gaussian vector, it vanishes. And, E{g ( )

2 [r]w K [m]w K [m ]} = 1 K K-1 ∑ k=0 E w k [m r + a r ]w k [m r ]w K [m]w K [m ] -δ m r +a r ,m r E w K [m]w K [m ] = 1 K K-1 ∑ k=0
δ k+m r +a r ,K+m δ k+m r ,K+m + δ k+m r +a r ,K+m δ k+m r ,K+m ≤ 2 K .

III. APPLICATION TO A MULTICOMPONENT CARDIAC SIGNAL

We consider a cardiac signal, namely a photoplethysmogram (PPG) recording, sampled at 300 Hz. In addition to the cardiac cycle measurement, this signal contains a slow varying component, which is the respiratory component, known as respiration induced intensity variation (RIIV) [2]. The top of Fig. 11 displays an excerpt of the signal along with a 3-second-long extension obtained by the SigExt forecasting. The lower part of Fig. 11 shows the respiratory signal recording the concentration of CO 2 . This signal was recorded simultaneously with the PPG signal, and visually highlights the low-frequency respiratory component contained in the PPG signal. Indeed, the intervals where the concentration of CO 2 drops-highlighted by the bluish areas-coincide with the decreases in the PPG signal. Note also that the forecasting breaks the waveform of the oscillations because of its inability to forecast the high-frequency harmonics contained in the signal. Nevertheless, as long as the low-frequency harmonics of the components contained in the signal are preserved, the forecasting is sufficient to reduce boundary effects in the low-frequency TF domain.

The ordinary and boundary-free SSTs of this signal are displayed in Fig. 12. These TF analyses brings out both components-the fundamental frequency of the cardiac component, the most energetic one, indicated by the blue arrows, is located around 1.3 Hz, while the respiratory component, less visible, indicated by the green arrows, is located around 0.2 Hz and its multiples. Clearly, near the boundaries, BoundEffRed helps improve the quality of the TF representation. Besides, the performance index of this representation takes the value D = 0.491. This means that BoundEffRed has reduced the right-side boundary effects by about 51% with respect to the ordinary SST. This shows the ability of our algorithm to work on signals containing several nonstationary components.

  . The A. Meynard is with the Department of Mathematics, Duke University, Durham, NC, 27708 USA. H.-T. Wu is with the Department of Mathematics and Department of Statistical Science, Duke University, Durham, NC, 27708 USA; Mathematics Division, National Center for Theoretical Sciences, Taipei, Taiwan.

Fig. 1 .

 1 Fig. 1. A segment of PPG signal (top) and the right boundary of a TF representation determined by the SST without extension (bottom left), and the right boundary of a TF representation determined by the SST with the proposed boundary effect reduction algorithm by forecasting (bottom right). The window length for the SST is 12 seconds. Videos of these real-time SSTs are available online at https://github.com/AdMeynard/BoundaryEffectsReduction.

1 )

 1 Extend the signal by forecasting it. The aim is to use a simple dynamic model to predict the values taken by the measured signal outside the measurement interval. Then, once this operation is done, we have access to an extended signal defined on a larger interval I ∆ , where ∆ denotes the size of the extension on both boundaries of I.2) Run the local analysis tool on the extended signal.

  Real-time output: Signal representation F (N) x end while A. Step 1: Extension by Forecasting a) Dynamical Model: Establishing a dynamical model means determining the relation linking Y to X; that is, finding a function f so that

Fig. 2 .

 2 Fig. 2. Evolution of the experimental forecasting variance as a function of the noise variance for three different values the forecasting sample index .

3 Fig. 3 .

 33 Fig. 3. Evolution of the experimental forecasting variance in function of the dataset size for three different values the forecasting sample index .

Fig. 4 .

 4 Fig. 4. Extended respiratory signal (blue) obtained by the SigExt forecasting (first panel), the EDMD forecasting (second panel), the GPR forecasting (third panel), and the TBATS forecasting (fourth panel), superimposed with the ground truth signal (red-dashed line).

Fig. 5 .

 5 Fig. 5. Extensions of a segment of the respiratory signal (left) where SigExt is outperformed by the EDMD and Symmetric extensions. Corresponding SST (top-right) and boundary-free SST obtained with SigExt (bottom-right).

Fig. 6 .

 6 Fig. 6. Two-component respiratory signal (top, yellow), and the associated SigExt extension (blue); the simultaneously recorded ECG signal is below. The vertical magenta dotted lines indicate the cardiac cycles in the respiratory signal.

Fig. 7 .Fig. 8 .

 78 Fig. 7. Ordinary SST (top-left) and boundary-free SST (top-right) of a two-component respiratory signal. Middle and bottom: zoom on the areas delimited by the red dashed rectangles; these images are rotated 90 degrees clockwise. The window length for the SSTs is 20 seconds. The instantaneous frequency of the respiration is indicated by the green arrows, while the instantaneous frequency of the cardiogenic artifact is indicated by the blue arrows.

A

  . Meynard is with the Department of Mathematics, Duke University, Durham, NC, 27708 USA. H.-T. Wu is with the Department of Mathematics and Department of Statistical Science, Duke University, Durham, NC, 27708 USA; Mathematics Division, National Center for Theoretical Sciences, Taipei, Taiwan.

  ) Thus, given the bounds (38) on |E{g[r]g[r ]}| and (48) on the second derivative of f ( ) m , we have:

  ) Thus, combining results (64) and (65) into expression (63) gives the following bound:|E{g[r]w K [m]}| ≤ C z σ K . (66)Thus, given the bound (47) on the first derivative of f ( ) m , and from (62) we have:| 3 [ ]| ≤ σM 2 (M + 1) d 1,z,M, c (3) ∆ = M 2 (M + 1)d 1,z,M, C z .Thus, combining bounds (59) on 1 , (61) on 2 and (67) on 3 gives the following bound of the bias:|µ[n]| ≤ c (1) σ 2 + of the Covariance γ.Outside the measurement interval I, denote by the index such that n = N -1 + . Then, given that h ( ) = α ( )α

η 3 [|η 4 1 ∑

 341 [n]| ≤ σC z A 0 max Mm,m =0 E{h ( ) [m]w K [m ]} But, E{h ( ) [m]w K [m ]} ≤ ∂g[r] g=0 E{g[r]w K [m ]} + o 1 K

Fig. 10 .

 10 Fig. 10. Extended ECG (blue) obtained by the SigExt forecasting (first panel), the EDMD forecasting (second panel), the GPR forecasting (third panel), and the TBATS forecasting (fourth panel), superimposed with the ground truth signal (dashed red).

Fig. 11 .

 11 Fig. 11. PPG signal (top, yellow), and the associated SigExt extension (blue). The simultaneously recorded concentration of CO 2 is below. The bluish areas show the synchrony between the drops in CO 2 concentration and the decreases in the PPG signal.

Fig. 12 .

 12 Fig. 12. Ordinary SST (left) and boundary-free SST (right) of the PPG signal. The window length for the SSTs is 18 seconds. The green arrows indicate the instantaneous frequency of the respiratory component. The blue arrows, on the other hand, indicate the instantaneous frequency of the blood pressure component. The red dashed boxes highlight the areas where boundary effects may appear.

TABLE II RESPIRATORY

 II SIGNAL: PERFORMANCE OF THE EXTENSION METHODS AND THE ASSOCIATED BOUNDARY-FREE TF REPRESENTATIONS

	Extension method	MSE	Performance index D (mean ± SD) STFT SST RS ConceFT
	SigExt Symmetric EDMD GPR	0.292 ±4.438 0.044 ±0.111 0.026 ±0.112 0.331 ±4.858	0.370 ±0.623 ±0.436 ±0.879 0.408 0.866 1.162 1.173 1.022 ±0.893 ±0.886 ±0.281 0.359 0.422 0.828 ±0.266 ±0.282 ±0.248 0.391 0.411 0.897 ±0.853 ±0.406 ±1.140	0.423 ±0.344 1.144 ±0.579 0.449 ±0.296 0.430 ±0.364

TABLE IV ECG

 IV : PERFORMANCE OF THE BOUNDARY-FREE TF REPRESENTATIONS ACCORDING TO THE EXTENSION METHOD

	Extension method	Median performance index D STFT SST ConceFT RS
	SigExt Symmetric EDMD GPR	0.584 1.199 0.538 0.639	0.630 1.354 0.558 0.588	0.462 1.427 0.496 0.485	0.642 0.943 0.714 0.616

Thus,

To conclude, we combine the expressions (69), (70), (71), ( 73), (74), (75), and (76) to determine an upper bound on the variance γ [n, n]. It follows:

a) If n ≥ N: When n ≥ N and n ≥ N, applying the Cauchy-Schwarz inequality, we obtain the following bound on the covariance γ[n, n]:

A first-order Taylor expansion of this bound as

b) If n ∈ I: When n ≥ N and n ∈ I, equation [START_REF] Williams | A discrete wavelet transform without edge effects using wavelet extrapolation[END_REF] shows us that:

where

Besides, thanks to the bound (72) we have:

Besides,

Besides, identically to the bound (76) on η 7 , we obtain

Finally, we combine expressions (77), (78), and (79) to determine an upper bound on the variance γ[n, n ]. It follows:

II. APPLICATION TO AN ELECTROCARDIOGRAM

We provide here an additional implementation of BoundEffRed, applied to an electrocardiogram (ECG) dataset. The dataset is constructed from a 500-second-long ECG, sampled at f s = 200 Hz, cut into 10 segments of 50 seconds each. Fig. 10 depicts the right boundary of one of these subsignals, together with the 6-second extensions estimated by SigExt (first panel), or EDMD (second panel), GPR (third panel), or TBATS (fourth panel). These extensions are superimposed to the ground-truth extension, plotted in red. The sharp and spiky ECG patterns make the AHM model too simplistic to describe this type of signal. Consequently, the forecast produced by SigExt is moderately satisfactory. Note that TBATS is the only one that seems to accurately capture the locations of QRS complexes after 37 seconds. We will explore this long-term prediction capability in the future work.

Table IV contains the median performance index D of the boundary-free TF representations, over the N subsignals evaluated, according to the extension method. As a result of the fair quality of the forecasts, the reduction of boundary effects is less significant than for PPG signal. Nevertheles, the results show that BoundEffRed has the same efficiency when the SigExt extension, the EDMD extension or the GPR extension is chosen. Indeed, t-tests performed under the null hypothesis that the mean are equals, with a 5% significance level, show no statistical significant difference between SigExt and EDMD or GPR, regardless of the representation considered. This justifies the choice of SigExt for real-time implementation.