Noise Reduction

This tutorial covers the denoising capabilities of AutoSignal. Noise can be removed from data using Fourier, Eigen, and Wavelet decomposition and reconstruction techniques and for low frequency data, the time-domain Savitzky-Golay smoothing filter is offered.

Generating a Noisy Test Signal

Select the Generate Signal option in the Edit menu or Main toolbar.

For this tutorial, we will create a sinusoid and add Gaussian noise whose standard deviation is exactly equal to the sinusoid's standard deviation. The Signal to Noise (S/N) ratio is 1.0 or 0 dB.

Click Read and select the file tutor6a.sig from the Signals subdirectory.

The following signal expression is imported:

SRATE=5000

NYQ=SRATE/2

AMP=100

FREQ=NYQ*0.1

PHASE=PI

Y=AMP*SIN(2*PI*X*FREQ+PHASE)

The X (time) values vary from 0 to 0.2047 with a 0.0002 sample increment. The Nyquist frequency is thus 2500 (half the 5000 sampling rate). White noise is added at a 100% level. The aim is to explore which of the denoising procedures can best recover the underlying sinusoid (frequency=250, amplitude=100, phase=p).

Click OK to process the current signal.

An AutoSignal graph is presented containing the 1024 point generated data. With the noise equal to the data, it is hard to discern the presence of a single sinusoid.

Click OK to accept the generated data. Click Yes when asked to update the main data table with the revised data.

Time Domain Denoising: Savitzky-Golay

The first approach to denoising will be to apply a Savitzky-Golay smoothing filter. Time domain smoothing procedures are limited to low frequency signals since they cannot interpolate a wave from a handful of points. Rather, the shape of the wave must be sufficiently defined.

Select the Savitzky-Golay Smoothing Filter option from the Time menu or toolbar. Set the one-sided size of the moving window (Win n) to 5. Set the Order to 4 for quartic fitting and set the sequential pass count (Passes) to 3. Be sure Data is selected.

When smoothing oscillatory data, an order of 4 or higher is recommended. Also, three sequential passes of the filter is about the optimum in trading off noise reduction against amplitude attenuation.

The Equivalent Noise % is based upon a time-domain cubic polynomial interpolation that uses the two points on each side of a given point, and then computes the difference between the actual data and predicted data at the central point. This is an approximation and valid only when the signal consists of low frequency information. In this case, the noise estimate should be in the vicinity of 100%. The values you see in this tutorial will be highly dependent on the random noise sequence since it has such a dominant presence in this data set. For the data generated in the examples shown, these Savitzky-Golay settings resulted in an estimated white noise level of 1.53%. Only 1.41% of the noise originally present remains.

To assess how closely the sinusoid was recovered in the smoothing, we will explore the Fourier spectrum of the output data and we will non-linearly fit a single sinusoid model to the output data. This will enable a separate assessment of noise reduction as well as offering an indication of how accurately the parameters are recovered.

Fourier Analysis and Non-Linear Sinusoidal Fit

Click OK to close the Savitzky-Golay procedure and answer Yes to update the main data table.

Select the Fourier Spectrum with Data Window option from the Spectral menu or toolbar. Select the Kaiser-Bessel window, set the window's one-sided spectral bin width (adjust) to 4. Use the Best Exact N algorithm with Nmin set to the data length 1024. Change the plot to dB Norm, and set the signal count in the sig field to 1.

Note that the central frequency is readily recovered and that a -25 dB noise floor exists at the lower frequencies. The filtered data hardly generate a clean spectral peak.

Click the Non-Linear Optimization button. Accept the default settings and click the OK button. At the conclusion of the iterative fit, click the Review Fit button.

Click the Numeric Summary button in the Non-Linear Optimization Review. Inspect the fit results.

Fitted Parameters

r² Coef Det      DF Adj r²      Fit Std Err      F-value

0.81085336      0.81029705      33.3885192      2188.46419

  Data Power       Model Power       Error Power

1203.5147860      976.11767392      227.64077441

 Comp      Type       Frequency       Amplitude       Phase       Damping       Power       %

    1      Sine      249.876652      97.7059488      3.25985841            975.479450      100.000000

Although the frequency was accurately recovered by the non-linear sinusoid fit, there is a small attenuation of the amplitude, and a minor shift in the phase. The most important values are the error power and r² goodness of fit values. Although this procedure produced a very smooth data stream, the filtered data are not well fitted by a sinusoidal model. Although the procedure did not distort the underlying sinusoidal parameters that can be recovered by non-linear modeling, the smoothing did produce a very distorted sinusoid.

Close the Numeric Summary and close the Non-Linear Optimization Review.

Close the Fourier spectrum window with the OK button since we want to have these same settings available for subsequent analyses.

Click the Reset XY Data button in the main toolbar to restore the original unfiltered data.

The Savitzky-Golay procedure can produce very high levels of smoothing. The best sinusoidal fit achieved with this data set occurred when a window width of 10 was used. The estimated white noise in the filtered signal was only 0.32%

Fitted Parameters

r² Coef Det      DF Adj r²      Fit Std Err      F-value

0.87988447      0.87953119      23.5908541      3739.57496

  Data Power       Model Power       Error Power

946.11496329      832.71767699      113.64309894

 Comp      Type       Frequency       Amplitude       Phase

    1      Sine      249.868497      90.2434929      3.26624311

The error power is halved and the goodness of fit values have increased. The filtered data set thus looks more like a sinusoid. At this level of smoothing, however, the amplitude is already seeing appreciable attenuation. Although the Savitzky-Golay is probably the most commonly used and the most powerful of the time-domain noise reduction methods, it is not typically a good algorithm for harmonic signals.

Frequency Domain Denoising: dB Thresholding

Select the Fourier Smoothing and Denoising option in the Process menu or toolbar. Select a dB threshold. Left click the mouse just below the first sidelobes adjacent the main peak. The mouse should be at about -16 dB.

Only the bin containing the main peak and the two adjacent bins should be active.

Here, the filtered signal appears sinusoidal and with very little variation from the target amplitude of 100. For the data in this example, the estimated noise is 0.295%; only 0.27% of the original noise remains. The power in this example was diminished by about half.

Click OK to close the Fourier denoising procedure and answer Yes to update the main data table.

Select the Fourier Spectrum with Data Window option.

A very clean spectral peak is now present.

Click the Non-Linear Optimization button. Click the OK button. Click the Review Fit button at the conclusion of the fit.

Click the Numeric Summary button in the Non-Linear Optimization Review. Inspect the fit results.

Fitted Parameters

r² Coef Det      DF Adj r²      Fit Std Err      F-value

0.96723597      0.96713960      12.3581081      15070.6104

  Data Power       Model Power       Error Power

951.83655223      920.65054905      31.186003175

 Comp      Type       Frequency       Amplitude       Phase       Damping       Power       %

    1      Sine      249.595685      94.8718693      3.44237876            920.156884      100.000000

The goodness of fit indices and error power readily confirm a much more sinusoidal signal. The unfortunate side is that the parameters are not well recovered. Although this approach is an intuitive one, it is not necessarily a wise course. The amplitude is attenuated, a phase shift is present, and the frequency is shifted downward.

Close the Numeric Summary and close the Non-Linear Optimization Review.

Close the Fourier spectrum window.

Click the Reset XY Data button in the main toolbar to restore the original unfiltered data.

Frequency Domain Denoising: dB Thresholding

Select the Fourier Smoothing and Denoising option in the Process menu or toolbar. Select a Frequency threshold. Left click the mouse at a frequency of about 300.

Note that all of the noise in the higher frequencies is being zeroed.

The signal appears more variable.

Click OK to close the Fourier denoising procedure and answer Yes to update the main data table.

Select the Fourier Spectrum with Data Window option.

It is apparent that frequency thresholding is an incomplete noise removal approach.

Click the Non-Linear Optimization button. Click the OK button. Click the Review Fit button at the conclusion of the fit.

Click the Numeric Summary button in the Non-Linear Optimization Review. Inspect the fit results.

Fitted Parameters

r2 Coef Det      DF Adj r2      Fit Std Err      F-value

0.89896613      0.89866897      23.1340959      4542.26113

  Data Power       Model Power       Error Power

1081.6676056      972.62923165      109.28506187

 Comp      Type       Frequency       Amplitude       Phase

    1      Sine      249.850258      97.5296645      3.27792823

The parameters are now less distorted, but the error power has increased.

Close the Numeric Summary and close the Non-Linear Optimization Review.

Close the Fourier spectrum window.

Click the Reset XY Data button in the main toolbar to restore the original unfiltered data.

Eigendecomposition Denoising

Fourier denoising depends on the ability to attribute all spectral information beyond a certain frequency or below a certain spectral threshold to noise. Because of spectral leakage from edge effects and components failing to center within a frequency bin, the signal is not likely to be conveniently confined to certain bins while noise is confined exclusively to the remainder. In practice, the signal/noise separation is incomplete.

An eigendecomposition uses adaptive non-parametric basis functions, partitioning a data stream by signal strength rather than frequency. The highest power components distribute amongst the initial eigenmodes while the noise is found in the subsequent eigenmodes. Denoising using eigendecomposition procedures is usually the best way to remove noise from stationary data. Further, since the eigendecomposition is non-parametric and adaptive, non-stationary signal components undergoing continuous and smooth change can often be mapped as well.

Two eigenmodes are needed to capture an oscillation of a given frequency. That oscillation can be a sinusoid, an exponentially damped sinusoid, a sawtooth, a square wave, or any other shape. If the amplitude of the oscillation varies in accord with a continuous function of time, the two eigenmodes will still map this non-stationary oscillation. On the other hand, variable frequency oscillations such as chirps may require a large number of eigenmodes to fully map the non-stationary data trend.

Select the Eigendecomposition Smoothing and Denoising option from the Process menu or toolbar. Be sure the CovM FB (forward-backward covariance matrix) algorithm is selected and set the Order of the eigendecomposition matrix to 50. To accommodate the single oscillatory component, set the Signal Space to 2.

The estimated noise remaining after reconstructing the first and second eigenmodes is about 0.31% or 0.29% of the noise originally present. The filtered curve appears to be very close to an amplitude 100 sinusoid.

Click OK to close the Eigendecomposition denoising procedure and answer Yes to update the main data table.

Select the Fourier Spectrum with Data Window option.

While the spectrum lacks the sharp spectral peak that can be attained using Fourier dB thresholding, the noise level falls off rapidly to the third or fourth decimal place and the sinusoidal frequency is determined to a very good accuracy.

Click the Non-Linear Optimization button. Click the OK button. Click the Review Fit button at the conclusion of the fit.

Click the Numeric Summary button in the Non-Linear Optimization Review. Inspect the fit results.

Fitted Parameters

r² Coef Det      DF Adj r²      Fit Std Err      F-value

0.97579947      0.97572829      10.9204046      20584.0817

  Data Power       Model Power       Error Power

1006.2557375      982.23773674      24.351919422

 Comp      Type       Frequency       Amplitude       Phase

    1      Sine      249.862078      98.0108563      3.27012808

This is the best noise removal thus far encountered in this exercise. The error power is the lowest and the goodness of fit values are the highest. The sinusoidal frequency is accurately recovered and the amplitude shows the least attenuation. A minor phase variation is still present.

Close the Numeric Summary and close the Non-Linear Optimization Review.

Close the Fourier spectrum window.

Click the Reset XY Data button in the main toolbar to restore the original unfiltered data.

Higher orders can sometimes offer modest benefits since additional eigenmodes are available to more fully capture the noise, and less noise leaks into the signal space eigenmodes. The following results are for a matrix order of 100:

Fitted Parameters

r² Coef Det      DF Adj r²      Fit Std Err      F-value

0.98553429      0.98549174      8.38522847      34779.8439

  Data Power       Model Power       Error Power

992.53466372      978.52164368      14.357721952

 Comp      Type       Frequency       Amplitude       Phase

    1      Sine      249.849401      97.8250880      3.27638141

Although the parameters are essentially unchanged, the error power has been cut over 40%. The next results are for a matrix order of 200:

Fitted Parameters

r² Coef Det      DF Adj r²      Fit Std Err      F-value

0.99209635      0.99207311      6.13716693      64079.9474

  Data Power       Model Power       Error Power

973.11497638      965.73326483      7.6911558252

 Comp      Type       Frequency       Amplitude       Phase

    1      Sine      249.880309      97.1851563      3.25687527

The error power is halved once again and the goodness of fit values further improved. With some patience, the following results can be achieved with AutoSignal's maximum matrix order of 300:

Fitted Parameters

r² Coef Det      DF Adj r²      Fit Std Err      F-value

0.99659047      0.99658044      3.95021420      1.4922e+05

  Data Power       Model Power       Error Power

934.54977462      931.67800056      3.1863760497

 Comp      Type       Frequency       Amplitude       Phase

    1      Sine      249.922908      95.4588491      3.22616749

The error power is now minimal and the goodness of fit is exceptional. The signal is now almost completely free of noise.

Note that SVD matrix decomposition is a particularly slow N³ process. Whereas an order of 100 can require a few seconds, an order of 300 can require a minute or more of processing time. Also, bear in mind that such large matrices consume considerable memory. If insufficient physical memory is available and the computations must use disk based memory, the processing can become extremely slow.

Denoising by Parametric Modeling of Harmonics

When a signal is composed of one or more sinusoids in noise, there is often no need to denoise the data prior to parametric fitting. The parametric model is itself a noise-free representation of the signal.

Select the Parametric Interpolation and Prediction option in the Process menu or toolbar. Be sure the algorithm is set to AR,Data FB, that the order is set to 40, that the Signal Subspace is set to 2, that the Undamped model is selected, and that Enable is checked for NL Optimization. Zoom-in a portion of the reconstructed signal in the lower graph.

Click the Numeric Summary button and inspect the fit results.

Fitted Parameters

r² Coef Det      DF Adj r²      Fit Std Err      F-value

0.49362519      0.49213586      70.1667999      497.646522

  Data Power       Model Power       Error Power

1985.3952949      980.28782485      1005.3541579

 Comp      Type       Frequency       Amplitude       Phase       Power

    1      Sine      249.879086      97.9145969      3.25809353      979.645852

Although the goodness of fit is poor due to the equal measure of noise present, the parameters are recovered with the same accuracy as the denoised fits. You can think of this approach as a least-squares noise filter for signals consisting of sinusoidal components.

Note also that this error power of 1005 represents the noise power for unfiltered data. The highest order eigendecomposition denoising cut this noise power from 1005 to 3.19 (99.61% of the noise was removed).

Close the Numeric Summary.

Close the Parametric Interpolation and Prediction procedure.

We will now explore the denoising capabilities of AutoSignal with non-stationary data.

Generating a Chirp Signal

Select the Generate Signal option in the Edit menu or Main toolbar.

Click Read and select the file tutor6b.sig from the Signals subdirectory.

The following signal expression is imported:

RATE=2.5

Y=SIN(2*PI*X*(25+50*(1-EXP(-RATE*X)))+PI/2)

The X (time) values vary from 0 to 3 with a 0.003 sample increment. White noise is added at a 50% level.

Click OK to process the current signal.

This chirp consists of a sinusoid ramping up in frequency at a rate based upon first order kinetic formation. At the beginning of the sample the sinusoid's frequency is 25. The frequency ramps up to 75, reaching steady state in the middle of the data sequence. The generated data sequence consists of 1000 values.

Click OK to accept the generated data. Click Yes when asked to update the main data table with the revised data.

Wavelet Denoising

Removing noise using wavelet thresholding is not as straightforward as doing so with Fourier filtering. The key to successful denoising using wavelet decomposition is visualizing the signal/noise threshold. AutoSignal simplifies this visualization.

Select the Wavelet Smoothing and Denoising option in the Process menu or toolbar. Be sure the wavelet is set to Morlet and that the Complex box is checked. Set the Morlet's adjustable parameter (Adj) to 12. Set the dB lim value to 48. Select the dB threshold and enter a value of 48.

Click the Modify Contour Properties button in the upper graph's toolbar. Be sure the grid is set to 60 and that the 48 Spectrum gradient is selected. Click OK.

When the dB lim value is set at 48, all spectral information from 0.004% amplitude (-48dB) to maximum amplitude (0dB) is rendered within the gradient. All spectral information below -48dB is rendered at the base color. In a spectrum plot, this will be a dark red color. In this case where a 48dB range is plotted across 48 color gradients, each color represents a 1dB range in the spectrum.

The chirp is clearly evident in the plot, but so is a great deal of noise across all frequencies and times.

Set the dB lim value to 12.

Only a trace of the noise remains in the plot. This type of visualization establishes that -12dB is a good signal/noise threshold. Note that the dB lim adjustment impacts only the plot. In order to adjust the wavelet threshold, the Threshold value must also be set.

Set the dB Threshold value to 12.

Although the signal's +1 to -1 amplitude range is not perfectly recovered, most of the noise has been removed.

This initial release of AutoSignal supports only sinusoid and damped sinusoid models for non-linear parametric fitting. To confirm the effectiveness of this noise removal, the above output was fitted in TableCurve 2D to a user-defined function consisting of this chirp model with the amplitude, starting frequency, frequency delta, and phase as adjustable parameters:

r² Coef Det      DF Adj r²      Fit Std Err      F-value

0.9625197162      0.9623690417      0.1295811089      8525.9905561

 Parm      Value      Std Error      t-value      95% Confidence Limits

 a      0.927103186      0.005794401      159.9998393      0.915732552      0.938473821

 b      25.09529394      0.042548631      589.8026185      25.01179869      25.17878918

 c      49.90976419      0.040488523      1232.689180      49.83031159      49.98921679

 d      1.477775140      0.035085410      42.11936315      1.408925333      1.546624947

The target values are 1.0, 25, 50, and 1.57. Although the amplitude and phase parameters weren't recovered as well as the two frequency parameters, the high goodness of fit values confirm a very successful denoising.

Click OK to close the wavelet denoising procedure and answer No to update the main data table.