Reza Sameni
email: rsameni@shirazu.ac.ir

Online Filtering Using Piecewise Smoothness Priors: Application to Normal and Abnormal Electrocardiogram Denoising

Keywords: Tikhonov regularization, forward-backward filtering, electrocardiogram filtering, wavelet denoising, Lipschitz regularity, Wiener smoothing

In this work, a block-wise extension of Tikhonov regularization is proposed for denoising smooth signals contaminated by wide-band noise. The proposed method is derived from a constrained least squares problem in two forms: 1) a block-wise fixed-lag smoother with smooth inter-block transitions applied in matrix form, and 2) a fixed-interval smoother applied as a forward-backward zero-phase filter. The filter response is maximally flat and monotonically decreasing, without any ripples in its pass-band. The method is also extended to smoothness of multiple smoothness orders, and its relationship with Lipschitz regularity and block-wise Wiener smoothing is also studied.

The denoising of normal and abnormal electrocardiogram (ECG) signals in different stationary and non-stationary noise levels is studied as case study. While most ECG denoising techniques benefit from the pseudo-periodicity of the ECG, the developed technique is merely based on the smoothness assumption, which makes it a powerful method for both normal and abnormal ECG. The performance of the method is assessed by Monte-Carlo simulations over three standard normal and abnormal ECG databases of different sampling rates, in comparison with bandpass filtering, wavelet denoising with various parameters, and Savitzky-Golay filters using Stein's unbiased risk estimate shrinkage scheme.

Introduction

Signal denoising is a prominent branch of signal processing, which is based on a mixing model and prior assumptions on the signal and noise structure. While the most widely studied mixing models are additive mixtures of signal and noise, the signal/noise priors are totally subjective. Deterministic vs. stochastic, temporal vs. frequency domain or spatial, static vs. dynamical, linear vs. periodic, are among the different priors used for designing denoising schemes [START_REF] Candy | Model-Based Signal Processing[END_REF].

In absence of reliable priors regarding the signal and noise structure, generic filtering techniques such as bandpass filtering, wavelet denoising [START_REF] Donoho | Adapting to unknown smoothness via wavelet shrinkage[END_REF], Tikhonov This is the author's version of the manuscript accepted for publication in Elsevier Signal Processing, October 2016. The paper can be found in its final format at: http://dx.doi.org/10.1016/j.sigpro.2016. [START_REF] Zhang | Bioelectric signal detrending using smoothness prior approach[END_REF].019 regularization and smoothing splines [START_REF] Golub | Tikhonov regularization and total least squares[END_REF][START_REF] Boyd | Convex optimization[END_REF], Savitzky-Golay filters [START_REF] Candan | A unified framework for derivation and implementation of Savitzky-Golay filters[END_REF], Stein's unbiased risk estimate (SURE) regularization [START_REF] Eldar | Generalized SURE for Exponential Families: Applications to Regularization[END_REF], blind or semi-blind source separation (for multi-channel signals) [START_REF] Comon | Handbook of Blind Source Separation: Independent Component Analysis and Applications[END_REF], or even empirical techniques [START_REF] Boudraa | EMD-Based Signal Filtering, Instrumentation and Measurement[END_REF] have been used, which are merely based on prior assumptions such as signal smoothness, time or transform-domain sparsity, or stochastic independence. Specifically, Tikhonov regularization has found great interest in various fields of science and engineering (including biomedical engineering) for smoothing noisy data and signal detrending [START_REF] Tarvainen | An advanced detrending method with application to HRV analysis[END_REF][START_REF] Boyd | Convex optimization[END_REF][START_REF] Zhang | Bioelectric signal detrending using smoothness prior approach[END_REF]. Despite its simplicity and good performance over short stationary data blocks, Tikhonov regularization becomes computationally demanding for long data. Moreover, it is an offline procedure and does not adapt to signal and noise nonstationarities.

Herein, a generic filtering scheme is proposed based on smoothness priors and a constrained least squares problem, which extends Tikhonov regularization to a block-wise formulation that is computationally efficient and applicable to nonstationary scenarios. Since the method is applied over rather short segments of data, it is shown that it is computationally more efficient than the classical Tikhonov regularization and it can be implemented as an online fixed-lag smoother with smooth inter-segment transitions, for both stationary and nonstationary mixtures of signal and noise. The smoother is also studied in the frequency domain and is shown to be equivalent to an adaptive filter, which tunes its bandwidth according to the signal/noise contents of each data segment. It is shown that the proposed method bridges between various concepts and filtering schemes, including smoothing splines, Lipschitz regularization [START_REF] Mallat | A Wavelet Tour of Signal Processing: The Sparse Way[END_REF], Wiener filtering, and closely related to SURE regularization [START_REF] Eldar | Generalized SURE for Exponential Families: Applications to Regularization[END_REF], and Savitzky-Golay filters [START_REF] Candan | A unified framework for derivation and implementation of Savitzky-Golay filters[END_REF].

The denoising of Electrocardiogram (ECG) signals is used as a case study for the proposed method. The problem of ECG denoising is of major importance in biomedical signal processing. Due to the large diversity of ECG waveforms, the current gold standard for ECG acquisition system front-ends are based on rather simple bandpass filters [START_REF] Kligfield | Recommendations for the standardization and interpretation of the electrocardiogram: part I: the electrocardiogram and its technology[END_REF]. More effective techniques for ECG denoising, such as wavelet denoising, are based on the concept of wavelet shrinkage using signal smoothness priors [START_REF] Kabir | Denoising of ECG signals based on noise reduction algorithms in EMD and wavelet domains[END_REF]. Considering the broad range of parameters, which can be used for wavelet denoising, various studies have been carried-out to study the best set of wavelet denoising parameters for ECG signals [START_REF] Singh | Optimal selection of wavelet basis function applied to ECG signal denoising[END_REF]. More recently, Savitzky-Golay filters using the SURE for optimal order selection have also been used for ECG denoising and compared with various waveletbased schemes [START_REF] Krishnan | On the selection of optimum savitzky-golay filters[END_REF][START_REF] Sadasivan | An optimum shrinkage estimator based on minimum-probabilityof-error criterion and application to signal denoising[END_REF].

More advanced (yet sophisticated) methods for ECG denoising have used the pseudo-periodicity of the ECG for constructing model-based filters using the extended Kalman filter and its extensions [START_REF] Sameni | Filtering Electrocardiogram Signals Using the Extended Kalman Filter[END_REF][START_REF] Sameni | A nonlinear bayesian filtering framework for ECG denoising[END_REF][START_REF] Sameni | Model-based Bayesian filtering of cardiac contaminants from biomedical recordings[END_REF]. The idea has been further utilized for ECG segmentation [START_REF] Sayadi | ECG denoising and compression using a modified extended Kalman filter structure[END_REF], decomposition [START_REF] Kheirati Roonizi | A signal decomposition modelbased bayesian framework for ecg components separation[END_REF], anomaly detection [START_REF] Sayadi | Robust detection of premature ventricular contractions using a wave-based bayesian framework[END_REF], fiducial point extraction [START_REF] Akhbari | Ecg denoising and fiducial point extraction using an extended kalman filtering framework with linear and nonlinear phase observations[END_REF], and fetal ECG extraction [START_REF] Sameni | Extraction of Fetal Cardiac Signals from an Array of Maternal Abdominal Recordings[END_REF][START_REF] Niknazar | Fetal ECG extraction by extended state Kalman filtering based on single-channel recordings[END_REF]. Despite the effectiveness of these highly modelbased filtering schemes (and their outperformance versus wavelet denoising [START_REF] Sameni | A nonlinear bayesian filtering framework for ECG denoising[END_REF]) they are mainly applicable to normal ECG with pseudo-periodic behavior and their performance degrades in cardiac anomalies with highly abnormal beat morphologies. However, in the hereby studied case studies, it is shown that the proposed method, which is merely based on smoothness priors (similar to wavelet denoising), is equally applicable to normal and abnormal ECG contaminated by broad-band noise, without requiring the ECG R-peaks. Due to its minimal assumptions on the ECG and noise structure, it can be potentially used in a broad range of ECG-based applications, even in very low signal-tonoise ratios (SNR).

The organization of the paper is as follows. In Section 2, a piecewise smooth matrix form of the proposed algorithm together with the parameter selection details are presented. The proposed scheme is related to the concept of Lipschitz regularity in Section 3. In Section 4 using a sample-wise forward-backward smoothing formulation, a frequency domain interpretation is presented, which gives a Wiener smoothing perspective to Tikhonov regularization and the proposed method. A generalization of the developed framework to multipleorder smoothness is presented in Section 5; followed by a detailed case study for ECG denoising, in Section 6. The paper ends with a summary and promising future directions, in Section 7.

Piecewise smooth Tikhonov regularization

Discrete smoothing operator

Consider d d , a length M finite impulse response approximation for the dth order derivative operator. As a simple realization of d d , one can use the following recursion, which is based on the first order difference approximation of the derivative operator:

d d = d d-1 * d 1 , d ≥ 2 (1)
where

d 1 ∆ = (1, -1)
is the first-order difference impulse response and * is the convolution operator. Throughout the paper, the second-order difference impulse response

d 2 = (1, -2, 1) is of special interest.
For integers T > M, Dd ∈ R (T -M+1)×T is defined as the Toeplitz matrix form of d d . As a special case we have

D2 =                     1 -2 1 0 • • • 0 0 1 -2 1 . . . . . . . . . . . . . . . . . . . . . 0 0 • • • 0 1 -2 1                     (2)
The left multiplication of Dd in any column-wise signal of length T is equivalent to the convolution of the signal with d d , which is its dth order approximate derivative in our case. The effectiveness of finite difference approximations of the derivative depends on the sampling frequency of the data and the desired accuracy. In practice, (2) can be replaced by discrete-time differentiators with a desired impulse response or even a weighted difference (cf. [START_REF] Candan | A unified framework for derivation and implementation of Savitzky-Golay filters[END_REF] for a Savitzky-Golay perspective).

Data model

Let us consider a signal θ(t) with a decaying Fourier transform Θ(ω), i.e.,

+∞ -∞ |Θ(ω)|(1 + |ω| ν )dω < +∞ (3)
where ν > 0 determines the rate of decay. Equation (3) implies that θ(t) is globally d = ν times continuously differentiable (cf. [START_REF] Mallat | A Wavelet Tour of Signal Processing: The Sparse Way[END_REF]Theorem 2.5]). Further consider θ[m] = θ(mT s ) as discrete-time samples of θ(t), x[m] (0 ≤ m < L) as noisy samples of θ[m] contaminated by uncorrelated additive noise, and a strictly monotone sequence of user-defined (or automatically selected) knot points

K = {p 0 , p 1 , • • • , p N }, p k ∈ [0, L).
The knots are used to segmentize the noisy signal into N non-overlapping segments. Accordingly, the segment S k (0 ≤ k < N) is defined in vector form as

x k ∆ = (x[p k ], x[p k + 1], • • • , x[p k+1 -1]) T (4)
with length n k ∆ = p k+1p k . We fix p 0 = 0 and p N = L to assure that each sample of x[m] falls in exactly one segment.

In each segment S k , the following signal model is considered

x k = θ k + η k (5)
where θ k ∈ R n k is the noiseless (desired) signal, and the noise vector η k ∈ R n k is assumed to be uncorrelated with θ k and zero-mean with variance η T k η k /n k = δ 2 k (a segment dependent variance).

In order to consider the inter-segment boundary behavior of the proposed method, two vectors containing q neighboring samples of S k from segments S k-1 and S k+1 are defined:

α k ∆ = (θ[p k -q], θ[p k -q + 1], • • • , θ[p k -1]) T ∈ R q β k ∆ = (θ[p k+1 ], θ[p k+1 + 1], • • • , θ[p k+1 + q -1]) T ∈ R q (6
) Exceptionally, for S 0 and S N-1 , which do not have any left and right neighbors, respectively, α 0 and β N-1 are defined to be empty vectors.

The knot-sequence and segment boundaries are illustrated in Fig. 1. Accordingly, the augmented vector

θT k ∆ = [ α T k θ T k β T k ] (7)
contains (n k + 2q) successive samples of the desired signal from the end of S k-1 to the beginning of S k+1 . In the following, we set q = M -1 to guarantee smooth transitions between successive segments.

q q α k β k x k θ k Sk-1 Sk Sk+1 pk pk+1 S k S k+1 p k+1
Figure 1: Illustration of the data model for each segment and its boundaries Using these definitions, the objective is to find the least square error (LSE) estimate of θ k with controllable degree of smoothness over S k and its boundaries with neighboring segments. The problem can be formulated as a LSE problem with an inequality constrain:

θk = arg min θ k x k -θ k 2 , s.t. Dd θk 2 ≤ 2 k ( 8
)
where Dd ∈ R (n k +M-1)×(n k +2M-2) is the order d difference operator, • represents the L2-norm and 2 k is a predefined upper bound for the signal's "steepness". The Lagrangian form of ( 8) is as follows:

θk = arg min θ k x k -θ k 2 + λ k Dd θk 2 (9)
where λ k is a positive regularization factor that is tuned to minimize the cost, while satisfying the inequality constraint in [START_REF] Boudraa | EMD-Based Signal Filtering, Instrumentation and Measurement[END_REF]. Equation ( 9) is in the form of Tikhonov regularization problems [START_REF] Golub | Tikhonov regularization and total least squares[END_REF][START_REF] Tarvainen | An advanced detrending method with application to HRV analysis[END_REF], and similar to the regularized SURE cost function [START_REF] Eldar | Generalized SURE for Exponential Families: Applications to Regularization[END_REF]. As noted in [START_REF] Tarvainen | An advanced detrending method with application to HRV analysis[END_REF], a nice interpretation for the constrained LSE is to draw the unconstrained LSE solution toward the null space of the matrix Dd , resulting in a smooth estimate of θ k .

If Dd is partitioned into Dd 1) and D d ∈ R (n k +M-1)×n k , the argument of (9) can be expanded as

= [ E d D d F d ], where E d , F d ∈ R (n k +M-1)×(M-
J k = x k -θ k 2 + λ k D d θ k + b k 2 ( 10
)
where b k [START_REF] Zhang | Bioelectric signal detrending using smoothness prior approach[END_REF] is in the form of a classical LSE problem with a quadratic constraint. It is straightforward to show that the solution is as follows

∆ = E d α k + F d β k is a boundary vector indepen- dent of θ k . Equation
θk = (I + λ k D T d D d ) -1 [x k -λ k D T d b k ] (11)
In the special case that the boundary smoothness constraint is not required (the boundary conditions are not available or the data has only one segment), b k = 0 and (11) reduces to the Tikhonov regularization method used in previous studies [START_REF] Tarvainen | An advanced detrending method with application to HRV analysis[END_REF]:

θk = (I + λ k D T d D d ) -1 x k (12)
We refer to (I

+ λ k D T d D d ) -1
, as the smoothing matrix. This matrix is shown in Fig. 2 for T = 100 and λ = 50, 200.

In Appendix A, it is shown that an exactly similar solution is obtained by seeking "maximally smooth estimates under an upper bound on the noise variance of each segment," using a segment-wise regularization factor γ k = 1/λ k .

For the hereby proposed method (and the alternative formulation stated in Appendix A), there are practical considerations regarding the selection of the regularization factors λ k (or γ k ), and the boundary vectors α k and β k , which are discussed below.

Regularization factor selection

Depending on whether the upper bound of signal roughness k is know or not, the regularization factor λ k can be fixed or calculated per-segment. Both cases are studies below.

Known (presumed) smoothness bound

When the smoothness bound is presumed or known (e.g., from prior knowledge of the signal form), an optimal λ k can be calculated using k . For this, the singular value decomposition (SVD) of D d is performed:

D d = UΣV T ( 13
)
where U ∈ R (n k +M-1)×(n k +M-1) and V ∈ R n k ×n k are orthonormal matrices and

Σ ∈ R (n k +M-1)×n k is a diagonal matrix, with diagonal entries (σ 1 , σ 2 , • • • , σ n k ). Next, letting θk ∆ = V T θk , xk ∆ = V T x k , bk ∆ = U T b k , (14)
it is straightforward to show that the block solution in [START_REF] Mallat | A Wavelet Tour of Signal Processing: The Sparse Way[END_REF] reduces to a sample-wise update:

θik = xik -λ k σ i bik 1 + λ k σ 2 i , (i = 1, • • • , n k ) ( 15
)
where θik , xik , and bik are the ith elements of the n kvectors θk , xk , and bk , respectively. Moreover, the mean square error and the smoothness bound are 16)

e(λ k ) = x k -θk 2 = xk -θk 2 = n k i=1 ( σ 2 i xik + σ i bik 1/λ k + σ 2 i ) 2 (
c(λ k ) = D d θk + b k 2 = Σ θk + bk 2 = n k i=1 ( σ i xik + bik 1 + λ k σ 2 i ) 2
(17) It is seen that e(λ k ) and c(λ k ) are monotonically increasing and decreasing with λ k , respectively. Therefore, the minimum of ( 8) is obtained on the boundary of its constraint region.

When the smoothness bound k is known (or assumed), the optimal regularization factor, denoted by λ * k , is numerically calculated by solving the following equation for λ * k :

c(λ * k ) = n k i=1 ( σ i xik + bik 1 + λ * k σ 2 i ) 2 = 2 k ( 18
)
Due to the positive and monotonic decay of (18) toward zero by increase of λ * k , it has a unique solution if 2 k ≤ c(0). In order to find λ * k , we use a variant of the Bisection algorithm [START_REF] Press | Numerical Recipes in C: The Art of Scientific Computing[END_REF], detailed in Appendix B. This optimal value is next used in [START_REF] Mallat | A Wavelet Tour of Signal Processing: The Sparse Way[END_REF] or [START_REF] Krishnan | On the selection of optimum savitzky-golay filters[END_REF] to calculate θ k . For 2 k > c(0), the boundary of the constrain region is no longer achieved and λ * k = 0, resulting in θk = x k .

Unknown smoothness bound

When k is unknown, or no assumption is available on its upper bound (or equivalently on the noise variance δ 2 k , as shown in Appendix A), methods such as the L-curve have been used for finding optimal regularization factors [START_REF] Gander | Least squares with a quadratic constraint[END_REF][START_REF] Golub | Matrix Computations, 3rd Edition[END_REF][START_REF] Hansen | The L-curve and its use in the numerical treatment of inverse problems[END_REF][START_REF] Calvetti | Tikhonov regularization and the l-curve for large discrete ill-posed problems[END_REF][START_REF] Boyd | Convex optimization[END_REF]. In [START_REF] Oraintara | A method for choosing the regularization parameter in generalized tikhonov regularized linear inverse problems[END_REF], a maximum a posteriori (MAP) interpretation has been presented for the constrained least squares problem and its optimal regularization factor selection using the Lcurve. In our case, the L-curve is the plot of x kθ k versus D d θ k + b k , as λ is swept from zero to infinity. The λ corresponding to the corner of the L-curve (the point at which x kθ k + D d θ k + b k is minimized) can be empirically chosen as the optimal regularization factor, which provides the best balance between the signal smoothness and minimum mean square error. The results of this method versus the aforementioned method, which is based on prior knowledge of the signal noise variance are detailed in the results section.

Depending on the problem, one may have additional priors regarding the signal/noise structure, which can improve the regularization factor estimate. For instance, if the signal and noise are both nonstationary from one block to another (like the ECG example studied in Section 6.6), but the noise is known to have rather slow inter-block variations, the regularization factor λ * k can be made more consistent and its inter-segment fluctuations can be damped by using the following block-wise updates:

λ * k = κλ * k + (1 -κ) λ * k-1 (19)
where 0 < κ ≤ 1 is a forgetting factor, which adds a degree of freedom for block-wise smoothing of λ * k . This approach is further justified in Section 4, where we present a more physical interpretation for λ * k (beyond a simple regularization factor).

Knot sequence selection

In theory, the knot sequence K is basically arbitrary. In practice, in order to minimize the impact of smoothing at the segment boundaries, one can customize the knot sequence on a signal-dependent basis. For instance, for ECG signals, one can detect the fiducial points of the ECG and set the knots at intermediate points, which do not convey major clinical information. Alternatively, for online filtering, the knots can be equally spaced throughout the signal. Our experiments over biosignals such as the ECG, show that by introducing the segment smoothness margins α k and β k in [START_REF] Mallat | A Wavelet Tour of Signal Processing: The Sparse Way[END_REF], the boundary effects are negligible, and the fixed spaced knots is more appropriate in practice (imposing no processing overhead).

Boundary conditions

Recalling [START_REF] Mallat | A Wavelet Tour of Signal Processing: The Sparse Way[END_REF], it is seen that for estimating the signal in each segment, the noiseless samples (or an estimate) of the boundaries of its left (α k ) and right (β k ) segments have been assumed to be known; while in practice they are not known. For this, we use a two-run procedure for estimating the boundaries.

Starting from a given knot sequence

K = {p 0 , p 1 , • • • , p N } (p k ∈ [0, L
], p 0 = 0 and p N = L), in the first run, a boundary-less solution is found for each block according to [START_REF] Kligfield | Recommendations for the standardization and interpretation of the electrocardiogram: part I: the electrocardiogram and its technology[END_REF]. Due to the local behavior of the smoothing matrix shown in Fig. 2, the resulting estimates are smooth in their intermediate points; but not necessarily at the segment boundaries.

In the second run, a secondary knot sequence K

= { p0 , p1 , • • • , pN } is defined: p0 = p 0 pi = (p i + p i-1 )/2, 0 < i < N pN = L (20)
which is a 50% shifted and overlapped version of K (except for the first and last points, which are fixed at 0 and L, respectively), and the signal is re-estimated from [START_REF] Mallat | A Wavelet Tour of Signal Processing: The Sparse Way[END_REF], using the boundary vectors α k and β k , from the first run.

The block-wise algorithm

Finally, by putting together the details in this section, the overall block-wise algorithm is summarized in Algorithm 1.

Algorithm 1 Block-wise smoothing algorithm

Require: Noisy signal x[m] (0 ≤ m < L) Require: Knot sequence K = {p 0 , p 1 , • • • , p N } Require: Upper bound of segment steepness 2 k (or up- per bound of segment noise variance δ 2 k ) 1: Segmentize x[m] into N segments x k using K 2: for k = 0, • • • , N -1 do First Run 3: Calculate λ * k for x k , using 2 k in Algorithm 2 4:
Estimate θk from [START_REF] Kligfield | Recommendations for the standardization and interpretation of the electrocardiogram: part I: the electrocardiogram and its technology[END_REF], without boundary overlaps 5: end for Re-estimate θk from [START_REF] Mallat | A Wavelet Tour of Signal Processing: The Sparse Way[END_REF], using the boundary vectors 12: end for It should be noted that this algorithm may be simply implemented in real-time with one segment of lag. In this case, the first and second runs of the algorithm can be interlaced over successive blocks.

6: Calculate K = { p0 , p1 , • • • , pN } from (20

Local Lipschitz regularity

There is a close relationship between the proposed block-wise scheme and the concept of global versus local regularity. In fact, under the Fourier decaying assumption in (3), θ(t) is uniformly Lipschitz-ν, over R [START_REF] Mallat | A Wavelet Tour of Signal Processing: The Sparse Way[END_REF]Ch. 6], which means that there exists a polynomial φ u (t) of degree d = ν and a constant > 0 (independent of ν), for which

∀t ∈ R : |θ(t) -φ u (t)| ≤ |t -u| ν (21)
The parameter ν ≥ 0 is the Lipschitz (or Hölder) exponent and the smallest , which satisfies [START_REF] Kheirati Roonizi | A signal decomposition modelbased bayesian framework for ecg components separation[END_REF]

[m] = θ(mT s ) (m ∈ Z).
On the other hand, by selecting a "sufficiently small" sampling time T s > 0, the sequence obtained from applying the approximate difference operator d d to the sample sequence θ[m] has the same order as θ(t) (d) and is bounded by the same upper bound in [START_REF] Kheirati Roonizi | A signal decomposition modelbased bayesian framework for ecg components separation[END_REF].

While ν and are global constants, if we confine the study to a single segment S k , the upper bound in [START_REF] Kheirati Roonizi | A signal decomposition modelbased bayesian framework for ecg components separation[END_REF] can become tighter and segment dependent. In other words, one can have local regularity Lipschitz exponents ν k and Lipschitz constants k , which vary from one segment to another. The block-wise scheme proposed in the previous section, implicitly seeks the upper bound 2 k through the segment dependent regularity factor λ k . This is the same idea behind the zooming property of wavelet transforms [START_REF] Mallat | A Wavelet Tour of Signal Processing: The Sparse Way[END_REF]Ch. 6].

Frequency domain analysis: a block-wise Tikhonov-Wiener smoother

The block-wise matrix-form scheme developed in the previous section can also be stated in terms of a zero-phase, forward-backward linear time-invariant (LTI) smoother. Although, the transfer function of classical Tikhonov regularization has been already derived in previous studies [START_REF] Kirkeby | Fast deconvolution of multichannel systems using regularization[END_REF][START_REF] Dumuid | A comparison of filter design structures for multi-channel acoustic communication systems[END_REF], our formulation additionally provides a frequency domain interpretation for the classical Tikhonov regularization in terms of a Wiener smoother. To show this, the sample-wise form of the data model in ( 5) is used. In each data segment k, we have

x k [m] = θ k [m] + η k [m], (0 ≤ m < n k ) ( 22
)
where x k [m] is the noisy signal, θ k [m] is the desired signal and η k [m] is the zero-mean observation noise, with

E{η k [m] 2 } = δ 2 k .
Using the difference operator defined in Section 2.1, the dth order approximate derivative of θ k [m] is found as follows:

w k [m] ∆ = d d [m] * θ k [m] = d d [0]θ k [m] + • • • d d [M -1]θ k [m -M + 1]
(23) Therefore, (9) can be rewritten as

θk [n] = arg min θ k [n] n k -1 m=0 (x k [m] -θ k [m]) 2 + λ k n k +M-1 m=M w k [m] 2 (24)
where λ k ≥ 0 is the regularization factor 1 . By setting the derivative of ( 24) with respect to θ k [n] equal to zero and after some simplifications, we find:

(δ[n] + λ k d d [n] * d d [-n]) * θk [n] = x k [n] ( 25
)
where δ[•] is the Kronecker delta. If the first term in ( 25) is invertible, we arrive at

θk [n] = (δ[n] + λ k d d [n] * d d [-n]) -1 * x k [n] (26)
which simply represents LTI non-causal smoothing of the noisy signal using the following transfer function:

ψ k (z) ∆ = 1 1 + λ k d d (z)d d (z -1 ) (27)
Equation ( 27) can be regarded as the transfer function of Tikhonov regularization (smoothing).

It is insightful to rewrite [START_REF] Gander | Least squares with a quadratic constraint[END_REF] in the following form

ψ k (z) ∆ = 2 k d d (z)d d (z -1 ) 2 k d d (z)d d (z -1 ) + δ 2 k = S θ (z) S θ (z) + S η (z) (28)
which is exactly in the form of an optimal Wiener smoother for a random process obtained from an autoregressive all-pole model with a spectra denoted by S θ (z), driven by input white (innovation) noise with variance 2 k , in presence of white observation noise with variance δ 2 k and a white spectra denoted by S η (z). In this case, the parameter λ k = δ 2 k / 2 k is the equivalent regularization factor defined throughout this work.

Fig. 3, illustrates the analogy between Tikhonov regularization and optimal Wiener smoothing using the sample-wise formulation. This shows that the method proposed in the previous section can be envisaged as a piecewise-stationary Wiener smoother with smooth inter-block transitions. The frequency response of each block is controlled by the regularization factor λ k . Using a real-valued impulse response for d d , ψ k (z) is all-pole, lowpass and non-causal with conjugate symmetric pole pairs. Additionally, due to its invariance with respect to z and z -1 , its poles are reciprocal with respect to the unit circle. The only exception is the possible poles lying on the unit circle (which is studied below). In other cases, ψ k (z) can always be factorized as the cascade of a minimum-phase system ψ min k (z) (using the poles inside the unit circle) applied as a causal filter, and a maximum-phase system ψ max k (z) (using the poles outside the unit circle) applied as an anti-causal (reverse) filter. Therefore, a forward-backward smoothing scheme can be used to implement the smoother. Due to the reciprocity of the poles of ψ k (z) inside and outside the unit circle, only ψ min k is required for the filter implementation. The procedure is as follows: first find . The whole procedure is rather similar to the filtfilt function in Matlab and Octave, using ψ min k (z) as the transfer function. It should just be noted that according to the two-run scheme described in Section 2.6, during the forward and backward filtering phases, the output initial conditions of the convolutions in each segment k are taken from the outputs of the neighboring segments S k-1 and S k+1 , obtained from the first run. circle, the reciprocity of the poles is only attained when the number of poles is a multiple of four, i.e., M is odd. For even values of M, (at least) two of the conjugate symmetric poles are forced to be on the unit circle, and the remaining poles are reciprocal inside and outside the unit circle. Therefore, the factorization of ψ k (z) in terms of the cascade of a minimum-phase and a maximum-phase filter (without any unit-circle roots) is only achievable for odd values of M, resulting in a stable forward-backward smoother. The system is unstable for even values of M. As a special case, using the frequency form of the difference operators defined in (1), the frequency response of the forward-backward smoother is found as follows

v k [n] = ψ min k [n] * x k [n]. Next, v k [n] is time reversed to ṽk [n] = v[-n],
ψ k (e jω ) = 1 1 + λ k (2 -2 cos ω) d = 1 1 + λ k (2 sin ω 2 ) 2d (29
) which resembles the amplitude response of a discretized Butterworth filter of order d, in low frequencies (|ω| = |2π f / f s | 1). The filter response is monotonically decreasing over ω ∈ [0, π] and is maximally flat in its passband, in the sense that its first 2d -1 derivatives are zero at ω = 0 without any pass-band or stop-band ripples. Therefore, the formulation is also attractive from the filter design perspective. The root loci and the corresponding filter frequency response of ψ k (z) are shown in Figs. 4 and5 for different values of d and λ k . Accordingly, the block-wise scheme presented in the previous section can be interpreted as a smoother, which adapts its bandwidth according to the noise variance and the signal steepness controlled by λ * k .

Regularization factor selection

The properties of the LTI model under the variation of the parameter λ k is similar to the block-wise scheme: the solution becomes smoother as λ k increases from zero to infinity. From [START_REF] Hansen | The L-curve and its use in the numerical treatment of inverse problems[END_REF], it is noticed how λ k is inversely related to the equivalent lowpass filter bandwidth. If the upper-bound of the observation noise variance δ 2 k is known (or assumed), λ k has an optimal value, which can be numerically found by solving:

1 2π 2π |X(e jω )(1 -ψ k (e jω ))| 2 = Lδ 2 k ( 30
)
which is equivalent to solving e (γ * k ) = Lδ 2 k for γ * k , in the alternative formulation stated in Appendix A and then setting λ * k = 1/γ * k .

Generalization to multiple order smoothness

An interesting generalization of the proposed method is to penalize the norms of multiple derivatives of the estimated signal at the same time, but with different regularization factors. As compared to the formulation proposed in Section 2, this approach provides additional control over the estimated waveforms, as it draws the results towards a smaller subset of the null space of D d . This feature can be useful for extracting smooth signals of order d, with bounded lower order derivatives. In this case, the cost function in [START_REF] Zhang | Bioelectric signal detrending using smoothness prior approach[END_REF] can be rewritten as follows:

J M k = x k -θ k 2 + d i=0 λ (i) k D i θ k + b (i) k 2 ( 31
)
where b (i)

k ∆ = E i α k + F i β k
is the multiple order extension of the boundary vector and λ (i) k is the ith derivative regularization factor. It is straightforward to show that in this case, the block-wise solution in [START_REF] Mallat | A Wavelet Tour of Signal Processing: The Sparse Way[END_REF] is generalized as follows θk = (I

+ d i=0 λ (i) k D T i D i ) -1 [x k - d i=0 λ (i) k D T i b (i) k ] (32)
Similarly, the sample-wise solution in [START_REF] Niknazar | Fetal ECG extraction by extended state Kalman filtering based on single-channel recordings[END_REF], its corresponding transfer function in [START_REF] Gander | Least squares with a quadratic constraint[END_REF], and the frequency response in ( 29) become

(δ[n] + d i=0 λ (i) k d i [n] * d i [-n]) * θk [n] = x k [n] (33) ψ k (z) = 1 1 + d i=0 λ (i) k d i (z)d i (z -1 ) (34)
ψ k (e jω ) = 1 1 + d i=0 λ (i) k (2 sin ω 2 ) 2i (35)
The discussion regarding the stability is as before (cf. Section 4.1). As compared to the single order model, it is apparently more difficult to find multiple optimal regularization factors λ (i) k . However, they can be simply fixed by intuition or be found by Monte Carlo simulation over desired datasets.

Case study: electrocardiogram denoising

Motivation

The electrocardiogram (ECG) is a good example of non-stationary signals. While the normal ECG has a pseudo-periodic behavior (due to the cardiac rhythm), abnormal and ectopic ECG beats are highly irregular.

In modern digital ECG acquisition systems, in order to enable better visual interpretations and to cover the vast inter-subject variations, the ECG is commonly oversampled far above the minimum required rate. For instance, although the major diagnostic information of the ECG is below 100 Hz (in adults), current ECG acquisition standards recommend a sampling frequency of 250-500 Hz for low and moderate quality ECG, and 1-10 kHz for high quality recordings [START_REF] Kligfield | Recommendations for the standardization and interpretation of the electrocardiogram: part I: the electrocardiogram and its technology[END_REF]. On the other hand, apart from the baseline wander, other ECG noises and artifacts can span across the entire Nyquist band. Therefore, in moderate and high-quality digitized ECG, the desired waveform is rather smooth (the adjacent samples are highly correlated), while the noise samples are rather uncorrelated. Based on these facts, in this section the proposed denoising scheme is used for removing wide-band noise from ECG recordings.

It should be noted that in previous research, filtering techniques such as extended Kalman filters based on the pseudo-periodicity of ECG signals have resulted in very good results for normal ECG [START_REF] Sameni | A nonlinear bayesian filtering framework for ECG denoising[END_REF]. However, the performance of these methods significantly degrade for abnormal ECG, which makes them inapplicable for denoising highly irregular and ectopic beats. Therefore, LTI filters with fixed impulse responses are still considered as the gold standard of ECG acquisition system front-ends.

Database

Three standard databases from the PhysioBank are used for evaluation:

1. The MIT-BIH normal sinus rhythm database, which contains 18 long-term two-channel normal ECG recordings, sampled at 128 Hz [START_REF]The MIT-BIH Normal Sinus Rhythm Database[END_REF]. 2. The MIT-BIH arrhythmia database, which contains 48 half-hour two-channel ambulatory ECG recordings with annotated anomalies, sampled at 360 Hz [START_REF]The MIT-BIH Arrhythmia Database[END_REF]. 3. The PTB diagnostic database (PTB), which contains 549 fifteen channel ECG recordings, sampled at 1 kHz [START_REF]The MIT-BIH PTB Diagnosis Database[END_REF]. Due to the huge size of this database, only the first ninety subjects are used in this study.

These data are selected to show the performance of the proposed methods in low, moderate and high sampling frequencies for normal and abnormal ECG.

Algorithm parameters

In the following case studies, the smoothness order is set to d = 2, d 2 = (1, -2, 1), and the smoothing operator D2 is as defined in [START_REF] Donoho | Adapting to unknown smoothness via wavelet shrinkage[END_REF]. The boundary vectors α k and β k have a length q = M -1 (=2 in this case). The knot spacing has been fixed to 100 ms, which is a fraction of a typical ECG beat length and the ECG is roughly smooth up to its second order difference. The accuracy threshold of Algorithm 2 is set to C = 10 -8 , and the forgetting factor in ( 19) is set to κ = 1, to show the block-wise variations of the regularization factor λ * k without any smoothing.

In the later presented results, several modes of the proposed algorithm are studied:

1. Single block version (Section 4), using a global optimal λ 0 found by a naive grid search. 2. A known input noise variance (δ 2 k ) 3. Underestimated input noise variance (0.8δ 2 k ) 4. Overestimated input noise variance (1.2δ 2 k ). The motivation for the two latter cases is to check the performance when the pre-assumption on (or estimation of) the input noise variance is inexact.

Benchmark algorithms

To date, numerous filters have been used for ECG denoising. In order to keep the comparisons fair and generic, we select benchmarks, which are not data dependent, do not require prior assumptions (such as pseudo-periodicity), the detection of the R-peaks or other ECG fiducial points, and have been reasonably successful for ECG denoising in previous work.

Linear time-invariant (LTI) lowpass filter

Although simple, linear-phase LTI filters are currently the standard choice for ECG system front-ends, as they impose little assumptions on the signals and preserve their phase contents. Therefore, the first benchmark is a naive linear-phase lowpass Butterworth filter with variable cutoff frequency. The -3dB cutoff frequency was swept linearly over 100 values ranging from 5 Hz to the Nyquist frequency (half the sampling frequency), in three Butterworth filter orders 2, 4, and 6. The designed filters were applied to the data using the filtfilt function in Matlab (to have zero-phase lag and provide similar comparison conditions with the proposed method). According to the simulation procedure described in the sequel, it was found that the filter order 2 slightly outperformed higher orders.

Wavelet denoiser (WD)

The well-known wavelet shrinkage scheme developed by Donoho et al. was originally developed using signal smoothness priors [START_REF] Donoho | Adapting to unknown smoothness via wavelet shrinkage[END_REF]. The optimal performance of this method using appropriate mother wavelets and different wavelet shrinkage schemes has been excessively studied in the past two decades. WD is offline and block-wise and highly benefits from the signal properties in various decomposition levels, which are much more sophisticated (and perhaps more effective) than the simple smoothness priors used in this work. Nevertheless, we use it for performance comparison with the proposed method. Due to the diversity of the literature on WD performance on ECG signals, an exhaustive Monte Carlo simulation was performed on the three databases using various combinations of WD parameters, listed in Table 1. Overall, 6960 different parameter combinations have been studied, which is broader than our previous WD combination set reported in [START_REF] Sameni | A nonlinear bayesian filtering framework for ECG denoising[END_REF].

Each method was tested in presence of input white noise in variable SNR ranging from -15 dB to +30 dB. The average and standard deviation of SNR improvement obtained after filtering was used as a measure of performance.

Among the different tested parameter combinations, the SURE shrinkage rule, with a single-level rescaling and soft thresholding strategy resulted in the highest performance over all three databases. The multi-level rescaling performance was rather close (but lower) than the single-level; but the no-rescaling method always had poor results, except in very high input SNR (24dB and above). Among the tested mother wavelets, Coiflets 2, 3, 4, 5, and Symlets 5, 6, 7, 8, Daubechies 2 commonly gave superior results, without a significant difference (considering the standard deviation of the results). The optimal decomposition level for WD is known to be sampling frequency dependent. Throughout the Monte Carlo simulations, the optimal level was found to be 7, 5, and 9, for the normal, arrhythmia, and PTB databases, respectively. In the results section, the WD results are reported in two cases:

1. The best WD decomposition level adapted per each input SNR (named WAVELET 0 ) 2. A fixed decomposition level optimized over all input SNR (named WAVELET).

The first case is only used as a corner case, as the exact input SNR is not known in practice and the optimal decomposition level is unknown.

Order-optimized Savitzky-Golay filter (SGF)

The Savitzky-Golay filter using the SURE shrinkage rule for optimal order selection has been recently used for ECG denoising [START_REF] Krishnan | On the selection of optimum savitzky-golay filters[END_REF]. Three variants of the SGF have been reported in [START_REF] Krishnan | On the selection of optimum savitzky-golay filters[END_REF]:

1. SGF of order 3 with variable length impulse response and using SURE for finding the optimal impulse response (SGBW). 2. SGF of order 3 with variable length impulse response and using regularized SURE for finding the optimal impulse response (SGBWR). 3. SGF of variable order from 1 to 5 using SURE for finding the optimal impulse response length (SGO).

All three schemes were applied to the studied databases 2 , with the order parameter swept from 5 to 75 in steps of 2 (as required by the algorithm in [START_REF] Krishnan | On the selection of optimum savitzky-golay filters[END_REF]). It was noticed that the optimal filter order was inversely related to the input SNR and the performance highly relied on this parameter. However, as with all the benchmark methods, the optimal order was selected as the one 0 0.2 0.4

Figure 6: A sample noisy ECG segment (in gray), the smoothed signals after the first stage of Algorithm 1 without any boundary constraints (in blue), and after the second stage using boundary constraints (in red).

that maximized the average SNR improvement over the entire range of input SNR.

6.5. Stationary white noise removal 6.5.1. Visual inspection In order to show the impact of considering the boundary conditions between segments, a short ECG segment adopted from the PTB diagnostic database and the results obtained from the first and second runs of Algorithm 1 are shown in Fig. 6. Accordingly, the transitions between blocks have become totally smooth after the second stage of the smoother (in Algorithm 1). Both runs give rather similar estimates in their intermediate points, which shows that the segment boundary constraints have local impacts and do not affect the interior samples.

In Fig. 7, a sample ECG segment adopted from the MIT-BIH arrhythmia database is shown, which contains ectopic beats. The signal was artificially contaminated by white Gaussian noise in three input SNR: 0dB, 6dB and 12dB. The denoising results using the block-wise and sample-wise schemes are shown in Fig. 7. In order to calculate the optimal λ * k , the alternative formulation of the algorithm described in Appendix A was used and the noise variance δ 2 k was given to the algorithm as a prior. For the sample-wise smoother (Section 4), a single regularization factor λ was fixed to be the maximum of λ * k obtained from the block-wise scheme. The optimal penalty factor γ * k = 1/λ * k is shown in the last row of Fig. 7. Accordingly, although the noise variance was uniform throughout the whole signal, γ * k has significant block-wise variations, due to the local ECG amplitude changes over each beat.

The L-curves corresponding to the data in Fig. 7 are shown in Fig. 8. Accordingly, for all the plotted L-curves, due to the fixed signal SNR, the optimal λ * k (shown by circles) are located at a fixed error level x kθ k , which by visual inspection are rather close to the corner of the L-curve. Therefore, for real scenarios in which δ 2 k (or 2 k ) are unknown, according to our primary model in [START_REF] Boudraa | EMD-Based Signal Filtering, Instrumentation and Measurement[END_REF] an upper bound of these parameters can be used instead. Choosing worst case upper bounds for δ 2 k (or 2 k ) leads in sub-optimal results. Alternatively, a sweep of the regularization factor for obtaining the corner of the L-curve gives an empirical estimate for λ * k .

Monte Carlo simulation

In order to compare the overall performance of the proposed method, it was applied to the three standard databases described in Section 6.2. The total number of 10 s data segments were 1116, 2016, and 1150, for the Normal, Arrhythmia, and PTB databases, respectively. White noise was synthetically produced and added to each segment in different SNR ranging from -15 dB to 30 dB in 3 dB steps, and the performance was studied for three difference operator orders d =2, 4, and 6. High SNR (above 24 dB) cases are considered very good quality ECG; intermediate SNR (ranging between 18 dB to -6 dB) have been used as a typical range for noise stress testing studies [START_REF] Moody | A noise stress test for arrhythmia detectors[END_REF], and very low input SNR (-6 dB and below) are only of interest for extremely noisy cases demanding rather coarse ECG analysis, such as R-peak detection from poor quality data (for instance for noninvasive fetal R-wave detection from highly noisy data [START_REF] Sameni | A Review of Fetal ECG Signal Pro-cessing; Issues and Promising Directions[END_REF]).

For this study, the proposed methods were tested using the parameters described in Section 6.3. For the block-wise scheme the noise variance was given as prior and the optimal regularization factor was calculated according to the γ-based formulation, detailed in Appendix A. For the sample-wise forward-backward smoothing scheme, an optimal global regularization factor γ 0 was found for each data segment by using a bisection search between 10 -12 to 10 6 over 100 iterations. The optimal regularization factor γ * 0 , was selected as the one which resulted in the maximum SNR improvement (∆SNR(dB) = SNR out (dB) -SNR in (dB)) at the filter's output. The average and standard deviation of SNR improvements are shown in Fig. 9, for the three studied databases, in comparison with the benchmark methods. Accordingly, in low input SNR, WAVELET 0 has the best average SNR improvement, which was anticipated due to the scenario used for implementing WAVELET 0 , as noted in Section 6.4.2. In high input SNR, the proposed method using fixed λ 0 and underestimated input noise variance (0.8∆ 2 k ) outperform other methods. For the WD, as expected, optimizing the decomposition level per input SNR (WAVELET 0 ) gives better results as compared with a global decomposition level (WAVELET). The simple lowpass filter and the SGF highly degrade in high input SNR. Overall, the SGF does not show any superior results (neither in terms of average nor standard deviation), except for the SGBWR applied to the PTB database in input SNR ranging from 0dB to 10dB.

For the proposed method, the difference operator order d = 2, outperformed higher order factors. Moreover, the SNR improvement obtained in higher sampling frequencies (specifically for the PTB database) was better than the results obtained in lower sampling frequencies. This result was already anticipated, since in higher sampling rates, the noise power, which spans over the total Nyquist band has less overlap with the desired signal. In all cases, the optimal regularization factor γ * 0 is inversely related to the SNR. Moreover, the results show that it is better to overestimate the input noise variance in low input SNR and to underestimate it for high input SNR. Using the exact input noise variance has intermediate results.

Considering the standard deviation (SD) results in Fig. 9, all variants of the proposed method are more robust (have smaller SD) than the previous techniques, including the WD. In fact, while the SD of the proposed methods remain bounded with the variation of the input SNR, the SD of other methods tend to increase in higher input SNR.

Non-stationary noise removal

The experiment described in the previous section was repeated for non-stationary white noise modulated by sinusoidal envelopes, in three input SNR: -3dB, 0dB and 3dB (Fig. 10). This noise resembles non-stationary noise due to electrode displacements and muscle artifacts, which are commonly wider in bandwidth as compared with the contaminated ECG; but overlap in the time and frequency domains. In this case, due to the variations of the noise variance, using a fixed prior noise variance δ 2 k is inappropriate as it leads to sub-optimal performance in each block. For this case, we used the L-curve method to find the optimal regularization parameter λ * k . From the last row of Fig. 10, it is noticed how the variations of γ * k = 1/λ * k inversely related to the envelope of the non-stationary noise, which is in full agreement with the interpretation of λ * k , in Section 4 (Fig. 3), as the ratio of the observation noise variance divided by the innovation noise variance, i.e., γ * k is inversely proportional to the observation noise variance.

Alternatively, for such scenarios (with non-stationary variable noise variance), a block can be added to the proposed algorithm, which estimates the block-wise noise variance using complementary priors. This extension is not covered in this study.

Summary and future work

In this paper, a block-wise extension of the wellknown Tikhonov regularization was proposed for denoising smooth signals contaminated by wide-band noise. The proposed method can be implemented in two forms: 1) as a block-wise fixed-lag smoother with smooth inter-block transitions applied in a matrix form, 2) as a fixed-interval smoother applied as a forwardbackward zero-phase filter. The filter response is maximally flat and monotonically decreasing, without any ripples. The method was also extended to smoothness of multiple derivative orders and its relationship with Lipschitz regularity and Wiener smoothing was also studied.

ECG denoising was studies as a typical case study and compared with conventional bandpass, WD and SGF filters. All Matlab source codes of the proposed method are online available in the open-source electrophysiological toolbox (OSET) [START_REF] Sameni | The Open-Source Electrophysiological Toolbox (OSET)[END_REF].

The proposed method is very general with potential applications in various fields of signal processing. Some of the promising extensions of the hereby proposed schemes include: • The data model described in [START_REF] Candan | A unified framework for derivation and implementation of Savitzky-Golay filters[END_REF] was based on the assumption that the observed data x k were simply noisy observations of θ k . For applications such as model identification and signal compression, θ k can be further considered to have a structure, such as θ k = Hs k , where s k are some latent variables (hidden features) and H is a mixture matrix or even some basis functions.

• The same framework can be extended to other priors, such as temporal or spectral priors. This extension is simply achieved by replacing the operator D d with other filter responses, e.g., the original Savitzky-Golay convolution coefficient set [START_REF] Savitzky | Smoothing and differentiation of data by simplified least squares procedures[END_REF], or recent methods proposed for the optimal order selection of these filters [START_REF] Krishnan | On the selection of optimum savitzky-golay filters[END_REF]. The filter operator can also be time-varying, such as signal-dependent weighted averages.

• The ECG case study shows that WD using SURE shrinkage rule has very good results for offline applications, while the proposed framework is more robust (in terms of SD) and is applicable online. In future studies, the combination of the two methods may result in a piecewise online WD schemes with smooth inter-block transitions. The concept of Lipschitz regularity studied in Section 3 can serve as a common ground for relating wavelet zooming and wavelet denoising schemes with the proposed framework.

• The Wiener equivalent formulation stated in [START_REF] Golub | Matrix Computations, 3rd Edition[END_REF], suggests the possible extension of the proposed method for colored spectral noise (non-constant S η (z)), having spectral overlap with the desired signals. In case of non-stationary colored noise, spectral factorization methods can be used to derive dynamical models for the proposed filters, which can also be used to construct a Kalman filter (smoother) implementation of the proposed schemes.

• The extension to multichannel spatially correlated data is another promising direction, which can provide a joint spatio-temporal smoothing technique, related to semi-blind source separation techniques using smoothness priors.

• In the current formulation, no assumptions were made on the stochastic distribution of the signal or noise; since the constrained least squares scheme does not rely on the signal distributions. However, if any such assumptions exist in a specific scenario (such as Gaussian assumption on the signal or noise), the formulation can be extended to a Maximum Likelihood or Maximum a Posteriori framework, which provide additional statistical interpretations for the minimized costs. λ * k = (λ max + λ min )/2; 23: end if

3

 3

Figure 2 :

 2 Figure 2: The smoothing matrix of Tikhonov regularization, for T = 100 and λ = 50 (left) and λ = 200 (right). This matrix has a local effect on the intermediate samples, when left-multiplied in a signal vector.

) 7 : 9 :

 79 Re-segmentize x[m] into N segments xk using K 8: for k = 0, • • • , N -1 do Second Run Calculate λ * k for xk , using Algorithm 2 10: Set the boundary vectors α k and β k from the output of the first run 11:

Figure 3 :

 3 Figure 3: An equivalent Wiener smoother scheme for Tikhonov regularization

  and given to the minimum-phase filter y k [n] = ψ min k [n] * ṽk [n]. The final result is θk [n] = y k [-n]. Unity DC gain is achieved by normalizing the results of the forward and backward filters by n ψ min k [n]

4. 1 .Figure 4 :

 14 Figure 4: Root locus of the forward-backward smoother defined in (27), for d = 2. Each line segment is the track of one of the four conjugate symmetric and reciprocal poles, as λ k varies from 0 to ∞ (in directions of arrows).

Figure 5 :

 5 Figure 5: The frequency response of the forward backward smoother defined in (27), (a) for d = 2, as λ k varies from 0 to ∞ (from top to bottom), (b) for λ k = 1 and d = 2, 4, 6, 8, 10, 12, 14 (from smoothest to steepest). The maximally flat pass-band property of the filter is seen here.

Figure 7 :

 7 Figure7: A sample ECG segment containing an ectopic beat in three input SNR: 0dB (left), 6dB (middle) and 12dB (right). From top to bottom, the signals correspond to: the original signal overlaid with the noisy signal, the sample-wise smoother scheme with a fixed λ, and the block-wise matrix form. The last signal is the optimal penalty factor γ * k = 1/λ * k corresponding to each 100 ms data segment.

Figure 8 :

 8 Figure8: The L-curves corresponding to the data in Fig.7as λ is swept from 0.04λ * k to 25λ * k , where λ * k is the optimal regularization factor located at the circles.

  dB) ∆ SNR STD (dB)

Figure 9 :

 9 Figure 9: The results of Monte Carlo simulations over three standard databases. The average SNR improvement (first row) and SNR improvement standard deviation (second row). The left column corresponds with the Normal Sinus database, the middle column shows the Arrhythmia database and the right column is the result of the PTB database.

Algorithm 2

 2 Calculation of λ * k using Bisection Require: Bisection stopping accuracy threshold C Require: Maximum number of bisection iterations I max 1: if c(0) < 2 k then while |λ maxλ min | > C & i < I max do 14: i = i + 1; 15:λ mid = (λ max + λ min )/2; 16:if [c(λ min ) -2 k ][c(λ mid ) -2 k ] ≤ 0 then 17:λ max = λ mid ;

  is known as the Lipschitz constant. Equation (3) implies that θ(t) is d-times continuously differentiable (globally) and the dth order derivative θ(t)(d) is uniformly Lipschitz-µ, where 0≤ µ = (νd) < 1. At each u the polynomial φ u (t)is uniquely defined. Specifically, since θ(t) is dtimes continuously differentiable in a neighborhood of u, then φ u (t) is the Taylor expansion of θ(t) around u [11, Ch. 6, P. 206]. This result holds for all t ∈ R, including uniformly sampled values θ

Table 1 :

 1 Wavelet denoising parameter combinations tested over the databases

	Parameter	Values
	Mother wavelet	Daubechies 1 (Haar), 2, 3,• • • , 10, 12, 16 Coiflets 1,• • • , 5
		Symlets 2,• • • , 8
		Biorthogonal 1.5, 2.6, 2.8, 5.5, 6.8
	Shrinkage rule	SURE, heuristic SURE, universal, minimax
	Thresholding strategy hard, soft
	Rescaling approach	no-scaling, single-level, multi-level
	Decomposition level	1,• • • , 10

The lower and upper bounds of the second summation in[START_REF] Sameni | Extraction of Fetal Cardiac Signals from an Array of Maternal Abdominal Recordings[END_REF] are due to the fact that in extending[START_REF] Tarvainen | An advanced detrending method with application to HRV analysis[END_REF] to the sample-wise form, the first and last M -1 samples of w k [m] correspond with the boundary vectors and only the middle n k correspond to θ k [n].

The source codes for this benchmark were kindly provided by the authors of[START_REF] Krishnan | On the selection of optimum savitzky-golay filters[END_REF].

Acknowledgment

The author would like to thank Dr. Sunder Ram Krishnan and Dr. Chandra Sekhar Seelamantula for providing the Matlab source codes of the order-optimized Savitzky-Golay filter [15], used in our benchmarks.

Figure 10: A sample ECG segment containing an ectopic beat contaminated by non-stationary noise resembling electrode displacements and muscle artifacts, in three input SNR: -3dB (left), 0dB (middle) and 3dB (right). From top to bottom, the signals correspond to: the original signal overlaid with the noisy signal, the fixed-lambda sample-wise smoother scheme, and the block-wise scheme with λ * k found by the L-curve method. The last row is the optimal penalty factor γ * k = 1/λ * k corresponding to each 100 ms data segment.

Appendix A. An alternative formulation: maximal smoothness with bounded noise variance

In practice, the knowledge of the segment measurement noise variance δ 2 k can be more accessible than the signal's steepness bound 2 k . In this case, an alternative formulation for ( 8) is to seek the smoothest estimate of θ k under an upper bound for the mean square error variance:

In this case, the Lagrangian and the matrix solutions are respectively as follows

3) is identical to [START_REF] Mallat | A Wavelet Tour of Signal Processing: The Sparse Way[END_REF], when the error penalty γ k = 1/λ k . Using the SVD matrix definitions from Section (2.3), and defining θ k ∆ = V T θ k , the samplewise update, the mean square error and the smoothness bound are found as follows

where it is again noticed that the optimal solution is on the boundary of the constraint region, if n k δ 2 k ≤ e (0).

In this case, the optimal error penalty γ * k is found by numerically solving the monotonically decreasing function

using the modified bisection scheme described in Algorithm 2. As before, there is a unique solution on the constraint boundary for n k δ 2 k ≤ e (0). For n k δ 2 k > e (0), we have γ * k = 0. The discussion on possible temporal smoothing of γ * k is similar to Section 2.3.

Appendix B. Modified bisection for optimal regularization factor selection

The modified bisection algorithm used for calculating the optimal regularization factor described in Section 2.3 is detailed in Algorithm 2.


C L Bottasso

L Trainelli

IMPLEMENTATION OF EFFECTIVE PROCEDURES FOR UNILATERAL CONTACT MODELING IN MULTIBODY DYNAMICS

The paper discusses the treatment of general contact/impact problems in multibody mechanical systems, with particular attention to the implementation of such procedures in commercial codes. We describe a unilateral contact model based on a flexible parameterization of the contacting surfaces, a geometrical minimum distance problem and a constitutive model of the interaction forces. The resulting numerical procedures are demonstrated by means of examples.

Introduction

The modeling of contact/impact processes represents a key topic in advanced multibody simulation, and is crucial for an effective use of this technology in many industrial problems. Apart from mechanisms with variable topology that represent the majority of applications, intermittent contact can also be caused by the presence of clearance at the joints, possibly due to wear or imperfections. The applicability of multi body dynamics to the virtual prototyping of an ever increasing range of problems in multiple areas of engineering makes the modeling of unilateral contacts an important asset of any software tool that aims at being general.

Several commercial multibody codes are available, and some of them are widely used by industry. While research codes are invaluable for testing new algorithms and procedures, they usually can not compete with industrial software for number of implemented elements and features, graphic capabilities, maturity of the user interface and integration with other CAD/CAM/CAE tools. Furthermore, most commercial codes are expandable and customi:z;able through general user-defined subroutines that allow to add new features.

However, to date the most widely used multibody packages offer only primitive contact/impact modeling capabilities. Therefore, it would seem worthwhile to explore the idea of implement-ing a contact/impact model in a commercial package, rather than in a research code. In fact such a tool would be applicable to a vast class of impact problems, either representing normal operative situations or non-regular working conditions which nevertheless have to be analyzed in the design process to ensure safety and/or efficiency.

A numerical procedure with unilateral contact modeling capabilities needs to address the following problems: a) give a mathematical representation of the geometry of the contacting parts; b) solve the minimal distance problem between the curves or surfaces where contact occurs, in order to evaluate the position of the candidate contact points; c) apply, at the contact points, the interaction forces defined by a model that is appropriate for the phenomenon being investigated.

Scope of this work is not the development of new methodologies for dealing with unilateral contacts, but the critical analysis of various possible alternatives proposed in the literature, with the specific goal of implementing such models in existing commercial software through user-definable features, i.e. without accessing the program source code. While a mathematical model should be primarily judged for its appropriate description of the phenomenon being investigated and for its performance, it seems that its ease of integration into existing procedures has not been previously considered. In our analysis, we will assume that user defined routines only permit the coding of non-linear holonomic constraints and of forces that are function of the system state.

The paper is organized in the following manner. First, we discuss the two main categories of unilateral contact schemes, the impulsive and the continuous force models, and we point out their strengths and deficiencies in view of what has been pointed out above. Next we consider the problem of description and parameterization of the contacting curves and surfaces and of the determination of the minimum distance between them. These are standard problems thoroughly covered in the literature, and are therefore only briefly sketched here. We then discuss the practical implementation of the selected unilateral contact model in a commercial package, pointing out a few general "tricks" that allow to bypass common limitations of the existing codes. Finally, we demonstrate the implemented procedures with the help of numerical simulations of multibody systems characterized by intermittent motion using two popular codes, ADAMS 1 and LMS DADS

Impulsive and Continuous Force Models

The various approaches to the modeling of unilateral contact conditions fall into two main categories. The first considers an impact as an impulsive phenomenon of null duration [4,5]. The configuration of the system is "frozen" during impact, and an appropriate model is used for relating the states of the system immediately before and immediately after the event.

There are two alternative flavors to this theory: Newton's method and Poisson's method [START_REF] Pfeiffer | Multi body Dynamics with Unilateral Contacts[END_REF]. The former relates the relative normal velocities of the contacting bodies through the use of an appropriate restitution coefficient. The latter divides the impact in two phases: an initial compression phase brings the normal relative velocity of the bodies to zero through the application of an impulse at the contact location; then, an expansion phase applies an impulse of opposite sign. The magnitude of the second impulse is related to the magnitude of the first one through a restitution coefficient.

A second approach models the contact/impact condition as a finite duration event, and tries to describe the time history of the resulting interaction forces for the whole duration of the phenomenon (6,8,1). This is achieved by introd~cing a suitable phenomenological model for the contact forces, usually expressed as functions of the inter-penetration, or "approach", between the contacting bodies. The computation of the approach at each time instant of the simulation is achieved solving the same set of equations that express the minimum distance problem when the bodies are not in contact. As for all contact models, we have even in this case a complementarity problem: either the sum of the relative distance and the approach is greater than zero, and in this case the contact forces vanish, or the same sum is null, and the relative distance is equal and opposite to the approach (inter-penetration) with non-vanishing interaction forces. These models are well established in the literature, and need not be discussed in greater detail here. Furthermore, they both can deal with contacts with friction and rolling conditions. Although the impulsive model has been used with success for multibody contact/impact simulations [START_REF] Pfeiffer | Multi body Dynamics with Unilateral Contacts[END_REF], it is quite clear that it seems most suited to the description of events of very short duration. Furthermore, its implementation requires the interruption of the numerical integration of the equations of motion when an impact is detected, followed by the solution of a linear complementarity problem that implies the manipulation of entire descriptive matrices of the whole multibody system prior to the restarting of the integration process. The resulting procedure would not be easily implemented in proprietary software packages through the standard user-accessible channels, which are typically routines that allow to add constraint conditions and external forces. In particular, the requirement to access entire system matrices (e.g. the mass matrix) can not usually be easily met. Furthermore, stopping and restarting the computation requires some form of "driving" procedure that can pose further implementation problems.

The continuous force model does not rely upon a short duration assumption. In addition, its implementation seems simpler than in the previous case. In fact, a contact force is applied if the bodies are in contact on the sole knowledge of the inter-penetration and its time rate. These quantities are already available through the solution of the minimum distance problem.

Another important difference between the two methods is related to the accurate .determination of the contact event. Indeed, contact can only be captured within the accuracy of the current time step. This means that in general some time refinement will be necessary in the proximity of the events. With a continuous force model, contact will trigger sharp gradients in the computed system state variables. Assuming that the target software package is equipped with a reasonably reliable and robust time adaption scheme, the sharp gradients will automatically activate the selection of smaller steps. This will ensure the automatic determination of the contact conditions within a desired tolerance without the need to modify the code time marching procedures. On the contrary, time adaption with an impulsive force model could not rely on the same mechanisms, and some form of modification of the adaption algorithm would be necessary. This however might be not possible without accessing the source code.

In view of these considerations, for this work we select a continuous force model based on a phenomenological law that relates the magnitude of the approach to the force of contact. In general, the interaction force will be given by the sum of elastic and dissipative components. The elastic term can usually be based on the classical Hertz model. For the determination of the second term, we follow the classical work of Hunt and Crossley (7], where the dissipative forces are proportional to the elastic ones and are energetically equivalent to the Newton restitution coefficient model. This has the advantage of determining the dissipative component through a quantity that can be easily measured experimentally.

Parameterization of Contact Curves and Surfaces

Following the methodology just outlined, we have seen that we need to provide a mathematical description of the geometry of the contacting parts. This will allow the solution of the minimum distance problem, that in turn will trigger the activation of the interaction forces between the contacting parts. Therefore, some kind of parameterization of the geometric entities is required. It is clear that the parameterization adopted should be as general and flexible as possible, in order to• allow the greatest range of possible applications. One parameterization that is broadly used for its generality and flexibility is based on Non-Uniform Rational B-Splines (NURBS) [START_REF] Farin | Curves and Surfaces for Computer Aided Geometric Design[END_REF][START_REF] Piegl | The NURBS Book[END_REF], often used for CAD solid modeling, manufacturing, graphics etc. NURBS can represent everything from a simple 2D line or circle to the most complex free-form surface.

Based on the NURBS representation, it is possible to evaluate all the geometric attributes of the entities, for example tangents, normals and point locations. These geometric attributes become the ingredients of the minimum distance problem that is solved at each time instant of the simulation. A NURBS-based parameterization of curves and surfaces can be easily coded in a library of routines. This library is then linked with the rest of the solver, and provides the needed geometric functionalities to the user defined subroutines implementing the minimal distance problem.

The Minimum Distance Problem

Letting q be the minimal distance between two entities, we say that they are in contact if q < 0. Hence, at any integration time step, the distance q has to be evaluated and a contact force F must to be applied as

F = H(q)f(a, a), ( 1
)
where a = -q is the approach between the entities and H(q) is a Heaviside step function, null for positive q. To decide on the application of the impact force, the minimal distance between the entities must be determined. Let us formalize the problem that has to be solved at any time step to find the point eligible for contact in both the 2D and 3D cases [START_REF] Pfeiffer | Multi body Dynamics with Unilateral Contacts[END_REF].

For the 2D case, the regular curves Ci c lR 2 (i = I, II) are defined through the mapping [START_REF] Farin | Curves and Surfaces for Computer Aided Geometric Design[END_REF] where ~i is the curvilinear coordinate of the curve and Ui E C 1 (Ai)• The points eligible for contact are then obtained by solving the non-linear system:

nr(6) • tu(6r) = 0, dr-u(~r,~II) • tii(~II) = 0, (3) (4)
where t and n are the vectors representing the tangent and normal to the curves, respectively. Eq. ( 3) forces the normal vector on C 1 to be perpendicular to the tangent vector on Cu, while eq. ( 4) forces the distance vector dr-II to be perpendicular to the tangent vector on Cu. The non-linear system is solved by two values 6 and ~II corresponding to the minimal distance points on the curves.

For the 3D case, the regular surfaces S, C R. 3 ( i = I, II) are defined through the mapping ui :A, H R. 3 , where ~i and ' T/i are curvilinear coordinates on the surfaces, and u, E C 1 (B;).

The points eligible for contact are obtained by solving the following non-linear system of equations (

) (6) (7) 5
t and s are two orthogonal R. 3 vectors lying on the tangent planes and n is the R. 3 vector normal to the surfaces. Eqs. ( 5) and ( 7) force the normal vector on 8 1 to be perpendicular to the tangent vectors on SII; eqs. ( 6) and (8) force the distance vector to be perpendicular to the tangent vectors on SIJ. The non-linear system is solved by the values ~1 , 'T/I, (u, 'T/II corresponding to the minimal distance points on the surfaces.

Once the minimal distance points have been determined, in both the 2D and 3D case the minimal distance value q is given by the following expression q = dr-II • nr. [START_REF] Piegl | The NURBS Book[END_REF] Practical Implementation in General Proprietary Multibody Software Some care needs to be exercised when implementing the minimum distance problem in a proprietary code using the typical functionalities provided by the user defined subroutines. In fact, in general the point of application of user forces on a given body will need to be specified a-priori. This is clearly in contrast with the problem here considered, since the point of application of the contact forces depends on the location of the contact point, which is a problem unknown and can not be specified a-priori.

Implementation Through User Defined Constraint Subroutines

Two main approaches seem to be possible. In the first method, one treats the equations defining the minimum distance problem as a set of non-linear holonomic constraints that are added to the set of constraint conditions of the multibody system. This approach is similar to the way one would attack the problem in case of implementation in a research code. However, the problem needs to be slightly reformulated in order to be practically used.

First, one needs to introduce massless dummy parts that are constrained to move on the curves or surfaces that describe the geometry of the contacting bodies. The dummy parts will coincide at each instant with the candidate contact points, and will therefore allow the application of the contact forces. Clearly, these additional parts must be of null mass in order not to modify the dynamics of the system. For example, one can constrain the dummy .parts to move on the contact curves or surfaces with axes gw and g2D belonging to the tangent plane and axis g 3 v aligned with the local normal. These constraint relations are once again holonotnic constraint conditions that can be implemented using appropriate user constraint routines. The minimum distance problem is then defined using the axis system associated with the dummy parts and the interaction forces are applied to the same parts, rather than to the contacting bodies.

A total number of 14 scalar constraint equations is written for the 2D case, that determine the 6 degrees of freedom of each dummy part in addition to the ~~ and ~II parametric coordinates. For the 3D case, we have 16 scalar constraints, corresponding to the 6 degrees of freedom of each dummy part and to the 6, 'f/J, ~II, rm parametric coordinates. More precisely, we have in the 2D case: • eq. ( 10) states that the g 2 v axis of a dummy part must be orthogonal to the curve tangent;

• eq. ( 11) states that the g 3 v axis of a dummy part must be orthogonal to the curve tangent;

• eq. ( 12) states that the g3D axis of a dummy part must be orthogonal to bi-normal b to the curve; ,

• eq. ( 13) states that the normal to curve I must be orthogonal to the tangent to curve II;

• eq. ( 14) states that the relative distance between curve I and I I must be orthogonal to the tangent to curve I I;

• eq. ( 15) translates the point on curve constraint for a dummy part and its associated curve.

Eqs. (10-12) must be repeated for both curve I and //.

For the 3D case we have the same eqs. [START_REF] Pfeiffer | Multi body Dynamics with Unilateral Contacts[END_REF]11) and (13-15), plus the additional conditions:

g3D • S = 0, n1 • su = 0, d1-11 • su = 0, ( 16
) (17) ( 18
)
where eqs. (16-18) correspond to eqs. ( 10), (13), and ( 14), respectively, and impose the conditions for the second tangent to the associated surface. Eq. ( 16) must be repeated for both surface I and surface I I.

A second implementation detail that needs to be considered is related to the fact that constraint equations can usually only be expressed in terms of state variables of the model bodies. However, for the implementation of the minimum distance problem just outlined, the constraint relations also depend on the parametric coordinates which represent additional algebraic unknowns.

The problem can be easily solved introducing two (four in the 3D case) additional massless parts, called "trackers" in the following, associated with the parametric coordinates. The trackers are constrained to the ground by translational joints, in order to have one single degree of freedom left. The value assumed by this displacement at each time instant represents the value of a parametric coordinate along a curve or surface, this way identifying the contact point location.

Having introduced all the constraint conditions just discussed, we have that all the massless parts are completely determined: the dummy parts are constrained to slide on the curves or surfaces so as to be located at each time instant at minimum distance, while the tracker locations correspond to the values of the parametric coordinates at the candidate contact points. It is now possible, based on the solution of the minimum distance problem, to trigger the application of a contact force which can be easily programmed once again through user defined force routines.

Some multi body codes have features that allow the hiding of all the details of the proposed implementation from the user, therefore easing model preparation and verification. In particular, it is usually possible to render the generation of the dummy parts, trackers and constraint equations, fully automatic without user intervention. This will make for a cleaner interface and simpler model construction.

This approach was implemented in the ADAMS software package.

Implementation Through User Defined Force Subroutines

The previously discussed. approach gives a rigorous implementation of the minimum distance problem within a proprietary code using user defined constraint routines. As a matter of fact, the non-linear holonomic constraints that define the minimum distance problem are appended to the other constraint conditions of the multibody system and solved together with all the other equations. However, this comes to the price of introducing the dummy parts that move along the contacting curves or surfaces identifying the candidate contact points, and also requires the introduction of the trackers that represent the algebraic unknowns associated to the parametric coordinates. This somewhat increases the number of degrees of freedom of the system. Furthermore, the solver needs to be able to deal with massless bodies, a feature that not all commercial packages will support.

A possibly simpler strategy implements the minimqm distance problem within a user defined force routine, therefore avoiding the introduction of holonomic constraint equations. In practice, each time the user force routine is called by the solver, one solves within that routine the minimum distance problem given by eqs. {3) and {4) in 2D or by eqs. {5-8) in 3D. Once the non-linear problem has been solved, for example using a quasi-Newton method, one can evaluate the minimum distance q and, in case of contact, apply the corresponding contact force.

Since the contact points are known from the values of the parametric coordinates that satisfy the minimum distance problem, one can at this point evaluate the transport of the contact force to a known pre-determined point on the body, e.g. the center of mass. This way, one avoids completely the need to define massless dummy parts and furthermore one does not need to introduce trackers for letting the solver deal with algebraic unknowns, since these algebraic unknowns are dealt with internally by the user force routine. Consequently, the implementation leads to a smaller number of unknowns and equations. The price to pay for this is that a {small) non-linear problem must be solved at each call of the user force routine, and therefore we have a non-linear problem nested within the non-linear time marching problem. However, the convergence of the inner loop is usually extremely fast, since we can use the converged values at the previous call as a starting point for the iteration.

This second approach was implemented in LMS DADS and in ADAMS.

Numerical Examples Circuit Breakers

The first problem is concerned with the simulation of the closing operation of a circuit breaker. The system is composed of several rigid components and its configuration is depicted in Figure 1. The input energy is provided by a closing spring, that is responsible for both the quick activation time of the mechanism and for the recharging of the opening spring. Ten unilateral contact conditions are present in this example, that exhibits a very complex behavior with short duration impacts as well as prolonged contacts, together with considerable interaction forces.

The problem was ana1yzed using the ADAMS software package, equipped with the unilateral contact modeling capabilities described in this paper. Figure 2 gives a few snapshots of the system configuration during the simulation. Figure 3 presents a plot of the time history of the contact forces for the various unilateral constraints in last part of the simulation. A detailed description of the functionality of the system is beyond the scopes of the present work, however we can notice the typical behavior that is expected in the contact/impact analysis of these mechanisms, i.e. large interaction forces characterized by very sharp gradients that account for the changes in system topology. Close inspection of the plot shows multiple repeated contacts followed by separation of the bodies, or sticking conditions when the bodies remain together for prolonged periods of time.

A second circuit breaker is analyzed next. The mode of operation of the breaker is similar to the previous case, although the topology, geometry and details of the mechanism are different. In this case, we consider the opening of the breaker, induced by a pre-charged spring. During this operation, large intermittent contact forces arise on a number of stops placed both on the fixed frame of the breaker and on the various moving links, in order to restrain their relative motion within predefined values.

This problem was analyzed using the LMS DADS software. The overall configuration of the system is given in Figure 4. Figure 5 shows an animation of the breaker during the simulation (some parts of the mechanism are not represented for clarity). Even in this case we have several violent impacts with the development of large contact forces, as depicted in Figure 6. Figure 7 shows the time history of the minimum distance between a link and a fixed stop, illustrating the repeated contacts that take place in a few milliseconds.

Recharging Mechanism

The example deals with the transient simulation of a recharging mechanism. The model is composed of several rigid body components, in particular a wheel connected to a torsional spring, and an oscillating mechanism that engages the wheel geometry through a tooth.

When the oscillating part of the mechanism moves in one direction, the tooth pushes on the wheel profile and makes it turn. When the oscillating part starts moving in the opposite direction, the tooth slides on the wheel as the wheel itself is engaged by another arresting tooth that prevents it from spinning in the opposite direction.

The problem involves multiple impacts of the pushing and arresting teeth with the wheel, and sliding of the same teeth on the wheel profile. The complex shapes. of the wheel and of the teeth are easily handled through the NURBS parameterization. Figure 8 gives a few snapshots of the system configuration throughout the simulation, conducted using the ADAMS software.

Conclusions

In this work we have considered the problem of implementing unilateral contact modeling capabilities in commercial multibody codes, comparing two well known methods, namely the impulsive and the continuous force models. We have concluded that the latter is better suited for the implementation through standard user defined subroutines, i.e. without accessing the program source code, since it only necessitates of the definition of non-linear holonomic constraints and of forces that depend on the system state.

We have discussed the details of two possible implementations of the described algorithms.

In the first approach, one appends the non-linear holonomic constraints that define the minimum distance problem to the constraint equations already present in the multibody system. For this to work, one must introduce massless dummy parts that are constrained to move on the contact curves or surfaces and that identify the candidate contact points at each instant of the simulation. Furthermore, one has also to introduce massless trackers that define the values of the algebraic unknowns entering the problem. In the second approach, one avoids the above mentioned complications by solving the non-linear minimum distance problem within each call to a user defined force routine.

We have successfully implemented the described procedures in two popular commercial packages, namely ADAMS and LMS DADS. These additions to the codes allow to model intermittent motion of rigid and flexible multibody systems. For ensuring a wide spectrum of applicability to the method, we have equipped the procedures with a library implementing a parameterization of the contacting parts based on NURBS. Furthermore, we have been able to provide a user-friendly interface that hides all the unnecessary details that could increase the user's workload.

After debugging on a set of simple test cases, :we have tested the procedures with a number of problems of industrial relevance. The numerical simulations demonstrate the capabilities of the proposed approach. In particular, the tests conducted have shown that the time adaption algorithms of the target codes can deal effectively even with the strong requirements imposed by contact/impact analysis of multibody systems. The two proposed implementations are both rigorous by treating exactly the minimum distance problem, therefore they yield virtually identical solutions to the problems here analyzed.

  nr(~r,'TJr) •su(~u,'TJII) = 0, dr-II(6, 'TJI, 6r, 'TJII). SII(6r, 'T/II) = 0, nr(~r,'TJr) • tii(6r,'TJir) = 0, dr-II(6, 'TJI, (IJ, 'T/II) • tu(~II, 'T/II) = 0.

  v (i = 1, 3) are the dummy part unit vectors, and:

Figure 2 :Figure 3 :

 23 Figure 2: Snapshots of the circuit breaker configuration during the simulation.

Figure 4 :Figure 5 :

 45 Figure 4: View of the second circuit breaker.

ADAMS is a trademark of Mechanical Dynamics, Inc.

LMS DADS is a trademark of LMS International.

Acknowledgments

This research was supported by ABB Ricerca under contract CRPOLIMI 17089/1/AS/99 with Politecnico di Milano. This support is gratefully acknowledged.