Signal denoising based on the Schr\"odinger operator's eigenspectrum and a curvature constraint

Recently, a new Signal processing method, named Semi-Classical Signal Analysis (SCSA), has been proposed for denoising Magnetic Resonance Spectroscopy (MRS) signals. It is based on the Schr\"odinger Operator's eigenspectrum. It allows an efficient noise reduction while preserving MRS signal's peaks. In this paper, we propose to extend this approach to different signals, in particular pulse shaped signals, by including an optimization that considers curvature constraints. The performance of the method is measured by analyzing noisy signal data and comparing with other denoising methods. Results indicate that the proposed method not only produces good denoising performance but also guarantees the peaks are well preserved in the denoising process.


Introduction
The Semi-Classical Signal Analysis (SCSA) method decomposes the signal into a set of functions given by the squared eigenfunctions of the Schrödinger operator associated to its negative eigenvalues (2).Thus, and unlike traditional signal decomposition tools, the SCSA expresses the signal through a set of functions that are signal dependent, i.e these functions are not fixed and known in advance but are computed by solving the spectral problem of the Schrödinger operator whose potential is the signal to be analyzed.These eigenfunctions will capture more details about the signal and its morphological variations.
The SCSA has been successfully applied in many applications for signal representation, denoising, post-processing and feature extraction.For example, it has been used for arterial blood pressure waveform analysis in (2), ( 3), (4) and for Magnetic Resonance Spectroscopy (MRS) denoising (1) and MRS water suppression (5).Is has been also used for feature extraction in epileptic seizure detection using Magnetoencephalography (MEG) signals (6).
The SCSA has shown very good performance in analyzing pulse-shaped signals that can be found in many applications, in particular, in biomedical applications.For this type of signals, the peak shape as well as the peak position are of paramount importance (7).
However, data come with noise and error, with a variety of noise origins such as electronic, mechanical and optical interferences, causing signal spectrum to be noisy regardless of how careful the experiment is carried out.In addition, the sensitivity of the instruments tends to drop down with use and the signals tend to have strong interferences from the background noises.
Conventional methods for removing noise exist, such as improving accuracy by finding the source of noise and eliminating its effect at the data acquisition stage, or suppressing the noise by replicating the measurements.It can be seen that both of the approaches are practically not feasible.The first approach to find noise in a highly sophisticated instrument has hidden requirements of a significant amount of expertise in that field, while the second is not feasible for financial considerations if samples are of biological, clinical or pharmaceutical origins (4).Therefore, signal processing methods are always needed in such scenarios to denoise pulse-shaped signals.
In this paper, we propose to extend the SCSA method to a more general signal denoising framework and analyze the performance of this approach.To deal with more general pulse shaped signals, we propose to combine the SCSA with an optimization that computes the optimal values for the semi-classical parameter under some curvature constraints.However and unlike the previous contribution (1), which uses prior knowledge on the position and the peaks of the MRS spectrum and also the predominance of the noise in some specific areas allowing the computation of the SNR, in this paper we use a curvature constraint that makes the approach applicable for any pulse shaped signal.We refer to the algorithm introduced in (1) (α-SCSA) and the new proposed algorithm C-SCSA for Curvature-SCSA.
The paper is structured as follows.Section II provides a brief introduction of the SCSA method, and a description of the notation required in the rest of the paper.In Section III the major concepts of curvature denoising as well as conventional α-SCSA denoising are described.Consequently, a novel SCSA-based denoising strategy is presented along with other popular denoising techniques.The performance evaluation of the novel denoising technique C-SCSA is illustrated in Section IV.The performance of the different SCSA techniques and existing popular denoising methods are illustrated in Section IV, and the final conclusions are drawn in Section V.

SCSA for signal reconstruction
The SCSA method decomposes a real positive signal y(t) into a set of squared eigenfunctions through the discrete spectrum of the Schrödinger operator.The reconstructed signal y h (t) is presented by the squared eigenfunctions (refer to (2)) where λ nh = −κ 2 nh , with κ 1h > κ 2h > • • • > κ nh are the negative eigenvalues, and The L 2 2 -normalized eigenfunctions reconstructed from a pulse shaped signal are shown in Fig. 1.N h is the number of negative eigenfunctions and h is a positive parameter known as the semi-classical constant.It is found that when h tends to zero, the reconstructed spectrum y h converges to the true spectrum y.This matches the semi-classical properties of the Schrödinger operator where the number of negative eigenvalues thus the number of corresponding eigenfunctions increases when h decreases.One of the important characteristics is that eigenfunctions which correspond to large eigenvalues represent the profiles of the peaks, whereas the remaining functions characterize the noise details of these profiles.2 Find the negative eigenvalues and the associated eigenfunctions for the matrix above.
3 Normalize the eigenfunctions and compute y h using Equation (1).
4 If the reconstruction is accurate, stop; if not, decrease the value of h and go to step

SCSA for signal denoising
Now, let's consider the following noisy signal where y(t) is the noiseless signal and n(t) is the additive noise.The aim of digital signal denoising is to produce an accurate estimate of the original signal y(t).The noise variance of n(t), which can be known or unknown depending on different cases, is denoted as σ As in the reconstruction, the parameter h plays a key role in the denoising with the SCSA method.On one hand, when h tends to zero, the reconstructed spectrum y h converges to the noisy signal y δ .On the other hand, It is demonstrated that N h increases when h decreases and the squared eigenfunctions ψ nh are such that the number of oscillations in the eigenfunctions increases with the order n while their amplitude decreases.Therefore, the highest order eigenfunctions will mainly reconstruct the noise.
In a general sense, as N h value increases, both the original signal and noise will be gradually reconstructed, first the major peaks and details of the signal and then the noise, as shown in Fig. 2. A selection of optimal h has to be made in order to separate noise from the original signal.In simulation, h is initiated at a relatively small value at first and then gradually increased to discard the noise part.One can infer that the choice of the stop criterion is critical, since it sets the optimum h value, which leads to a reliable signal reconstruction and therefore to an accurate data analysis.It is found that the best stop criterion is a function not only of the minimum reachable distance between y δ (t) and reconstructed signal y h , but also of the Signal to Noise Ratio (SNR) of y h reconstructed spectrum (1): where M represents the number of peaks and α is a positive weight parameter, which allows to balance accuracy and denoising.SNR y h is the SNR of y h computed using the following formula where the interval [t 1 , t 2 ] is the interval where the noise is dominant, max and std represent the maximum and standard deviation of the function y h This type of cost function is quite standard, and keeps a balance between fidelity to the signal and denoising effect.However, this method has some pre-assumptions, in which signal peak localization is needed and also the interval where the noise is dominant needs to be known.

Curvature based SCSA Denoising
In this section we propose a new denoising algorithm based on the SCSA and some curvature  general principles of denoising is still the same, i.e. to come up with a standard of properly select h in order to reduce the noise of the signal, the cost function selection will need more mathematical intuition.
Therefore, inspired by the smoothing methods of ( 8), we propose a cost function J in the following form: where |k(t)| is a certain smoothness penalty term which operates on the reconstructed signal y h .|k(t)| dt describes the "wiggliness" of the reconstructed signal y h , µ is a non-negative smoothing parameter that needs to be properly selected.It depends on the characteristic of the input signal.Larger values of µ force y h to be smoother.In this paper we define smoothness penalty term k(t) to be the curvature.Without loss of generality, let's assume x m and w m to be jointly Gaussian and zero mean with variance σ 2 m , with their joint distribution defined as: where ρ m = COV {xm,wm} . If the signal are equally spaced with interval denoted as ∆, the curvature defined in Eq. ( 6) can be approximated as Eq. ( 8) at the mth sample : Proposition 1 (Expectation of the curvature) Given f (x m , w m ) and k m in Eq. ( 7) and (8), E{k m } can be approximated as: The proof of the theorem is given in the appendix.Since x m and w m are subtraction between two neighborhood samples, they resemble noise distributions in homogeneous regions of the signal (a common assumption for noise propagation analysis like ( 9)), whose distribution σ 2 m are of the similar distribution of the noise, i.e. σ 2 m 2σ 2 noise when noise is identically distributed.Therefore when noise level σ noise increases, the σ 2 m also increases.Given Eq. ( 9) one can easily infer that when σ 2 m increases, E{k m } also increases.Therefore, in order to reduce noise level, we use curvature term k in cost function J.Then, given N samples of the signal y δ , we propose a scanning method which iteratively scans h, to minimize the cost where y δ is the input noisy signal and y h is the signal reconstructed with SCSA method.N is the number of samples in y δ .µ is automatically adjusted given a specific type of signal ) in order to make sure that fidelity term and penalty term are of similar order of magnitude.The algorithm makes a step further comparing to (1) in that it doesn't need to locate the signal peaks during the denoising process, which in some cases is not trivial.The fidelity term in the cost function J ranges the whole input noisy signal.Also, the SNR is unknown in many cases, while curvature is a characteristic of all types of signals that can be numerically computed.
4 Simulation Results

Denoising methods
There are many existing methods for pulse-shaped signal denoising.With respect to traditional digital filters, moving average filter is the simplest method (10,11), with the drawbacks of temporal autocorrelation at a lag determined by the length of the moving window, thus shifting the peak position which will negatively affect the signal.In comparison, the Savitzky-Golay algorithm (SG) is more effective in smoothing noisy data obtained from spectrum data (for example (12)) and is currently the most commonly applied filter to eliminate the irrelevant information from noisy input data (13).Other well-known candidates for denoising are techniques based on wavelets (16).There are also other empirical methods such as Empirical Modes Decomposition (EMD) method ( 18) that we include in our comparison in this paper.The reason is that both EMD and SCSA share some common properties.Indeed, both of them adapts flexible basis during denoising: EMD approach thresholds the Intrinsic Mode Functions (IMFs) and the SCSA method thresholds the number of squared eigenfunctions.

Simulated Signal and Synthetic Noise
Spectral peaks can be modeled by Gaussian peaks, Lorenz peaks or their combination (14).
We choose multi-peak Gaussian signals as test data.The Gaussian peaks are generated by M is the number of peaks, A i , u i and σ i is the amplitude, position and width of peak i The noise added to the signal is Gaussian noise generated by rand() function in Matlab.In the simulation, M = 1, u 1 = 5, σ 1 = 15, A 1 = 2. Noise levels vary between 1% and 12% with 0.5% interval.
Peak preserving performance: As an evaluation of peak preserving performance of different methods, we use a simulated Gaussian peak to assess peak-preserving performance.
µ is selected according to SCSA method.Depending on the level of noise, µ is usually selected to be smaller when SNR is smaller.In the Savitzky-Golay method, the width of the sliding window is 29 points and polynomial degree is 4. In wavelet method, the function "wden" in Matlab toolbox was used to smooth the signal, the wavelet is Sym4 and decomposed level is 3.We compare the performance of peak-preserving which is described with peak height.The peak height relative error is determined as follows: where M h is the peak amplitude of denoised signal and M c peak amplitude of clean signal.
Denoting the peak width of denoised signal W h and peak width of clean signal W c , the formula of peak width relative error is given by: We compare the performance of peak-preserving. in different noise conditions.The results are shown in Fig. 3.One can easily see that the relative errors of peak heights increase as noise level increases.At the same time, one can also see that the SCSA method has better performance of peak-preserving than other methods in terms of peak preserving performance.

Denoising Performance and comparison:
The denoising performance analysis throughout this paper is assessed by the mean squared error (MSE): and the signal to noise ratio (SNR): , where y is the original signal with length N , y h is the related denoised signal.
In this experiment, a 5-peak simulated signal is generated to demonstrate the denoising process.To make the signal more general, peaks are with different widths and heights with certain peaks overlapping each other.For Savitzky-Golay filter, signal becomes increasingly smooth as the window size increases.On the other way, too broad of a window will reduce the effect of the resolution enhancement and distort the derivative spectra.The best parameters for the Savitzky-Golay method are selected usually by a trial-and-error method (13) (15).In terms of wavelet method, ( 16) presents a selection procedure of mother wavelet basis functions applied for denoising of the noisy signal in wavelet domain while retaining the signal peaks close to their full amplitude.The universal threshold selection by Donoho and Johnstone is applied with varying wavelet basis function, which will also be compared in the following section.
In order to select the best parameter for each method, we optimize each method's parameter at noise level 5%, where we iteratively optimize its parameter using the noisy signal and true signal.
For SG method, consider it having l filter length and r order polynomial, all possible combinations (l, r) are then tried to yield the best smoothing performance.For wavelet method, method is chosen according to (16), where different base functions are compared and selected.µ in SCSA cost function is also selected in a scanning manner.We then fix the parameter and try different noise levels.It can be seen from Fig. 4 that both wavelet method and SCSA based method is out performing the traditional digital filter, while SCSA and wavelet methods are giving comparable denoising results.Both wavelet and Savitzky-Golay methods are optimized as described.sym4 base function with decomposition level 3 is selected for wavelet method, and order 4 with filter length 17 is selected for Savitzky-Golay method.

ECG Signal Denoising
In this section, we will consider the ECG case corrupted by real noise.Real noise records are taken from the MIT-BIH noise stress test database (17).For ECG signal, the analog recordings were played back on a Del Mar Avionics model 660 unit during the digitalization process.The records selected were played back at real time using a specially constructed capstan for the model 660 unit.
The analog outputs of the playback unit were filtered to limit analog-to-digital converter (ADC) saturation and for anti-aliasing, using a passband from 0.1 to 100 Hz relative to real time, which is well beyond the lowest and highest frequencies recoverable from the recordings.The bandpassfiltered signals were digitized at 360 Hz per signal relative to real time.Let n ma (t) and n em (t) be the muscle artifact record and the electrode motion record, respectively.The total noise utilized to corrupt the original clean signal y(t) is as k 1 and k 2 are combined in a way that reaches the same initial SNR in Table 1, with a case example shown in Fig. 5.
Finally, the denoising test is repeated under the same circumstances with different records at different SNRs.The results are presented in Table 1 in terms of SNR after denoising with corresponding methods, with similar way of selecting parameters as in the simulation part.As can be observed here, the Savitsky-Golay method shows less ability to deal with real noise denoising than the SCSA and wavelet-based method.

General signal denoising
Apart from the piecewise-regular signal, three more representative test signals shown in Fig. 9 have been used for validation of the SCSA denoising techniques and others.
While it is proved that C-SCSA is also effective in denoising pulse-shaped signals compared to other popular methods, the C-SCSA real contribution lies in the fact that it can be applied to more general types of signals.Fig. 8 depicts an example of the well studied piecewise-regular signal corrupted by white Gaussian noise corresponding to 31.2446 dB signal-to-noise power ratio (SNR).
As can be seen, while C-SCSA provides an enhancement of the noisy signal, α-SCSA failed in denoising in this case.with four different sampling frequencies to generate four versions per signal, having 256, 512, 1024 and 2048 samples.As before, the results shown correspond to white Gaussian noise generalizations, and in all SCSA-based denoising methods, the penalty parameter is naturally set to make sure that fidelity term and penalty term are of the same order of magnitude.The adopted performance measure is the SNR after denoising, which corresponds to noise levels of 10% before denoising.The performance results for different methods shown correspond to Savitzky-Golay and EMD adaptive interval thresholding.The conclusions drawn from the results are that SCSA-C provide better denoising performance in most cases regardless of the changing of sampling frequencies.

Conclusion
The SCSA method combined with a curvature constraint is proposed as a general method to accomplish peak-preserving smoothing task.Details of the SCSA denoising algorithm and its implementation are given in this work.By designing a proper cost function, we can use the SCSA method to reduce noise while preserving peaks shape.The performance of the proposed method has been investigated and comparison with state of the art methods has been provided.The method not only produces good denoising performance but also guarantees the peaks are well preserved in the denoising process.

Acknowledgement
The research reported in this publication was supported by King Abdullah University of Science and Technology (KAUST) Base Research Fund, (BAS/1/1627-01-01).

Appendix
Let's define (y 1 , y 2 , • • • , y n ) to be input signal y with noise.Let x = y m+1 − y m and w = y m − y m−1 .
Since x and w are subtraction between two neighborhood samples, they resemble noise distributions in homogeneous regions of the signal.Without loss of generality, let's assume x and w to be jointly where ρ = COV {x,w} σ 2 . If the signal are equally spaced with interval denoted as ∆, the curvature defined in (4) can be approximated as: Let C = 1 2πσ 2 √ 1−ρ 2 , we then give the mathematic induction of E{k} as: Let α = x + w, β = x − w, then dαdβ = 2dxdw.Hence: Then Eq. ( 17) reduces to: where η = α/(2∆).

Fig. 1 (
Fig.1(b)  shows an example where {ψ 1h , • • • , ψ 4h } correspond to signal's major peaks.The SCSA analyzing process used in this paper is the standard one.According to this procedure, the SCSA algorithm decomposes the signal as follows.

Figure 1 :
Figure 1: SCSA method.(a) Input Signal with 4 major peaks (b) Squared eigenfunctions of the signal, with only few of them (4 in this case) corresponding to the major peaks and rest of them explaining the details constraints.We refer to this algorithm C-SCSA.The C-SCSA proposes a general solution to reduce signal noise without much pre-knowledge of the characteristics of it.While the h = 3, N h = 212

Figure 2 :
Figure 2: SCSA's application in signal denoising.Input signal with noise(in blue), SCSA spectrum (in red) and residual (in green) (a) Small N h value, not capable of reconstructing the signal.(b)(c) N h increases, the peaks of the signals are recovered, without recovering the rest part of the signal.(d)With N h continuing to increase, the non-peak areas are also recovered, while the noise is separated.The corresponding h is a suitable denoising coefficient.(e)(f ) High N h value will recover the whole signal, including the original signal and noise.

Figure 3 :
Figure 3: Peak preserving performance for single Gaussian peak .(a) Peak height preserving performance (b) Peak width preserving performance (Noise level ranges between 1% and 12% (with interval 0.5%) as shown in the horizontal axis).

Figure 4 :
Figure 4: Quantitative Denoising Performance of Different Noise Levels.Noise level ranges between 0.1% and 15% (with interval 0.1%) as shown in the horizontal axis.Mean Squared error .(b) SNR (dB).Both wavelet and Savitzky-Golay methods are optimized as described.sym4 base function with decomposition level 3 is selected for wavelet method, and order 4 with filter length 17 is selected for Savitzky-Golay method.

Figure 5 :
Figure 5: Enhancement for ECG signals.From top to bottom: (a) Real noise plot, with muscle artifacts (in red) and electrode motion artifacts (in blue).(b) Contaminated ECG signal (SNR = 9.0678 dB).(c) SCSA denoising method (SNR = 11.4315dB).In the last graphs, the reconstructed signal (in blue) and the original signal (in red) are superimposed for comparison purposes.

Table 1 :
Real noise experiment carried out over several records from the MIT BIH arrhythmia database.
d SNR value for the EMD-IT based method.