Xiao Hu
email: xiao.hu@recherche.enac.fr

Paul Thevenon

Christophe Mababiau

Improving reliability and efficiency of RTK ambiguity resolution using multiple rover receivers connected to the same antenna

Precise position information is significantly required for specific application scenarios like autonomous driving of vehicles or precise mobile mapping. Until recently, to satisfy the sub-meter or even centimeter-level demanding for this type of application, the GNSS carrier phase measurements appear mandatory. In this paper, we firstly use multiple rover receivers to enhance the performance of the RTK positioning in different environments by connecting them to the same antenna thus multiplying the number of observations by the number of receivers. We demonstrate through different simulation scenarios that we are able to improve some internal steps of precise position computation and the simulations results show that our multi-receiver RTK system is more robust to noise and degraded satellite geometry, in terms of ambiguity fixing rate, and get a better position accuracy under same conditions when comparing to the single receiver system. An analysis is also conducted by varying the correlation coefficient of the noise between the measurements to the same satellite performed by the different receivers to verify the advantage of our proposed system when dealing with highly correlated noise.

Introduction

GNSS receivers have become a mass-market device by millions of users every day. Current positioning accuracy is usually sufficient to lead the way of a car driver into an unknown area or provide a stable timing signal to a mobile cell tower. However, stand-alone positioning technique is not precise enough for applications requiring sub-meter to centimeter-level accuracy such as precise mobile mapping and vehicles automatic control. To achieve this kind of level accuracy of positioning, the GNSS carrier phase measurements appear compulsory. Nevertheless, they are also less robust and include a so-called integer ambiguity resolution (IAR) that prevents them to be used directly for positioning. In some harsh environments, severe code multipath, frequent cycle slips, and losses of the lock might also challenge the process of finding the correct IAR thus reduce the accuracy.

To increase confidence and speed up the IAR process by limiting the search space, constraints can be established by using an array of two or more receivers with prior known and fixed geometry, including the length of the baseline vectors between the receivers of the array and the orientation of the vectors.

Recently, several studies have focused on the use of multiple receivers for attitude determination [START_REF] Iafrancesco | GPS/INS Tightly coupled position and attitude determination with low-cost sensors Master Thesis[END_REF], [START_REF] Henkel | Attitude determination with low-cost GPS/ INS[END_REF] and calibration of magnetic field sensors [START_REF] Henkel | Calibration of magnetic field sensors with two massmarket GNSS receivers[END_REF], Daniel et al. [START_REF] Medina | On the Kalman filtering formulation for RTK joint positioning and attitude quaternion determination[END_REF] developed a method for the recursive estimation of the positioning and attitude problems using GNSS carrier phase observations from an array of receivers, but they calculated the position of each receiver separately thus they didn't take advantage of the known geometric constraints. Fan et al. [START_REF] Fan | Precise and Robust RTK-GNSS Positioning in Urban Environments with Dual-Antenna Configuration[END_REF] proposed a dual-antenna constraint RTK algorithm, which combines GNSS measurements of both antennas by making use of the geometric constraint between them, the reception diversity improves the availability and geometric distribution of GNSS satellites thus improves the positioning accuracy. However, their algorithm is severely restricted by the environmental conditions. Zheng et al. [START_REF] Zheng | Enhancing Ambiguity Resolution Performance Using Attitude Determination Constraints[END_REF] presented a methodology for integrating carrier phase attitude determination and positioning systems by considering one of the receiver pairs in the attitude determination system also used as the rover for the relative positioning system. Nevertheless, their positioning systems remained independent which did not much ameliorate the success rate of IAR for the RTK positioning. Nandakumaran et al. [START_REF] Nadarajah | Multi-GNSS PPP-RTK: From Large-to Small-Scale Networks[END_REF] provided a numerical insight into the role taken by the multi-GNSS integration in delivering high-precision solutions, however, they focused on the PPP solution rather than the RTK processing which requires long convergence time to reach cm-level accuracy.

As far as we know, the use of an array of receivers is not used to improve the accuracy of the array position or for some internal steps of precise position computation. The concept of RTK positioning and attitude determination using an array of GNSS receivers with known geometry has already been introduced in [START_REF] Hu | Improvement of RTK performances using an array of receivers with known geometry[END_REF]. In [START_REF] Hu | Cycle-slip Detection and Repair Using an Array of Receivers with Known Geometry for RTK Positioning[END_REF], we discussed the improvement of cycle-slip detection and repair for RTK processing by using an array of two receivers with known geometry.

In this contribution, we will focus on one specific scenario which is the improvement of RTK processing by using an array of two single-frequency receivers with zero baseline constraint which can be realized by connecting the receivers to the same antenna. The approach presented here aims at developing a method that takes full advantage of the measurement redundancy and diversity, along with the geometric constraints, between two closely placed rover antennas to improve the accuracy of the float ambiguity estimation and, consequently, the AR success rate.

An analysis is also conducted by varying the correlation coefficient of the noise between the measurements to the same satellite performed by the different receivers to verify the advantage of our proposed system when dealing with highly correlated noise between multiple rover receivers.

Mathematical Models

In this work, a Kalman Filter (KF) based position determination algorithm is developed. The scenario considers a moving vehicle, on which two low-cost, single-frequency receivers are connected to one single antenna, whereas the third one serves as a reference station, with its position static and known.

In the following, the state and the measurement vectors are described along with the statetransition and measurement models.

State transition model

The state transition or state-space model describes how the states or parameters of the system vary over time based on a specific linear model. In our KF modeling, the state parameter transition between subsequent epochs is given by:

= + 1
where: refers to the state vector at epoch refers to the transition matrix of epoch -1 refers to the so-called system noise vector at epoch

Together with the process noise vector, one can define the process noise covariance matrix as:

= T 2 3
This matrix has then the variances of the state parameter's estimates based on the system model.

The estimated parameters are collected inside the state vector. The state vector collects 3 vehicle state parameters and 2*( -1) satellite state parameters which are: the 3D position of GNSS receiver 1 relative to GNSS receiver 3, double-difference phase ambiguity of the GNSS receiver pair 1-3, and the double-difference phase ambiguity of the receiver pair 2-3. which is given by: = 3

where:

refers to the vector of the double-difference GPS integer ambiguity between receiver and refers to the relative position vector between reference receiver 3 and rover receiver 1

Transition model for position related state parameters

In our KF modeling, for the position related state parameters, we suppose the speed components are a zero-mean Gaussian process, which can be defined as:

, = , + , 4
where: is a centered Gaussian vector with a covariance matrix

Transition model for satellite-related state parameters

In the case of the satellite-related parameters, they are assumed all as constant over subsequent epochs with a very small noise compared to the position related state parameters.

, =

, is a centered Gaussian variable with a covariance matrix = is the identity matrix of size The resulting state transition matrix is then given by a unit matrix and different values of process noise variance are added to complete the model.

= 1 ⋯ 0 ⋮ ⋱ ⋮ 0 ⋯ 1 6
And the corresponding process noise matrix is given as follows:

= × × × × 7

Measurement model

The measurement model describes how the individual sensor measurements are related to system states. For every epoch n, the measurement vector , which contains all measured values, can be described as a function of the state vector as: = ( ) + 8 with is the function that relates one or more states with each measured value and is the measurement noise vector, which describes the expected Gaussian noise of every measured value.

As for the process noise covariance matrix, the definition of the measurement noise covariance matrix follows as:

= T 9
In our model, the measurement vector comprises the following measured values: Doubledifference (DD) code phase measurement vector of receiver 1, DD code phase measurement vector of receiver 2, DD carrier phase measurement vector of receiver 1, and DD carrier phase measurement vector of receiver 2.

= ( ) ( ) ( ) ( ) 10 In this measurement model, the position of receiver 2 is assumed to share the same position as receiver 1. The individual double-difference corrected pseudo-range and phase GPS measurement can be modeled as:

= ( ) + , = ( ) + , = ( ) + + , 11 = ( ) + + , where:
is the double-difference code phase measurement of receiver pair and satellite pair is the double-difference carrier phase measurement of receiver pair and satellite pair is the difference between the Line of Sight vector of satellite and l is the double-difference integer ambiguity of receiver pair and satellite pair , , , is the noise measurement of the double-difference code and phase measurement respectively

The LOS vectors are modeled by using the azimuth and elevation of the satellites with respect to the corresponding receiver, thanks to the following formula who define the LOS vector between satellite k and the receiver in the ENU frame: = cos cos cos sin sin 1 2

A fixed measurement noise variance between all satellites is defined to complete the measurement model, defining the measurement covariance matrix .

As the measurements are coming from signals received by the same antenna, it is safe to consider a certain level of correlation between these noise measurements. For example, the multipath error may be similar on the measurement performed by the 2 receivers from the same satellite.

The covariance coefficient between the 2 receivers is given by

= ( 13
)
Where is the covariance of the noise affecting the measurement performed by receivers 1 and 2 on satellite is the correlation coefficient of the noise affecting the measurement performed by receivers 1 and 2 is the standard deviation of the noise affecting the measurement performed by receiver on satellite

We consider a simplifying assumption that the noise standard deviation does not depend on the satellite or the receivers ( = ), and that the code measurement noise standard deviation and the phase measurement one are related by a simple factor = * , with a value k equals to 1/100 to account for the much better accuracy of carrier phase measurements. The measurement covariance matrix for the single difference measurements will have the following shape: The measurement covariance matrix for the double-difference measurements can then be deduced using the following formula:

= , , (16)
Where

= -1 0 ⋯ 0 1 0 -1 ⋱ ⋮ 1 ⋮ ⋱ ⋱ 0 ⋮ 0 ⋯ 0 -1 1 (17)
is the single-differencing matrix used for computing the double-difference.

Two alternating steps which are the state prediction step and the state update step are then conducted to complete the KF algorithm.

RTK processing

In the previous process, we obtain a float estimation of the double-difference integer ambiguity. The accuracy of the position state estimate is further improved by fixing the DD ambiguities to integer number by using the well-known LAMBDA algorithm [START_REF] Teunissen | The least-squares ambiguity decorrelation adjustment: a method for fast GPS integer ambiguity estimation[END_REF], [START_REF] Buist | The Baseline Constrained LAMBDA Method for Single Epoch, Single Frequency Attitude Determination Applications[END_REF]. This algorithm has shown its superiority in terms of both performance and processing time when compared to other algorithms in [START_REF] Baroni | Analysis of Attitude Determination Methods Using GPS Carrier Phase Measurements[END_REF].

One selects the integer candidates based on the sum of squared errors to get a fixed solution.

The candidate with the lowest error norm is chosen once the ratio of the Maximum A Posteriori error norm between the second-best candidate and the best candidate is bigger than a threshold. It is a pre-defined threshold or the critical value that the squared norm of ambiguity residuals of the best and second-best candidates should overpass to validate the integer estimation. In our paper, an empirical fixed value 3 is taken as in [START_REF] Teunissen | On the Foundation of the Popular Ratio Test for GNSS Ambiguity Resolution[END_REF].

Once the IAR process is declared successful, a new position is computed using the DD carrier phase measurements corrected by the validated DD integer ambiguities. This final position is the fixed solution.

Simulation of RTK performances using multiple rover receivers connected to the same antenna

In this section, the proposed precise position determination algorithm is verified with simulated measurements. The reason why we use the simulated measurements is that we want to analyze the advantage of multiple receivers connected to the same antenna when varying the amplitude of the measurement noise, the level of noise correlation between the rover receivers' measurements and the satellite geometry.

The vehicle is assumed to make a uniform circular motion around the fixed reference receiver 3 with a constant radius (R = 100 m) in all the scenarios. The speed of the vehicle is 10 m/s, the duration of the simulated trajectory is 1000 s and the sampling rate is 1 Hz.

In our simulator, the satellites are simulated with a fixed position. We consider 7 satellites with a fairly good geometric condition, i.e. well-spaced in azimuth and elevation. This would correspond to a typical case of a single constellation in an open sky environment. Finally, the DD carrier ambiguities are taken as fixed values separated by 10 cycles. This means that the DD ambiguities go from 10 to 120.

Typical results of the state estimations including the estimation of the integer ambiguity by using LAMBDA are given in fig. 1. We can see from the fig. 2 that the same convergence period is present in all states for the fixed solution. In the case of the float solution, the East, North orientations converge rapidly to the true value, however, due to the observability issue, the altitude takes more time to converge. What's more, once they converge, the error between the estimated result and the true value is extremely small especially in the case of the fixed solution. From the figure, one can notice that the system fixes the IAR in the 64th epoch and once the IAR is fixed, the positioning results become much more accurate compared to the float solution.

To validate the influence of the using of an array of receivers, several scenarios are conducted.

Results are compared in terms of fixed rate and 3D position error statistics (mean, standard deviation, and 95% bound). The numbers have been computed by running 100 simulations with the different random realization of noise and initialization values and taking the average of the investigated indicators.

Scenario 1: Geometry of satellite effect Analyses

A favorable satellite geometry is essential to get an accurate positioning result. That is why in the first scenario, we want to analyze the influence of the system when the geometry of the satellite changes. In this scenario, both the correlation level ( = 0 ) and the amplitude of the noise measurement ( = 1 m) are fixed and by degrading the geometry of the satellite (by removing one satellite at a time), we observe the robustness against degraded satellite geometry of multiple rover receivers RTK systems. The different scenarios are shown in fig. The analysis of Table 1 shows that the use of the receiver array improves the fixed rate for all considered satellite geometry, besides, the dual receiver array provides always better performance than the single receiver RTK solution under the same satellite geometry situation which confirms the advantage of our approach.

Scenario 2: Noise measurement effect Analyses

In the second scenario, the correlation level of the noise measurement ( = 0 ) and the satellite geometry ( = 7) are fixed, by varying the amplitude of the noise measurement, we observe the robustness against the noise of the multi-receiver RTK system.

We recall that the phase measurement noise standard deviation is scaled by a factor 1/100, to account for the improved accuracy of the carrier phase measurement, compared to the code measurement. The analysis of Table 2 shows that the use of the receiver array improves the fixed rate for all considered values, thus demonstrating the interest of such an approach.

Scenario 3: Noise measurement Correlation effect Analyses

In probability theory and statistics, correlation (or correlation coefficient) is used as a measure of the linear correlation between two variables and . The most widely used is the Pearson product-moment correlation coefficient (PPMCC). It is defined as the covariance of two variables divided by the standard deviation of the two variables as shown in ( 17)

, = cov( , ) = ( -)( -) (17)
The significance of the magnitude of the correlation coefficient is currently inconsistent in the statistical world, but it is usually considered as follows: correlation of these noise measurements to make it closer to the real situation. As an example, the multipath error may be similar on the measurements performed by the 2 receivers from the same satellite.

Consequently, in this scenario, we want to verify what level of correlation is acceptable to keep an improvement compared to a single receiver RTK system. The correlation level of the noise measurement is not fixed, by varying the noise correlation coefficient while keeping the same noise measurement amplitude ( = 1 m) and a favorable satellite geometry ( = 7). The different scenarios of the magnitude of the correlation coefficient are shown in Fig. The result is shown in Table 4, from the result one can conclude that the dual receiver RTK system provides always better performance than the single receiver RTK solution and the position accuracy sightly improves as the magnitude of the correlation coefficient decreases. In addition, the fixing rate and positioning accuracy are not significantly reduced even when dealing with highly correlated noise, which proves that our system is exceedingly robust and stable when handling high correlation noise situations. This may be due to the fact that the LAMBDA algorithm takes into account the correlation of the measurements in the IAR process, which is assumed perfectly known in our simulation.

Conclusion

In this contribution, the method using multiple receivers connected to the same antenna (zerobaseline configuration) for the RTK positioning in different environments is studied. Taking advantage of the doubled quantity of the observations, we are able to improve some internal steps of precise position computation. We demonstrate through simulations results that our multi-receiver RTK system is more robust to noise and degraded satellite geometry, in terms of ambiguity fixing rate, and get a better position accuracy under the same conditions when comparing with the single receiver system. Furthermore, the mitigation of high correlation noise between receivers is permitted by using our system.

Figure 1 .

 1 Figure 1. Illustration of state estimation including IAR Fig. 2 gives the error comparison between the float and fixed solution. Here the fixed solution is defined by two parts, for every epoch, if the IAR process is declared successful, a new position is computed using the DD carrier phase measurements corrected by the validated DD integer ambiguities. This position then is used to update the fixed solution. If the IAR process is not declared successful, we keep the float solution.

Figure 2 .

 2 Figure 2. Illustration of float vs fixed solution

3

 3

Figure 3 .

 3 Figure 3. considered satellite geometry with = 7 (left), = 5 (middle) and = 4 (right)

4

 4

Figure 4 .

 4 Figure 4. considered noise correlation distribution

Table 1 -Performance comparison for different satellite geometry,

 1

	= ,	= m

Table 2 -Performance comparison for different measurement noise,

 2

	= ,	=

Table 3 -Significance of the magnitude of the correlation coefficient

 3 As we discussed in the previous section, in our model, the measurements are coming from signals received by the same antenna, thus it is paramount to consider a certain level of

	0.00-± 0.30 slightly correlated	± 0.30-± 0.50 medium correlation
	± 0.50-± 0.80 significantly correlated	± 0.80-± 1.00 highly correlated

Table 4 -Performance comparison for different measurement noise,

 4 For verifying the sensitivity of the LAMBDA algorithm to a correct correlation coefficient, the final simulation has been conducted, where the correlation coefficient considered in the algorithm is different from the one used to generate the measurement ( = 0.4 for all the simulations), a similar result is obtained as shown in the last line of table 4.

	=	,	=


Dr Christian Laugier
email: christian.laugier@inria.fr

Increasing Driving Safety Perception & Assessment of Collision Risks with other Road Users

•©•

  Technological breakthrough & Numerous AV experimentations in real traffic conditions Tesla Autopilot L2 (ADAS) • 32 millions km covered since 2009 • 1 st US self-driving taxi service L3 (Phoenix, since Dec 2018) Autonomous Shuttle service • Millions of miles driven last decade… but SAFETY is still not fully guaranteed ! => Several benign or fatal accidents involving AVs (Perception failure & Takeover difficult) Robot taxi Tesla Autopilot level 2 (May 2016) => Perception & Decision-making technologies have still to be improved for mixt traffic !!! Self-driving Uber L3 (March 2018) C. LAUGIER -Confiance numérique: le cas du véhicule autonome AG IRT-Nanoelec -Table ronde "Une filière pour une société résiliante", Minatec Grenoble, Sept 8 th 2020 3 Dynamic & Open Environments, Incompleteness & Uncertainty, Sensors limitations, Real-time + Validation  Mixed traffic (Human in the loop) => Human Aware Decision-making process Taking into account Interactions + Behaviors + Social & Traffic rules Road Safety Campaign, France 2014 Situation Awareness & Decision-making => Safe intentional navigation (using semantics) Dealing with unexpected events => Avoiding upcoming collisions with "something"  Exploiting the Dynamic information for a better understanding of the scene !!!!  Reasoning about Uncertainty & Time window => Past & Future predicted events  Bayesian Sensors Fusion + Scene interpretation using Contextual Semantic information 1 st Paradigm: -Confiance numérique: le cas du véhicule autonome AG IRT-Nanoelec -Table ronde "Une filière pour une société résiliante", Minatec Grenoble, Sept 8 th 2020 5 © Laugier & Inria. All rights reserved  Predict environment changes on a given "time horizon t+d" => Using History & Motion models  Estimate the Probabilistic Risk of Collision at t+d (d = a few seconds ahead)  Make Driving Decisions by taking into account the Predicted behavior of all surrounding traffic participants (cars, cycles, pedestrians…) & Social / Traffic rules (traffic participants interactions) Yellow => time to collision: 3s • Orange => time to collision: 2s • Red => time to collision: 1s Several PhD theses 3 rd Paradigm: Models improvements using Machine Learning  Perception level: Construct "Semantic Grids" using Bayesian Perception & DL  Prediction & Decision-making level: Learn driving skills for Autonomous Driving  Open questions: Training step (Available Datasets limited), Real-time processing (difficult), Classification Errors (often not explainable), Domain adaptation (e.g. changing weather conditions) Step: Modeling Driver Behaviors using IRL o 2 nd Step:  Predict behaviors of surrounding vehicles (using Perception & learned Behavior models)  Make "safe & consistent" Driving Decisions for Ego Vehicle Ego vehicle Front cam Ego vehicle Back cam  White vehicle => Ego-vehicle (ground-truth)  Red box => Plan induced (Predicted trajectory)  Yellow boxes => Detected obstacles (CMCDOT) C. LAUGIER -Confiance numérique: le cas du véhicule autonome AG IRT-Nanoelec -Table ronde "Une filière pour une société résiliante", Minatec Grenoble, Sept 8 th 2020 7 © Laugier & Inria. All rights reserved Concluding remarks & Discussion  Increasing impact of AI + Real-time data processing capacity + Increased sensor performance + New Models & Embedded algorithms + Multiplication of tests in real conditions => The unmanned car is gradually becoming a technological reality  Safety is not yet fully guaranteed ! o Current Perception & Scene Understanding algorithms are not robust enough for complex & highly dynamic environments o Need to take better account of Interactions with other road users (using also AI approaches) o Need to develop Validation & Certification Tools and Methodologies => Realistic simulators, Realworld testing protocols, Formal methods (e.g. Enable-S3 EU project & future French project Prissma)  User confidence & Acceptance by the human society will be decisive to allow a real deployment (e.g. "cohabitation" with other users such as pedestrians, bicycles, scooters ...) o Autonomous vehicles a priori safer than cars driven by humans (inattention)…. but 0 tolerance in the event of a fatal accident involving an autonomous vehicle ! o Ethics & Responsibility issues must also be taking into account before any deployment

© Laugier & Inria. All rights reserved


Christian Laugier

Embedded Perception & Decision-making for Safe Navigation in Uncertain, Dynamic and Human-populated Environments

Keywords: Uncertainty, Partial Knowledge, World changes, Real

 Estimate Probabilistic Collision Risk at a given time horizon t+d (d = a few seconds)  Make Driving Decisions by taking into account the Predicted behavior of all the observed surrounding traffic participants (cars, cycles, pedestrians …) & Social / Traffic rules 2 nd Paradigm: Collision Risk Assessment & Decision-making Decision-making for avoiding Pending & Future Collisions Complex dynamic situation Risk-Based Decision-making => Safest maneuver to execute Alarm / Control Human-Aware Situation Assessment  Decision-making: Two types of "collision risk" have to be considered  Short-term collision risk => Imminent collisions with "something" (unclassified), time horizon <3s, conservative hypotheses  Long-term collision risk => Future potential collisions, horizon >3s, Context + Semantics, Behavior models

•©

  Strong involvement of Car Industry & GAFA + Large media coverage + Increasing Governments supports • An expected market of 515 B€ at horizon 2035 (~17% world automobile market, Consulting agency AT Kearney, Dec 2017 ) • But Validation & Certification issues are still unclear … idem for Legal & Regulation issues => Numerous experiments in real traffic conditions since 2010 (Disengagement reports  Insights on system maturity) => But still insufficient … Realistic Simulation & Formal methods are also under development (e.g. EU Enable-S3) "Self-Driving Taxi Service L3" testing in US (Uber, Waymo) & Singapore (nuTonomy)  Autonomous Mobility Service, Numerous Sensors +"Safety driver" during testing (take over in case  Uber: System testing since 2017, Disengagement every 0.7 miles in 2017 (improved now)  Waymo: 1 st US Self Driving Taxi Service launched in Phoenix in Dec 2018  Disengagement reports provide insights on the technology maturity EU CityMobil project & Inria Drive Me trials (Volvo, 2017) • 100 Test Vehicles in Göteborg, 80 km, 70km/h • No pedestrians & Separations between lanes Numerous EU projects in last 2 decades Cybus, 3 months experiment, low speed La Rochelle 2012 (France) Dense 3D mapping & Numerous vehicles 10 years R&D, 8 millions km covered since 2010 & 25 000 km/day Fatal accidents involving AVs -Perception failure  Tesla driver killed in a crash with Autopilot "level 2" active (ADAS mode) -May 2016 Tesla Model S Autopilot  The Autopilot failed to detect a white moving truck, with a brightly lit sky (Camera Mobileye + Radar)  The human driver was not vigilant & didn't took over  Self-driving Uber L3 vehicle killed a woman => First fatal crash involving a pedestrian Temple, Arizona, March 2018  Despite the presence of multiple sensors (lidars, cameras …), the perception system failed to detect the pedestrian & didn't disengaged  The Safety Driver reacted too lately (1s before the crash) C. LAUGIER -Embedded Perception & Decision-making for Safe navigation in Uncertain, Dynamic and Human environments Matinée Coboteam-Inria sur la navigation robotique, May 28 th 2020 4 © Laugier & Inria. All rights reserved Challenge 1: The need for Robust, Self-diagnosing & Explainable Embedded Perception Video Scenario: • The Tesla perception system failed to detect the barriers blocking the left side route • The driver has to take over and steer the vehicle away from the blocked route (for avoiding the collision) AVs have to face two main challenges Video source: AutoPilot Review @ youtube.com Challenge 2: The need for Understandable Driving Decisions (share the road with human drivers) Video scenario: Scene observed by the dash cam of a bus moving behind the Waymo AV • Waymo AV is blocked by an obstacle and it decides to execute a left lane change • The bus driver misunderstood the Tesla's intention and didn't yield • The two vehicles collided Video source: The Telegraph Human drivers actions are determined by a complex set of interdependent factors difficult to model (e.g. intentions, perception, emotions …)  Predicting human driver behaviors is inherently uncertain  AV have to reason about uncertain intentions of the surrounding vehicles C. LAUGIER -Embedded Perception & Decision-making for Safe navigation in Uncertain, Dynamic and Human environments Matinée Coboteam-Inria sur la navigation robotique, May 28 th 2020 5 © Laugier & Inria. All rights reserved Embedded Bayesian Perception & Decision-making Main Challenges & Required Technological Breakthrough => Robustness, Efficiency, Dynamic Human Environments (Safety is still not guaranteed) => Real-time integration of Perception & Motion planning & Control => Integration into Embedded Hardware & Software (future products) => Validation & Certification ADAS & Autonomous Driving C. LAUGIER -Embedded Perception & Decision-making for Safe navigation in Uncertain, Dynamic and Human environments Matinée Coboteam-Inria sur la navigation robotique, May 28 th 2020 6 © Laugier & Inria. All rights reserved Situation Awareness & Decision-making  Sensing + Prior knowledge + Interpretation  Selecting appropriate Navigation strategy (planning & control) Dynamic & Open Environments => Real-time processing & Reactivity (several reasoning levels are required)  Incompleteness & Uncertainty => Appropriate Model & Algorithms (probabilistic approaches)  Sensors limitations (no sensor is perfect) => Multi-Sensors Fusion  Hardware / Software integration => Satisfying Embedded constraints  Human in the loop (mixed traffic) => Human Aware Decision-making process (AI based technologies) Taking into account Interactions + Behaviors + Social rules (including traffic rules) for avoiding upcoming collisions with "something" => High reactivity & reflexive actions => Focus of Attention & Sensing => Collision Risk estimation + Avoidance strategy Dealing with unexpected events Road Safety Campaign, France 2014 C. LAUGIER -Embedded Perception & Decision-making for Safe navigation in Uncertain, Dynamic and Human environments Matinée Coboteam-Inria sur la navigation robotique, May 28 th 2020 7 Noisy data, Incompleteness, Dynamicity, Discrete measurements  Strong Embedded & Real time constraints  Our Approach: Embedded Bayesian Perception  Reasoning about Uncertainty & Time window (Past & Future events)  Improving robustness using Bayesian Sensors Fusion  Interpreting the dynamic scene using Contextual & Semantic information  Software & Hardware integration using GPU, Multicore, Microcontrollers…

© Laugier & Inria. All rights reserved


A Kheirati Roonizi

Christian Jutten

Band-stop Smoothing Filter Design

Arman Kheirati Roonizi * , and Christian Jutten, Fellow, IEEE Abstract-Smoothness priors and quadratic variation (QV) regularization are widely used techniques in many applications ranging from signal and image processing, computer vision, pattern recognition, and many other fields of engineering and science. In this contribution, an extension of such algorithms to band-stop smoothing filters (BSSFs) is investigated. For designing a BSSF, the most important parameters are the order and the cutoff frequencies. In this paper, we show that with the optimization approaches (smoothness priors or QV regularization), the cutoff frequencies are related to the regularized parameters and the order can be directly (and easily) controlled with the number of derivatives. We describe two ways to implement the BSSFs using these approaches. First, we present a parallel structure to BSSF and then illustrate why it is less than ideal. Next, we present a novel approach regarding parallel structure to produce BSSFs with very sharp transition bands for high-performance applications. An improved optimization-based approach to BSSF design is introduced. The performance of the new BSSFs is nearly ideal.

Index Terms-Smoothness priors, Quadratic variation regularization, Least-squares optimization, Band-stop smoothing filter, Parallel structure

I. INTRODUCTION

T HERE is no doubt that smoothing filters (i.e., smoothness priors and quadratic variation) have a respectable place in signal and image processing, computer vision and pattern recognition, statistics, time series analysis, economics and many other fields of engineering and science. Many papers and a number of books have appeared. In the following, some major contributions that have been developed in the design of smoothing filters are highlighted.

The conceptual structure of smoothing filter dates from 1899 with Bohlmann's work who studied the application of regularized method for time series smoothing [START_REF] Bohlmann | Ein ausgleichungsproblem[END_REF]. In 1923, Whittaker addressed the problem of estimating a smooth trend embedded in white noise and spoke about the method of graduating data [START_REF] Whittaker | On a new method of graduation[END_REF]. Therefore, his method is known as the method of graduating data [START_REF] Aitken | On the theory of graduation[END_REF]. At the same time, Henderson had presented a different solution which was widely extended in North America [START_REF] Henderson | A new method of graduation[END_REF]. So the method of graduating data is also known as Whittaker-Henderson (WH) graduation [START_REF]Mathematical theory of graduation[END_REF]. In 1973, Shiller worked on the same problem and introduced the notation of "smoothness priors" [START_REF] Shiller | A distributed lag estimator derived from smoothness priors[END_REF]. It seems that he was not aware of either Whittaker-Henderson graduation or Bohlmann's work. Subsequently, smoothness priors became popular in research communities. It was used in many applications.

In time series analysis, smoothness priors was used for modeling nonstationary time series [START_REF] Kitagawa | A smoothness priors-state space modeling of time series with trend and seasonality[END_REF], [START_REF]A smoothness priors time-varying ar coefficient modeling of nonstationary covariance time series[END_REF], generalizing a nonparametric estimator [START_REF] Shiller | Smoothness priors and nonlinear regression[END_REF], [START_REF] Eubank | A note on smoothness priors and nonlinear regression[END_REF] and exploring large-scale time series [START_REF] Kitagawa | Smoothness prior approach to explore mean structure in large-scale time series[END_REF]. Some developments, extensions and the leastsquares computational framework of smoothness priors in time series analysis are presented in [START_REF] Kitagawa | The Smoothness Priors Concept[END_REF], [START_REF] Brillinger | New Directions in Time Series Analysis[END_REF].

In signal and image analysis, smoothness priors was used for signal smoothing and detrending [START_REF] Tarvainen | An advanced detrending method with application to hrv analysis[END_REF]- [START_REF] Mao | Spatio-temporal signal recovery based on low rank and differential smoothness[END_REF], image smoothing and restoration [START_REF] Xu | Image smoothing via l0 gradient minimization[END_REF], [START_REF] Thompson | A study of methods of choosing the smoothing parameter in image restoration by regularization[END_REF], spectral estimation [START_REF] Kitagawa | A smoothness priors long ar model method for spectral estimation[END_REF] and spline smoothing [START_REF] Shiller | Smoothness priors and nonlinear regression[END_REF], [START_REF] Eilers | Flexible smoothing with b-splines and penalties[END_REF].

In pattern recognition, the second-order priors was used for surface reconstruction [START_REF] Terzopoulos | Multilevel computational processes for visual surface reconstruction[END_REF], [START_REF] Grimson | From Images to Surfaces: A Computational Study of the Human Early Visual System[END_REF], dense stereo [START_REF] Horn | Robot Vision[END_REF] and global stereo reconstruction [START_REF] Woodford | Global stereo reconstruction under second-order smoothness priors[END_REF].

Some other connections to smoothness priors are quadratic variation (QV) regularization [START_REF] Fasano | Baselinewander removal for bioelectrical signals by quadratic variation reduction[END_REF], Hodrick-Prescott filtering [START_REF] Hodrick | Postwar u.s. business cycles: An empirical investigation[END_REF], Savitzky-Golay filter [START_REF] Savitzky | Smoothing and differentiation of data by simplified least squares procedures[END_REF] and the ill-posed inverse problems of statistical Tikhonov regularization [START_REF] Golub | Tikhonov regularization and total least squares[END_REF].

The smoothing approach based QV regularization or smoothness priors is based on a penalized least-squares optimization which depends on a regularization factor that needs to be designed and weighted properly for getting good performance. There are several methods for designing the regularization factor and its weight (usually denoted as an hyperparameter). Some of them are briefly presented below.

• The L-curve method is a well-known heuristic method for choosing the regularization parameter for ill-posed problems by selecting it according to the maximal curvature of the L-curve [START_REF] Boyd | Convex Optimization[END_REF]. • Generalized cross validation is a popular tool for calculating the parameters in the context of inverse problems and regularization techniques. It was initially used to tune the smoothness parameter in ridge regression and smoothing splines [START_REF] Golub | Generalized cross-validation as a method for choosing a good ridge parameter[END_REF]. • Discrepancy principle is a posteriori parameter choice strategies for Tikhonov regularization for solving nonlinear ill-posed problems [START_REF] Engl | Discrepancy principles for tikhonov regularization of ill-posed problems leading to optimal convergence rates[END_REF]. • Maximum likelihood estimation (MLE) algorithm which is based on the concept of Bayesian likelihood, was first introduced by Akaike to determine the smoothness tradeoff parameter [START_REF] Akaike | Factor analysis and AIC[END_REF]. • Stein's unbiased risk estimate (SURE) regularization is another approach that uses the Jacobian matrix of the reconstruction operator with respect to the data [START_REF] Eldar | Generalized sure for exponential families: Applications to regularization[END_REF]. • Regularization factor design in terms of the cutoff frequency. Recently, we have derived a rather simple closed-form expression for the frequency response of the smoothness priors and proposed a closed-form expression for the regularization factor. The design parameter was calculated in terms of the cutoff frequency [START_REF] Kheirati Roonizi | Improved smoothness priors using bilinear transform[END_REF]. For extracting a signal within a predetermined frequency band, conventional linear filtering (digital filtering) techniques are commonly used [START_REF] Oppenheim | Discrete-Time Signal Processing[END_REF]. There are many papers and many ways for designing digital filters developed over the past decades. This technique is performed through linear time-invariant (LTI) systems which are characterized in time/frequency domain by their impulse/frequency response. The frequency response is the Fourier transform of the system's impulse response.

In [START_REF] Sameni | Online filtering using piecewise smoothness priors: Application to normal and abnormal electrocardiogram denoising[END_REF], [START_REF] Kheirati Roonizi | Forward-backward filtering and penalized least-squares optimization: A unified framework[END_REF], [START_REF] Unser | Recursive regularization filters: design, properties, and applications[END_REF], it is shown that smoothness priors or QV regularization can be designed for linear filtering or smoothing. In other word, there is an equivalence between penalized leastsquares optimization and zero-phase filtering. Especially, in a series of our recent papers, a general approach based on a regularized least-squares optimization framework has been introduced to signal denoising/smoothing [START_REF] Kheirati Roonizi | Improved smoothness priors using bilinear transform[END_REF], [START_REF] Kheirati Roonizi | Forward-backward filtering and penalized least-squares optimization: A unified framework[END_REF], [START_REF] Roonizi | A New Approach to ARMAX Signals Smoothing: Application to Variable-Q ARMA Filter Design[END_REF] that shows there is a strong connection between linear filtering techniques and the current approach based on regularized least-squares optimization. We note that the main difference between them is that the designed filters based on regularized least-squares optimization handles smoothing/denoising directly in the time domain while in the conventional linear filtering techniques, the desired filter is first designed in the frequency domain and then implemented in the time domain using convolution operator. In the following, we briefly present the main results regarding the connection between penalized least-squares optimization and filtering reported in [START_REF] Kheirati Roonizi | Improved smoothness priors using bilinear transform[END_REF], [START_REF] Kheirati Roonizi | Forward-backward filtering and penalized least-squares optimization: A unified framework[END_REF], [START_REF] Roonizi | A New Approach to ARMAX Signals Smoothing: Application to Variable-Q ARMA Filter Design[END_REF].

In [START_REF] Roonizi | A New Approach to ARMAX Signals Smoothing: Application to Variable-Q ARMA Filter Design[END_REF], a general framework based on a regularized leastsquares optimization has been introduced to signal smoothing when the signal is represented by an autoregressive moving average (ARMA) model. As an application, the framework was used for ARMA filter design. It was shown that the framework can be driven from a forward-backward filtering, which is accomplished through LTI system. In [START_REF] Roonizi | A New Approach to ARMAX Signals Smoothing: Application to Variable-Q ARMA Filter Design[END_REF], the design of a variable-Q1 notch filter design has been discussed. Recent comparisons with many classical methods such as simple notch filter or narrow-band notch filter using feedback structure proposed by Pei et. al. [START_REF] Pei | Narrowband notch filter using feedback structure[END_REF] showed the superiority of the current approach based on optimization with smoothing constraints.

In [START_REF] Kheirati Roonizi | Forward-backward filtering and penalized least-squares optimization: A unified framework[END_REF], a framework for unification of the penalized leastsquares optimization and forward-backward filtering scheme was presented. We showed that forward-backward filtering (zero-phase IIR filters) can be presented as instances of penalized least-squares optimization. Against conventional linear filtering that uses an IIR filter in its forward pass (followed by time reversing, filtering the reversed signal with the same IIR filter and finally flipping the result), the approach presented in [START_REF] Kheirati Roonizi | Forward-backward filtering and penalized least-squares optimization: A unified framework[END_REF] uses an FIR matrix filter in its forward pass which has the advantage of being inherently stable. Comparison with classical linear filtering showed the superiority of the regularized approach. Therefore, the regularized leastsquares optimization techniques (e.g., smoothness priors or QV regularization strategies) and zero-phase IIR filters can be conflated into one topic. In this paper, we propose a regularized least-squares optimization approach for BSSF design. The optimization problem used in this work is a linear leastsquares optimization (or a convex optimization) problem. It can be solved by equating its derivative with respect to each unknown variable to zero. The approach is different from linear time invariant (LTI) filter and matrix filter which can be used to design BSSF. For more details on matrix filters, we refer the interested readers to [START_REF] Smith | Introduction to Digital Filters: with Audio Applications[END_REF]. In [START_REF] Smith | Introduction to Digital Filters: with Audio Applications[END_REF], the design and implementation of causal recursive filters in terms of banded matrices are discussed. It is based on the fact that a discrete-time filter is described by the difference equation. In this paper, we consider the problem of filtering finitelength signals, since signal processing problems are generally formulated in terms of finite-length signals, and the developed algorithms are targeted for the finite-length case. In [START_REF] Kheirati Roonizi | Improved smoothness priors using bilinear transform[END_REF], we have shown that the regularized least-squares optimization techniques are also suitable for signal recovering even if the signal of interest is (approximately) restricted to a known frequency band. We introduced a new way to the design and implementation of smoothness priors and QV regularization and discussed its application to low-pass, high-pass and bandpass smoothing filter design. Building on that result, in this paper, we mainly focus on the design of band-stop smoothing filter (BSSF). In [START_REF] Kheirati Roonizi | Improved smoothness priors using bilinear transform[END_REF], we have proposed a simple parallel structure for band-pass smoothing filter (BPSF) or BSSF design. The simple band-stop smoothing filter is made out of a low-pass smoothing filter and a high-pass smoothing filter by connecting the two smoothing filter sections in parallel with each other. As it will be shown in Section III, the simple parallel structure has some limitations for BSSF design. Especially, its cutoff frequencies are not exactly mapped on the arbitrary cutoff frequencies and the amplitude at its center frequency is not exactly zero. Furthermore, it has a shallow frequency transition band which is another limitation of the simple parallel structure. It is known that the filters with steep frequency transition bands can better separate signal and noise components in adjacent frequency bands than filters with a shallow frequency transition band. To solve these issues, we propose an improved parallel structure to BSSF design. In the proposed approach, we consider a modified model and investigate the role of noise correlation which represents the link between two smoothing filter sections. Our results suggest that the noise correlation increases the steepness of the roleoff and boosts the performance. First we briefly review the simple method described in [START_REF] Kheirati Roonizi | Improved smoothness priors using bilinear transform[END_REF] for BSSF design and then show its limitations. Next, we present an improved parallel structure which yields nearly ideal BSSF: its cutoff frequencies are exactly mapped on the arbitrary cutoff frequencies and it has a sharp frequency transition band.

The remainder of the paper is organized as follows. In Section II, the relevant background on smoothness priors or QV regularization is reviewed. In Section III, we present a naive BSSF using smoothness priors or QV regularization approach and show its limited performance. Section IV presents an extension of such algorithms to BSSFs design, which can achieve nearly optimal performance. In Section V, the proposed method is used to electrocardiogram (ECG) signal denoising and compared with conventional filters such as zero-phase Butterworth filter. It is shown that the proposed BSSF outperforms band-stop Butterworth filter. Finally, some concluding remarks are presented in Section VI-B.

Throughout this paper, we use the following notations: boldface uppercase/lowercase letters are used to denote matrices/vectors and lowercase letters for scalars. The subscript k stands for discrete time index while (•) T , (•) -1 and (•) denote the matrix transpose, matrix inverse and deconvolution, respectively. Symbol D n denotes the n-th order derivative with respect to t, i.e., D n = d n dt n and * denotes the convolution.

II. BACKGROUND

Formally, the objective of univariate smoothness priors and QV regularization strategies is to find an estimate of signal x lp (t), where the subscript lp stands for low-pass signal, from its noisy measurement y(t) in the model

y (t) = x lp (t) + v lp (t), (1)
where v lp (t) is the unwanted additive noise signal, by solving the following optimization problem

xlp (t) = argmin x lp (t) [y(τ ) -x lp (τ )] 2 dτ +α 2 d n dt n x lp (τ ) 2 dτ, (2)
where d n dt n is the n-th order derivative and α denotes a smoothness tradeoff which balances the error variance and the output smoothness. The subscript lp is used as a lowpass in the above equations since the mentioned methods act as low-pass smoothing filter [START_REF] Kheirati Roonizi | Improved smoothness priors using bilinear transform[END_REF]. In order to digitally process the signal at discrete time k, one needs to consider y[k] = y[kT s ], the discrete-time samples of y(t), where T s is the sampling period, x lp [k] the sampled desired signal and v lp [k] the sampled observation noise:

y[k] = x lp [k] + v lp [k], k = 1, • • • , L.
The smoothing procedure (2) can be then expressed as

xlp [k] = argmin x lp [k] L j=1 (y[j] -x lp [j]) 2 + α 2 L j=1 (∇ n x lp [j]) 2 , ( 3
) for all k = 1, • • • , L, where ∇ n x lp [k]
is the n-th order difference approximation of the derivative. There are several methods that can be used for the difference approximation of the derivative operator. Some of them are forward, backward, central difference rule and bilinear transform. Without loss of generality, we consider the backward/forward difference rule for the difference approximation of the derivative. In this case, the first order difference is defined by

∇x lp [k] = x lp [k] -x lp [k ± 1]
, where one with "+" holds for forward difference rule and one with "-" holds for backward difference rule. The n-th order difference is then computed as

∇ n x lp [k] = ∇(∇ n-1 )x lp [k]. The n-th order difference, ∇ n x lp [k],
can also be computed as

∇ n x lp [k] = n i=0 (-1) i n i x lp [k ± i], (4)
where n i is the binomial coefficient. In the following, we use the backward difference rule. It then means that, from now, we used x lp [k -i] in (4). Hence, (3) can be written in the following matrix notation [START_REF] Kheirati Roonizi | Improved smoothness priors using bilinear transform[END_REF]

xlp = y -x lp 2 + α 2 D n x lp 2 , (5)
where

y = (y[1], • • • , y[L]) T , x lp = (x lp [1], • • • , x lp [L]) T ,
• 2 denotes the Euclidean norm and D n is the L×L Toeplitz matrix, built from the backward difference operator d n :

d n 1, - n 1 , • • • , (-1) n-1 n n -1 , (-1) n n n .
The difference operator d n can also be computed by the following recursion:

d 1 [+1, -1] n = 1 d n = d n-1 * d 1 n > 1 , (6)
where * denotes the convolution operator. Throughout the paper, the second-order smoothness priors is of special interest. As a special case, for data on L samples, D 2 is a L×L matrix defined as:

D 2 =           1 . . . 0 0 0 0 -2 1 . . . 0 0 1 -2 1 0 . . . 0 0 1 -2 1 . . . . . . . . . . . . . . . . . . . . . 0 0 . . . 0 1 -2 1           L×L . (7)
The solution that follows from the minimization of ( 5) is

xlp = (I + α 2 D T n D n ) -1 y. Any component xlp [k]
of the vector x lp can be written in the following convolution form [START_REF] Kheirati Roonizi | Improved smoothness priors using bilinear transform[END_REF]:

xlp [k] = δ[k] + α 2 d n [-k] * d n [k] * y[k], (8)
where (•) denotes the deconvolution. ( 8) can be implemented as

xlp [k] = Z -1 1 1 + α 2 d n ( 1 z )d n (z) * y[k],
where Z -1 denotes the inverse Z-transform. The above regularization method can be stated in terms of a zero-phase, linear time-invariant (LTI) smoothing filter (for more details see [START_REF] Sameni | Online filtering using piecewise smoothness priors: Application to normal and abnormal electrocardiogram denoising[END_REF], [START_REF] Kheirati Roonizi | Improved smoothness priors using bilinear transform[END_REF]). It can be analyzed in the frequency domain. To this purpose by substituting y[k] with δ[k] in [START_REF]A smoothness priors time-varying ar coefficient modeling of nonstationary covariance time series[END_REF], and taking the Z-transform, the frequency response of the smoothing filter is obtained as follows [START_REF] Kheirati Roonizi | Improved smoothness priors using bilinear transform[END_REF]

       H lp n (z) = 1 1 + α 2 d n (z)d n ( 1 z ) = 1 1 + α 2 (1 -z -1 ) n (1 -z) n H lp n (e jω ) = 1 1 + α 2 (2 sin ω 2 ) 2n . ( 9
)
The transfer function of the smoothing filter, i.e., [START_REF] Shiller | Smoothness priors and nonlinear regression[END_REF], is real and its phase is zero. It has a low-pass characteristic which makes it suited for extracting the low frequency components of the signal, hence the subscript lp in the previous equations.

The regularization factor α determines at which frequency the transition of pass-band to stop-band occurs. By equating the Fourier frequency response, H lp n (e jω ), to 1/2, the value of α is obtained as

α o = 1 (2 sin ωcut 2 ) n ,
where ω cut is the cutoff frequency of the smoothing filter. We refer to this approach as low-pass smoothing filter (LPSF).

Note that for designing a filter, the most important parameters are the order and the cutoff frequency. So, with the optimization approach, we can directly (and easily) control the order with the number of derivatives, and the cutoff frequency is related to the hyperparameter α. The method can also be used for extracting the high frequency components: by extracting the low frequency components and subtract it from the observed signal:

xhp [k] = y[k] -xlp [k]. (10)
In this case, the transfer function of the high-pass smoothing filter is

H hp n (z) = 1 -H lp n (z), i.e.,                H hp n (z) = (1 -z -1 )(1 -z) n (2 sin ω cut 2 ) 2n + (1 -z -1 )(1 -z) n H hp n (e jω ) = (2 sin ω 2 ) 2n (2 sin ω cut 2 ) 2n + (2 sin ω 2 ) 2n ,
where the superscript hp refers to high-pass smoothing filter (HPSF). The frequency responses of an LPSF with cutoff frequency ω 1 = 0.2 and an HPSF with cutoff frequency ω 2 = 0.8 for the first and second order are depicted in Figure 1. The transition band of LPSF and HPSF increases as the order of the smoothing filter (i.e., the order of the derivative) increases. In Section III, we present a simple approach to BSSF design which is based on the combination of an LPSF in parallel with an HPSF. But, we show that the simple BSSF is less than ideal. Next, in Section IV, we show a new approach to modify that smoothing filter, making its performance nearly ideal.

H lp n (e jω ) H hp n (e jω ) + y[k] xlp [k] xhp [k] xbs [k] LPSF HPSF

III. A SIMPLE DESIGN OF BAND-STOP SMOOTHING FILTER

Let us consider the problem of signal denoising where the desired signal x bs [k] (subscript bs stands for band-stop) contains a low-frequency and a high-frequency component. We assume that the noisy signal y[k] can be represented as

y[k] = x lp [k] + x hp [k] x bs [k] +v bs [k], (11)
where v bs [k] is the unwanted additive noise. We assume that the noise v bs [k] and the signal x bs [k] are on disjoint frequency ranges. A simple approach to recover x bs [k] is to extract low and high frequency components (i.e., x lp [k] and x hp [k]) separately and combine them to extract the desired signal (i.e., x bs [k]). Therefore, a simple approach to BSSF design is to combine an LPSF in parallel with an HPSF [START_REF] Kheirati Roonizi | Improved smoothness priors using bilinear transform[END_REF].

The block diagram of a simple BSSF is given in Figure 2. Suppose that a BSSF with cutoff frequencies ω 1 and ω 2 (ω 1 < ω 2 ) is desired. By extracting the low and high frequency components using an LPSF with cutoff frequency ω 1 and an HPSF with cutoff frequency ω 2 and combining them, one can extract a desired signal with frequency components in range

of ω ∈ [0, ω 1 ] ∪ [ω 2 , 1]
. Note that LPSF and HPSF have the following transfer functions

             H lp n (e jω ) = 1 1 + ( sin ω 2 sin ω 1 2 ) 2n H hp n (e jω ) = 1 1 + ( sin ω 2 2 sin ω 2 ) 2n . (12)
Therefore, the simplest BSSF that we can create has the frequency response function:

H bs n (e jω ) = 1 1 + ( sin ω 2 sin ω 1
2

) 2n + 1 1 + ( sin ω 2 2 sin ω 2 ) 2n .
The BSSF has a center frequency at smoothing filter that shows how fast the roll-off is, is defined by 3. We find that the first order BSSF shown in Figure 3 does not have cutoff frequencies at ω 1 and ω 2 . Another drawback is that the simple BSSF does not have zero gain at center frequency (ω c = 0.4) and its vicinity. Moreover, its transition band is not narrow. The problem even becomes worse when we need to design a BSSF with high Q-factor. These types of BSSF are known as notch smoothing filters (NSFs). For example, consider a simple NSF with center frequency ω c = 0.5, ω 1 = 0.47 and ω 2 = 0.53. The magnitude response plot of this NSF is shown in Figure 4. There, we find the worst case: its cutoff frequencies are not exactly the planned-for the desired cutoff frequencies, the amplitude at the center frequency is not zero and its transition band is not enough narrow to eliminate the unwanted frequency components and pass desired frequency components. Therefore, extracting low and high frequency components separately and combining them to construct BSSF is simple but naive and weak. We propose to establish the conditions for achieving optimal performance, in extracting the low and high frequency parts simultaneously and considering the interaction between LPSF and HPSF. In the following section, we propose an improved optimization-based approach to BSSF design.

ω c = √ ω 1 ω 2 and
Q = ω c |ω 2 -ω 1 | . If Q < 1,

IV. IMPROVED DESIGN OF BAND-STOP SMOOTHING FILTER

In order to construct a BSSF, in the previous section, we naively combined an LPSF with an HPSF. However, the HPSF was designed by subtracting the output of LPSF from the observed signal (see [START_REF] Eubank | A note on smoothness priors and nonlinear regression[END_REF]). As a result, its formulation depends on the LPSF. Therefore, it is not suitable for studying the interaction between the low-pass and high-pass smoothing filters. In the following, we propose a new approach to implement the regularization method (i.e., smoothness priors or QV regularization) which makes it suited for HPSF design.

To this end, by replacing the n-th order derivative using d n following (4) and ( 6), we represent (3) in the following regularized form:

xlp [k] = argmin x lp [k] L j=1 (y[j] -x lp [j]) 2 +α 2 L j=1 (d n [j] * x lp [j]) 2 , (13)
where

d 1 [k] = δ[k] -δ[k + 1] d n [k] = d n-1 [k] * d 1 [k]
The actual minimization problem [START_REF] Brillinger | New Directions in Time Series Analysis[END_REF] is related to the state model

d n [k] * x lp [k] = e lp [k] y[k] = x lp [k] + v lp [k] , (14)
According to [START_REF] Brillinger | New Directions in Time Series Analysis[END_REF], it can be realized that the traditional smoothness priors or QV regularization is an adaptive estimation algorithm which assumes that the state of the system at a time k evolved from the prior state at previous times according to the equation [START_REF] Tarvainen | An advanced detrending method with application to hrv analysis[END_REF]. The stabilizing kernel used in the QV regularization d n [k] is chosen to be some differential operator. It imposes smoothness constraints on the solution.

Therefore, the obtained filter acts as a low-pass filter. In the following, we suggest a new state space model for the signal which the stabilizing kernel u[k], the unit step function, is chosen to be some integral operation. Using the new state space model, we propose a modified approach to smoothness priors which is suited for HPSF design.

A. A new approach to HPSF design

Let us consider the following state space model for the desired signal:

b n [k] * x hp [k] = e hp [k] y[k] = x hp [k] + v hp [k] , (15)
where

b 1 [k] = u[k] b n [k] = b n-1 [k] * b 1 [k] .
Choosing u[k] as the stabilizing kernel is equivalent to imposing the irregularity constraints on the solution. In the following, we will see that when the stabilizing kernel u[k] is used, the obtained filter has a high-pass characteristic.

According to the state space model ( 15), we propose the following regularized approach:

xhp [k] = argmin x hp [k] L j=1 (y[j] -x hp [j]) 2 +β 2 L j=1 (b n [j] * x hp [j]) 2 , (16)
where β is the regularization factor. The regularization term in ( 16) is expressed with power of the integral of the signal2 , against the regularization term in (3) which was explicitly the smoothness expressed with power of signal derivatives. ( 16) can be expressed in the following matrix notation

xhp = argmin x hp y -x hp 2 + β 2 B n x hp 2 , (17)
where B n is defined as a lower triangular matrix and b n [k] is truncated to construct the matrix filter B n . For instance, B 1 is defined by

B 1 =       1 0 . . . 0 1 1 . . . . . . . . . . . . . . . 0 1 . . . 1 1      
, which represents a cumulative sum. B 2 is defined by

B 2 =       1 0 . . . 0 2 1 . . . . . . . . . . . . . . . 0 n . . . 2 1       . ( 18
)
The solution that follows from the minimization of ( 17) is

xhp = (I + β 2 B T n B n ) -1 y.
The parameter β can be computed in terms of cutoff frequency.

To this purpose, we write the component xhp [k] of the vector x hp in the following convolution form

xhp [k] = δ[k] + β 2 b n [k] * b n [-k] * y[k]. (19)
In the Z-domain, the amplitude response of ( 19) is

H hp n (z) = 1 1 + β 2 b n (z)b n ( 1 z ) = 1 1 + β 2 [(1 -z -1 )(1 -z)] n . ( 20
)
Eq. ( 20) becomes in the Fourier domain

H hp n (e jω ) = 1 1 + β 2 2 sin ω 2 2n . (21)
By equating ( 21) to 1/2, the value of β is obtained in terms of cutoff frequency:

β o = (2 sin ω c 2 ) 2n .
Finally, the amplitude response of ( 20) is expressed as

               H hp n (z) = [(1 -z -1 )(1 -z)] n (2 sin ω c 2 ) 2n + [(1 -z -1 )(1 -z)] n H hp n (e jω ) = (2 sin ω 2 ) 2n (2 sin ω c 2 ) 2n + (2 sin ω 2 ) 2n . (22)
Comparing ( 22) and ( 12), we conclude that the frequency response of the new regularization model is equivalent to the HPSF described in the previous section. It acts as an HPSF.

The advantage of the new HPSF is that its model does not depend on the LPSF model. Therefore, it can be used for studying the interaction between the LPSF and HPSF when we combine them to construct a BSSF. In the following section, the new HPSF is combined with the LPSF for designing the BSSF.

B. A new model to band-stop smoothing filter design

In order to overcome the drawbacks of the simple parallel approach discussed in Section III, in the following, we consider the interactions between the LPSF and HPSF. The proposed method is as follows:

We are interested in the estimation of the signal x bs [k] from its noisy measurements y[k] in the model [START_REF] Kitagawa | Smoothness prior approach to explore mean structure in large-scale time series[END_REF]. Figure 5 depicts our proposed scheme for BSSF, and the following equations give its mathematical representation

         d n [k] * x lp [k] = e lp [k] b n [k] * x hp [k] = e hp [k] y[k] = x lp [k] + x hp [k] x bs [k] +v bs [k]
This model assumes that the lp and hp part of the signal of interest are modelled by filtering a white noise. We assume

1 (1-z -1 ) n (1 -z -1 ) n + e lp [k] e hp [k] x lp [k] x hp [k] x bs [k] + v bs [k] y[k] H bs n (z) xbs [k]
E e lp [k] e hp [k] e * lp [m] e * hp [m] = σ 2 lp ρσ lp σ hp ρσ lp σ hp σ 2 hp δ k,m , (23)
where E{} refers to the expected value operation, σ 2 lp and σ 2 hp denote the noises variance and ρ is the noise correlation.

Considering the fact that the system noises can be correlated, we propose to estimate the desired signal xbs [k] by solving the following least-squares optimization problem

xlp [k], xhp [k] = argmin x lp [j],x hp [j] L j=1 (y[j] -x lp [j] -x hp [j]) 2 + α 2 L j=1 (d n [j] x lp [j]) 2 + β 2 L j=1 (b n [j] x hp [j]) 2 + 2ραβ L j=1 (d n [j] x lp [j]) (b n [j] x hp [j]) . (24)
The first regularization term is used to control the smoothness of x lp , the second regularization term controls the sharpness of x hp and the third regularization term is used to control the interaction between the smoothness of x lp and the sharpness of x hp .

Equation ( 24) can be expressed in the following matrix notation

xlp , xhp = argmin x lp ,x hp y -x lp -x hp 2 + α 2 D n x lp 2 + β 2 B n x hp 2 + 2ραβ (D n x lp ) T (B n x hp ) , (25)
where (•) T denotes the matrix transpose operator. Note that

D T n B n = B T n D n = I.
Hence, (25) can be expressed as [START_REF] Woodford | Global stereo reconstruction under second-order smoothness priors[END_REF] xlp , xhp = argmin

x lp ,x hp y -x lp -x hp 2 + α 2 D n x lp 2 + β 2 B n x hp 2 + 2ραβx T lp x hp .
By setting the derivatives of ( 26) with respect to x lp and x hp equal to zero, after some simplifications, we get the following system of equations for x lp and x hp :

M 11 x lp + M 12 x hp = y M 21 x lp + M 22 x hp = y , (27)
where

     M 11 = I + α 2 D T n D n M 12 = M 21 = (1 + ραβ) I M 22 = I + β 2 B T n B n Solving the set of equations (27), we find xlp = [M 11 M 22 -M 12 M 21 ] -1 [M 22 -M 12 ] y xhp = [M 11 M 22 -M 12 M 21 ] -1 [M 11 -M 21 ] y
Finally, the desired signal, xbs , is obtained as ( 28)

xbs = xlp + xhp = [M 11 M 22 -M 12 M 21 ] -1 [M 11 + M 22 -(M 12 + M 21 )] y.
The only unknown parameters are the regularization factors α and β. In the following section, the design parameters are calculated in terms of the BSSF cutoff frequencies.

C. Frequency domain analysis: parameter design based on cutoff frequencies

The BSSF developed in [START_REF] Hodrick | Postwar u.s. business cycles: An empirical investigation[END_REF], was derived from a regularized least-squares optimization in form of block-wise matrix formulation. It can also be stated in terms of a zero-phase, forward filtering backward smoothing scheme. Therefore, the regularization factors α and β can be calculated in terms of the BSSF cutoff frequencies. Any component xlp [k], xhp [k] and xbs [k] of the vector x lp , x hp and x bs can be written in the convolution form represented in [START_REF] Savitzky | Smoothing and differentiation of data by simplified least squares procedures[END_REF], on the top of the next page, where

     m 11 [k] = δ[k] + α 2 d n [-k] * d n [k] m 12 [k] = m 21 [k] = (1 + ραβ)δ[k] m 22 [k] = δ[k] + β 2 b n [-k] * b n [k]
The frequency response of ( 29) is

H bs n (z) = M 11 (z) + M 22 (z) -[M 12 (z) + M 21 (z)] M 11 (z)M 22 (z) -M 12 (z)M 21 (z) , (30)
where

M ij (z) is the Z-transform of m ij [k]:          M 11 (z) = 1 + α 2 d n ( 1 z )d n (z) M 12 (z) = M 21 (z) = 1 + ραβ M 22 (z) = 1 + β 2 b n ( 1 z )b n (z)
Equation ( 30) in Fourier domain becomes

H bs n (e jω ) = κ (αβ) 2 (1 -ρ 2 ) + κ , (31)
       xlp [k] = (m 11 [k] * m 22 [k] -m 12 [k] * m 21 [k]) * (m 22 [k] -m 21 [k]) * y[k] xhp [k] = (m 11 [k] * m 22 [k] -m 12 [k] * m 21 [k]) * (m 11 [k] -m 12 [k]) * y[k] xbs [k] = (m 11 [k] * m 22 [k] -m 12 [k] * m 21 [k]) * (m 11 [k] + m 22 [k] -m 21 [k] -m 12 [k]) * y[k] (29)
where

κ = α 2 (2 sin ω 2 ) 2n -2ραβ + β 2 (2 sin ω 2 ) 2n
.

The magnitude response of the BSSF, [START_REF] Boyd | Convex Optimization[END_REF], has two cutoff frequencies at ω 1 and ω 2 which are satisfying |H bs n (e jω )|= 1 2 . Substituting each of these two points into equation (31) results a system of equations which contains two equations that share two unknowns α and β:

α 2 (2 sin ω 1 2 ) 2n -2ραβ + β 2 (2 sin ω1 2 ) 2n = α 2 β 2 (1 -ρ 2 ) α 2 (2 sin ω 2 2 ) 2n -2ραβ + β 2 (2 sin ω2 2 ) 2n = α 2 β 2 (1 -ρ 2 ) ( 32
)
Subtracting the second equation from the first one, we find

β = α (2 sin ω 1 2 )(2 sin ω 2 2 ) n (33)
Substituting [START_REF] Engl | Discrepancy principles for tikhonov regularization of ill-posed problems leading to optimal convergence rates[END_REF] in the first equation of ( 32), the resulting equation is

(34) α 2 2 2n [(sin ω 1 2 ) 2n -2ρ(sin ω 1 2 sin ω 2 2 ) n + (sin ω 2 2 ) 2n ] = α 2 β 2 (1 -ρ 2 )
After some simplifications, we find 33) and [START_REF] Eldar | Generalized sure for exponential families: Applications to regularization[END_REF], the value of α is obtained as

β o = 2 n (sin ω1 2 ) 2n -2ρ(sin ω1 2 sin ω2 2 ) n + (sin ω2 2 ) 2n 1 -ρ 2 (35) From (
α o = β o (2 sin ω1 2 )(2 sin ω2 2 ) n ( 36
)
It means that by setting the value of α and β respectively to ( 36) and ( 35), the proposed optimization approach, i.e., [START_REF] Hodrick | Postwar u.s. business cycles: An empirical investigation[END_REF], acts as a BSSF with cutoff frequencies ω 1 and ω 2 . As an example, the amplitude response of the new BSSF for ω 1 = 0.2 and ω 2 = 0.8, n = 1, 2 and different values of ρ ∈ (-1, 1) is shown in Figure 6. It shows that the transition band of the BSSF increases as ρ becomes close to 1. We find that at ρ = 1 -≈ 1 ( is a small value), the BSSF has narrower transition band and its gain at ω = ω c is zero. Therefore, by setting the value of ρ close to 1, a nearideal BSSF is obtained. In Figure 6, we set to 10 -3 . The amplitude response of the BSSF, for ρ = 0.999, ω 1 = 0.2, ω 2 = 0.8 and n = 1, 2 is shown in Figure 7. Repeating our last example, for second order NSF with ω c = 0.5, ω 1 = 0.47, ω 2 = 0.53 and different values of ρ, we find that the frequency response of the improved NSF, shown in Figure 8 is nearly ideal when the value of ρ is sufficiently close to 1. The results of Figure 8 also show that the transition band of the NSF is related to ρ. In [START_REF] Roonizi | A New Approach to ARMAX Signals Smoothing: Application to Variable-Q ARMA Filter Design[END_REF], a general approach for ARMA with exogenous inputs (ARMAX) signals smoothing with application to variable-Q ARMA smoothing filter design was proposed. Specially, a narrowband notch smoothing filter with variable-Q was designed. The NSF proposed in [START_REF] Roonizi | A New Approach to ARMAX Signals Smoothing: Application to Variable-Q ARMA Filter Design[END_REF] employs an ARMA model. In Figure 9, we plot the amplitude response of the NSF proposed in [START_REF] Roonizi | A New Approach to ARMAX Signals Smoothing: Application to Variable-Q ARMA Filter Design[END_REF] and the first order NSF in the present paper with cutoff frequencies ω 1 = 0.07 and ω 2 = 0.13. We can see from the amplitude responses that the proposed NSF acts better in the pass band. The NSF proposed in [START_REF] Roonizi | A New Approach to ARMAX Signals Smoothing: Application to Variable-Q ARMA Filter Design[END_REF] does not have unit gain in the pass band ω ∈ [0, ω 1 ] while the proposed NSF has approximately unit gain in both pass bands, ω ∈ [0, ω 1 ] and ω ∈ [ω 2 , π]. We also plot the amplitude response of the first order Butterworth filter and Kaiser with order 40. Kaiser method has better amplitude response in its pass-band and stop-band, but compared to the proposed BSSF, its order is dramatically higher. It is notable that the design of notch filters is a classic problem in recursive (IIR) filter design. Especially, a second-order transfer system:

H n (z) = 1 -2 cos(ω 0 )z -1 + z -2 1 -2ρ cos(ω 0 )z -1 + ρ 2 z -2 (37)
with a pair of zeros on the unit circle and a pair of poles place inside the unit circle makes a simple and classical notch filter [START_REF] Pei | Narrowband notch filter using feedback structure[END_REF], [START_REF] Tseng | Stable IIR notch filter design with optimal pole placement[END_REF]. However, the classical notch filter has several drawbacks that has been discussed in [START_REF] Roonizi | A New Approach to ARMAX Signals Smoothing: Application to Variable-Q ARMA Filter Design[END_REF]. For avoiding duplication of results presented in [START_REF] Roonizi | A New Approach to ARMAX Signals Smoothing: Application to Variable-Q ARMA Filter Design[END_REF], we do not give theoretical derivations or algorithm descriptions that are easily found in [START_REF] Roonizi | A New Approach to ARMAX Signals Smoothing: Application to Variable-Q ARMA Filter Design[END_REF]. It has been shown that the approach presented in [START_REF] Roonizi | A New Approach to ARMAX Signals Smoothing: Application to Variable-Q ARMA Filter Design[END_REF] is more effective than the classical notch filter described in [START_REF] Pei | Narrowband notch filter using feedback structure[END_REF], [START_REF] Tseng | Stable IIR notch filter design with optimal pole placement[END_REF]. That is why we do not compare our proposed filter with filters designed using more classical methods.

Finally, a BPSF with cutoff frequencies ω 1 and ω 2 can be obtained by subtracting the BSSF output (28) from the observed signal (i.e., xbp = yxbs ).

D. A Causal filter

The designed band-stop smoothing filter is a non-causal filter as it takes all the measurements into account [see [START_REF] Grimson | From Images to Surfaces: A Computational Study of the Human Early Visual System[END_REF] and ( 25)]. So it can only be used for off-line applications. In order to implement it for online applications, in the following, we propose a causal filter.

From [START_REF] Hodrick | Postwar u.s. business cycles: An empirical investigation[END_REF], it can be seen that the matrix

Υ = [M 11 M 22 -M 12 M 21 ] -1 [M 11 + M 22 -(M 12 + M 21 )
] is a symmetric positive-definite matrix. Therefore, it has a unique Cholesky decomposition of the form Υ = L T L, where L is a lower triangular matrix with real and positive diagonal entries. Letting J denote the exchange matrix (square matrix with 1's on the anti-diagonal), L T is equal to JLJ. Therefore, (28) can be expressed as xbs = JLJLy,

As a result, we can implement it using forward-backward filtering scheme with an FIR filter (x f bs = Ly, where the superscript f stands for forward filter). For more details, see [START_REF] Roonizi | An efficient algorithm for manuvering target tracking[END_REF].

V. APPLICATION TO ECG DENOISING

In order to evaluate the performance of the proposed method, we employed it for ECG signal denoising. To this purpose, we applied the approach on ECG data from Phys-ioNet. The dataset contains 80 records, originally provided for the PhysioNet/Computers in Cardiology Challenge 2004 [START_REF] Moody | Spontaneous termination of atrial fibrillation: A challenge from physionet and computers in cardiology 2004[END_REF]. Each record, extracted from a two-lead, sampling frequency 128 Hz Holter ECG recording, 1 minute in length. In many . The amplitude response of the NSF proposed in [START_REF] Roonizi | A New Approach to ARMAX Signals Smoothing: Application to Variable-Q ARMA Filter Design[END_REF] and the first order NSF proposed in the present paper with cutoff frequencies ω 1 = 0.07 and ω 2 = 0.13. The proposed NSF has approximately unit gain in either pass band. We also plot the amplitude response of the first order Butterworth filter and Kaiser with order 40.

applications of signal processing such as ECG denoising, linear time invariant (LTI) filters (e.g., Butterworth filter) are still the standard choice for ECG system front-ends [START_REF] Widmann | Digital filter design for electrophysiological data -a practical approach[END_REF]. Therefore, we also employed a band-stop Butterworth filter for extracting the ECG signals and compared it with our proposed method. In the following case studies, the BSSF order is set to n = 2, and the operators D 2 and B 2 are as defined in [START_REF] Kitagawa | A smoothness priors-state space modeling of time series with trend and seasonality[END_REF] and [START_REF] Mao | Spatio-temporal signal recovery based on low rank and differential smoothness[END_REF], respectively. The Butterworth filters were applied to the data using the filtfilt function in Matlab to have zero-phase lag.

As a preliminary example, we plot a specific case (record b01m from PhysioNet/Computers in Cardiology Challenge 2004 [START_REF] Moody | Spontaneous termination of atrial fibrillation: A challenge from physionet and computers in cardiology 2004[END_REF]) at the top of Figure 10(a). At the bottom of that figure, we show its spectral representation.

At the top of Figure 10(b), we plot the same ECG, but corrupted with a synthetic band-limited noise (SN R = -15 dB), generated by the following equation:

v bs (t) = N -1 i=0 c i cos(2πf i t + ψ i ) (38)
where c i , f i and ψ i are respectively, the amplitude, frequency and phase of the i-th sinusoid. In this example, the frequencies of the noise were selected in range of f i ∈ (48, 52) Hz in [START_REF] Kheirati Roonizi | Forward-backward filtering and penalized least-squares optimization: A unified framework[END_REF]. Therefore, a BSSF with cutoff frequencies 48 and 52 Hz could be used to remove the noise and extract the original ECG. The results of applying the second order simple BSSF (the naive BSSF presented in Section III), the second order zero-phase Butterworth filter and the proposed BSSF (n = 2) on the noisy ECG as well as their spectral representation are reported in Figures 10(c)-10(e). The denoised ECG using the proposed BSSF is visually closer to the true ECG, compared to the denoised ECG by simple BSSF and zero-phase Butterworth filter. The spectral representation of the signals are also reported which confirms the superiority of the proposed BSSF visually (see spectrograms in Figure 10). For more details, we also report the values of NSR defined in [START_REF] Unser | Recursive regularization filters: design, properties, and applications[END_REF] obtained by each method. The NSR obtained by our proposed method is smaller than other two methods which indicates the superiority of our proposed BSSF.

For evaluating quantitatively the performance of the different BSSFs, we tested them with variable SNR in the range -20 dB to 20 dB, by controlling the power of additive noise v bs . The frequency range of the additive noise were selected such that f i ∈ [46 ± 4, 54 ± 4] Hz in [START_REF] Kheirati Roonizi | Forward-backward filtering and penalized least-squares optimization: A unified framework[END_REF].

The average and standard deviation of the SNR improvement obtained after filtering, is denoted SN R dif f = SN R out -SN R in and is used as a measure of performance. It is shown in Figure 11(a). The results, reported in Figure 11(a), show that the proposed BSSF outperforms the simple BSSF and zero-phase Butterworth filter. In Figure 11, we also report the result of the third order zero-phase Butterworth filter. In this case, also our proposed BSSF (with n = 2) outperforms the third order Butterworth filter. Other metrics that was used for evaluating the performance of the proposed method, are a classical ratio between the power of the reconstruction error and the power of the original signal (a noise-to-signal ratio denoted NSR), given by [START_REF] Kheirati Roonizi | Improved smoothness priors using bilinear transform[END_REF], [START_REF] Kheirati Roonizi | A Signal Decomposition Model-Based Bayesian Framework for ECG Components Separation[END_REF]

NSR = k (x[k] -x[k]) 2 k x 2 [k] , (39)
and the ratio between the power of the reconstruction error and the power of noise [START_REF] Kheirati Roonizi | Improved smoothness priors using bilinear transform[END_REF], [START_REF] Kheirati Roonizi | A Signal Decomposition Model-Based Bayesian Framework for ECG Components Separation[END_REF] imp = -10

log 10 k (x k -x k ) 2 k (y k -x k ) 2 (dB). ( 40
)
The results of the reconstruction procedures using these metrics are reported in Figures 11(b) and 11(c). These metrics also show that the proposed BSSF outperformed simple BSSF and zero-phase Butterworth filter. Figure 11(b) shows that the error of reconstruction for all methods decreases as the SNR increases (i.e., the reconstruction improvement continues for all methods). The NSR obtained by the proposed BSSF is smaller than the NSR obtained by Butterworth BSSF and simple BSSF which indicates the superiority of the proposed method. Furthermore, Figure 11(c) shows that the SNR improvement obtained by the proposed BSSF is larger than other two methods for all SNR, which means the proposed method outperforms Butterworth and simple BSSF in signal reconstruction.

VI. DISCUSSION AND CONCLUSION

A. Discussion

In this paper, we used the backward difference rule for implementing the derivative and integral due to its simplicity. Actually, the implementation of the unit sequence u[k] is more simple and straightforward than trapezoidal integrator and other methods. In fact, we have already employed the Tustin's or trapezoidal method for improving the smoothness priors in our recent paper [START_REF] Kheirati Roonizi | Forward-backward filtering and penalized least-squares optimization: A unified framework[END_REF]. In [START_REF] Kheirati Roonizi | Forward-backward filtering and penalized least-squares optimization: A unified framework[END_REF], we have shown that implementing the derivative operator using bilinear transform (Tustin's method) improves the smoothness priors method. In the current paper, we mainly focused on proposing a new idea to band-stop smoothing filter design. However, we have included the Tustin's method and compared it with backward difference rule. Figure 12 shows the amplitude response of the proposed BSSF using Tustin and backward difference rule with cutoff frequencies ω 1 = 0.47 and ω 2 = 0.53 (top left ρ = 0.5, top right panel ρ = 0.75, bottom left ρ = 0.9 and bottom right ρ = 0.999). The results show that for small values of ρ, the proposed BSSF using Tustin's method has sharper transition band while for large values of ρ they are close to each other. It means that for higher values of ρ that the BSSF has the sharpest frequency transition band, both methods are close to each other. That is why we chose the backward difference rule in our paper to implement the BSSF. The first and second order derivatives are two main basic operations and blocks in signal filtering/smoothing. In fact, the integer order differential and integral operators are defined uniquely. They are local as they consider the values of the neighboring points to the point of interest. The problem formulation and algorithms described in this paper be implemented using other methods such as fractional order derivative/integral. The fractional calculus operators are non-local and larger neighborhoods must be considered in the computation. Thus, it results in long-term memory effect. In the literature, a plenty of definitions of the fractional order derivatives [START_REF] Ortigueira | An introduction to the fractional continuous-time linear systems: the 21st century systems[END_REF] can be found: forward Grunwald-Letnikov (GL), Riemann-Liouville (RL) fractional derivatives, Fourier domain fractional derivatives.

In a preliminary experiment, the proposed BSSF has been implemented using fractional calculus and the Fourier domain fractional derivatives was used to implement the frequency response of the proposed BSSF. To demonstrate the impact of fractional order on the shape of the smoothing filter, we plot the amplitude response of the BSSF at various values of n in Figure 13. Figure 13 shows that for lower values of n, the BSSF has a shallow frequency transition band. We will consider further investigations in future work.

In this paper, we compared the proposed BSSF with maximally-flat Butterworth BSSF. The reason for choosing the maximally flat Butterworth filter for comparison is that our proposed least squares optimization based BSSF design has also a maximally flat frequency response. But the difference between them is that our proposed BSSF is directly designed in time domain while the zero-phase Butterworth is first designed in frequency domain, following by polynomial root finding, and finally applying forward filtering backward smoothing which could be quite involved. Another method that has been used for digital filter design is based on optimization technique method [START_REF]Design of digital IIR filter: A research survey[END_REF]. Numerous optimization methods have been developed based on modern heuristics optimization algorithms: Simulated Annealing (SA) algorithm, genetic algorithm (GA), particle swarm optimization (PSO), artificial bee colony and etc.. However, for optimal filter, designed in the frequency domain by the optimization methods, the error function and hence the optimization problem are generally nonconvex with multiple local optimal points. As a result, while it demands a large amount of computations, the designed filters need to be converted into time domain [START_REF]Design of digital IIR filter: A research survey[END_REF]. In the literature, the deep learning (DL) method has been also used for signal denoising, e.g. for ECG denoising [START_REF] Arsene | Deep learning models for denoising ecg signals[END_REF]- [START_REF] Chiang | Noise reduction in ecg signals using fully convolutional denoising autoencoders[END_REF]. However, in these papers, even if results seem interesting at first glance, comparison with the best more classical methods are still missing. Moreover, one of the advantages of the method we proposed in this paper is to obtain a closed-form equation relating each cutoff frequency to one parameter used in the optimization algorithm: At first glance, such a cut-off frequency control seems impossible in DL. Anyway, recent theoretical works in explainable DL introduce the idea of algorithm unrolling [START_REF] Monga | Algorithm unrolling: Interpretable, efficient deep learning for signal and image processing[END_REF] which relates each iteration of an optimization algorithm to one layer in a deep neural networks. In further work, following this idea, it seems interesting to investigate how algorithms we proposed could be implemented in deep neural networks and compare the performance (accuracy, computational load, etc.) with respect to more classical implementations. Finally, the integration of the proposed optimization approach into advanced control methods, such as those in [START_REF] Sun | Fuzzy adaptive finite-time faulttolerant control for strict-feedback nonlinear systems[END_REF]- [START_REF] Sun | A novel finite-time control for nonstrict feedback saturated nonlinear systems with tracking error constraint[END_REF] and its sensitivity on the parameters and faults will be investigated in the future.

B. Conclusion

Quadratic variation (QV) regularization and smoothness priors are the most commonly used signal denoising techniques. In this paper, an extension of such algorithms to band-stop smoothing filters was investigated. We showed that with the optimization approaches (QV regularization or smoothness priors), the cutoff frequencies are related to the regularized parameters and the order can be directly (and easily) controlled with the number of derivatives. First, a simple (naive) approach was presented. The simple band-stop smoothing filter can be made out of a low-pass smoothing filter and a highpass smoothing filter by connecting the two smoothing filter sections in parallel with each other. The simple band-stop smoothing filter suffers from some shortcomings. Its cutoff frequencies are not exactly the planned-for the desired cutoff frequencies, the amplitude at the center frequency is not exactly zero and its transition band is not enough narrow to eliminate the unwanted frequency components and pass desired frequency components. In order to overcome these problems and further improve the naive approach, we presented a new approach to band-stop smoothing filter design.

In the proposed approach, we considered a modified model and investigated the role of noise correlation which represents the link between two smoothing filter sections. Our results suggest that the noise correlation increases the steepness of the role-off and boosts the performance. Moreover, positive noise correlation is the key to always get the best bandstop smoothing filter. Specifically, ρ ≈ 1, gives band-stop smoothing filter with the steepest frequency transition band. The relationship between the cutoff frequencies of the filter and the parameters used in the optimization algorithm was shown in the paper. More precisely, choosing the desired cutoff frequencies lead to choosing the hyperparameters, i.e., the weights of the regularization factors.

Fig. 1 .

 1 Fig. 1. The amplitude response of LPSF with cutoff frequency ω 1 = 0.2 (left panel) and HPSF with cutoff frequency ω 2 = 0.8 (right panel) for different values of n.

Fig. 2 .

 2 Fig. 2. The simple parallel structure to band-stop smoothing filter (BSSF) design.

Fig. 4 .

 4 Fig. 4. The amplitude response of a simple NSF with cutoff frequencies ω 1 = 0.47 and ω 2 = 0.53 for different values of n.

Fig. 5 .

 5 Fig. 5. The proposed scheme for band-stop smoothing filter design.

Fig. 6 .Fig. 7 .

 67 Fig.6. The amplitude response of the proposed BSSF with cutoff frequencies ω 1 = 0.2 and ω 2 = 0.8, for different values of ρ a) the first order b) the second order. In this experiment, we set to 10 -3 .

Fig. 8 .

 8 Fig. 8. The amplitude response of the second order NSF with cutoff frequencies ω 1 = 0.47 and ω 2 = 0.53, for different values of ρ. The transition band of the proposed NSF increases as the value of ρ increases.

Fig. 9

 9 Fig.9. The amplitude response of the NSF proposed in[START_REF] Roonizi | A New Approach to ARMAX Signals Smoothing: Application to Variable-Q ARMA Filter Design[END_REF] and the first order NSF proposed in the present paper with cutoff frequencies ω 1 = 0.07 and ω 2 = 0.13. The proposed NSF has approximately unit gain in either pass band. We also plot the amplitude response of the first order Butterworth filter and Kaiser with order 40.

  (a) ECG record.

Fig. 10 .Fig. 11 .

 1011 Fig. 10. ECG denoising for record b01m from PhysioNet/Computers in Cardiology Challenge 2004. a) ECG record b) noisy ECG (SNR = -15 dB) c) provided by second order simple BSSF designed in Section III (NSR = 2.1) d) second order zero-phase band-stop Butterworth filter (NSR = 0.1) e) proposed BSSF for n = 2 (NSR = 0.02). The spectral representation of the signals is also reported.

Fig. 12 .

 12 Fig.[START_REF] Kitagawa | The Smoothness Priors Concept[END_REF]. The comparison of the amplitude response of the proposed BSSF using Tustin and backward difference rule with cutoff frequencies ω 1 = 0.47 and ω 2 = 0.53 (top left ρ = 0.5, top right panel ρ = .75, bottom left ρ = 0.9 and bottom right ρ = 0.999). For small values of ρ, the proposed BSSF using Tustin method has sharper transition band while for large values of ρ they are close to each other.

Fig. 13 .

 13 Fig.[START_REF] Brillinger | New Directions in Time Series Analysis[END_REF]. The comparison of the amplitude response of the proposed BSSF using fractional calculus for ρ = 0.99. For large values of n, the proposed BSSF has sharper transition band.

The quality factor or Q-factor of a filter is defined as the ratio of its center frequency over its bandwidth.

In continuous time, we haveu(t) * x(t) = t -∞ x(τ )dτ and in discrete time, we have u[k] * x[k] = k j=-∞ x[j].

ACCEPTED PAPER FOR PUBLICATION IN IEEE TRANSACTIONS ON SIGNAL PROCESSING, DOI 10.1109/TSP.2021.3060619