Smart detrending for both real noisy signal and pure sinus

In my problem I want to detrend a signal. Today I am using the MATLAB smooth function with the Savitzky-Golay method. The detrending procedure looks like

n = length(u);
uS = smooth(u,n,'sgolay')';
uM = u - uS; % detrend
uM = NLM_1D(uM,1,smoothP,1); % denoise using NLM algorithm

This code works pretty well for real noisy signal. However it fails on pure sinus signal

u = u0 + Q*sin(w*t)

Not only it fails to recognize that the sinus needs only "minus constant" detrending (in the plot below the red curve should be a flatline uS=0), it also adds a trend so that the signal actually increases with time. This false trend impairs further calculations I do with the "detrended" signal.

Blue: pure sinus. Red: the trend obtained by smooth. Green: uM = u - uM, detrended and denoised

The following example show what I seek: an signal is given (blue, in the plot below is y=10*exp(-0.1*x)+sin(x) ), I want to find the baseline trend (red, 10*exp(-0.1*x) ) and remain with the oscillatory component (green, sin(x) ). In this example I know the formula for the input signal, but I want to do it numerically even for signal I don't know the formula for.

Is there a way to smart detrend a signal so it will also work well on pure sinus signal? If not, should I split the code (e.g. by an "if" condition) to work separately on sinus (how do I recognize a pure sinus numerically?) and separately on real noisy signal?

[Greg Dionne]

Could you post a sample of your data?

[Evenor]

u0 = -15, Q=5, w=0.05 or w=0.1 or w=0.01.

[Greg Dionne]

I see. In your case the sinusoid is dominating. Can you post (upload) your actual data as a .mat file? Otherwise all we can do is give you a few simple suggestions to try:

• If you know you just have one (pure) frequency of interference, you can use meanfreq in the Signal Processing Toolbox to estimate it. Once you have that, you can do a least squares fit sine/cosine pair of that particular frequency to remove it.
• Alternatively you could try a combination of fsst, tfridge, and ifsst if your frequency interference is non-constant. The example here has background on how to use them.
• Is your wanted signal in a particular frequency band? If so, you can try constructing a bandpass filter to extract it (taking out low frequency trends). You can also construct a notch filter if you know you have a particular problematic frequency range (see this example).
• If you can tweak your sample rate so that interfering frequencies are all an integral multiple of the sampling rate, you can subtract off a moving sum to filter out those components.

Let us know.

[Evenor]

On real data there is no problem. The problem is on the artificial signal generated by the formula u = u0 + Q*sin(w*t) with u0=-15, Q=5, w=0.1 and t=0:200 and t=0:3600.

Answer 1

[Image Analyst]

Try sgolayfilt() and use order 3 or 5 since the Taylor series of a sine is close to a polynomial. Don't have the frame width be too wide (like covering several humps) or else you won't get a polynomial to fit that well. I'd use around 10 - 50 percent of the width of a single hump. Experiment until it looks good.

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 25;

% Make a noisy sine wave signal
x = 1 : 300;
period = 50
y = sin(2*pi*x/period);
noiseAmplitude = 0.8;
y = y + noiseAmplitude * rand(size(y));
subplot(2,1,1);
plot(x, y, 'b-', 'LineWidth', 2);
grid on;
title('Noisy Signal', 'FontSize', fontSize);
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);
% Give a name to the title bar.
set(gcf, 'Name', 'Demo by ImageAnalyst', 'NumberTitle', 'Off') 

% Now smooth with a Savitzky-Golay sliding polynomial filter
windowWidth = 27
polynomialOrder = 3
smoothY = sgolayfilt(y, polynomialOrder, windowWidth);
subplot(2,1,2);
plot(x, smoothY, 'b-', 'LineWidth', 2);
grid on;
title('Smoothed Signal', 'FontSize', fontSize);

[Evenor]

Thank you for the answer. However, my problem is not denoising a signal but detrending it so only the oscillatory component remains. For example: for a pure sine u = u0 + Q*sin(w*t) the trend is a a flat-line y=u0; For u = A*exp(-t) + sin(w*t) the trend should be A*exp(-t) . For u = a1*exp(-b1*t) + ... + an*exp(-bn*t) + sin(w*t) is a1*exp(-b1*t)+...+an*exp(-bn*t).

Example: detrending y = 10*exp(-0.1*x)+sin(x) in blue, the trend is in red and the oscillatory component is in green.

I tried to find the trend by "smoothing out" the oscillations to remain only with the trend. After finding the trend I can subtract it from the signal.

[Image Analyst]

Evenor, you can easily do this with fitnlm(). Fit the data to an exponential decay, then subtract the fit from the actual data. See below for a demo:

% Initialization steps.
clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 20;

   % Uses fitnlm() to fit a non-linear model (an exponential decay curve) through noisy data.
  % Requires the Statistics and Machine Learning Toolbox, which is where fitnlm() is contained.
  % Initialization steps.
  clc;    % Clear the command window.
  close all;  % Close all figures (except those of imtool.)
  clear;  % Erase all existing variables. Or clearvars if you want.
  workspace;  % Make sure the workspace panel is showing.
  format long g;
  format compact;
  fontSize = 20;

% Create the X coordinates as 1000 samples from 0 to 100.
X = linspace(0, 100, 1000);

% Define function that the X values obey.
a = 10; % Arbitrary sample values I picked.
b = -0.1;
period = 6.6;
Y = a * exp(b * X) + sin(2 * pi * X / period);

% Now we have noisy training data that we can send to fitnlm().
% Plot the noisy initial data.
plot(X, Y, 'b-', 'LineWidth', 2, 'MarkerSize', 15);
grid on;

% Convert X and Y into a table, which is the form fitnlm() likes the input data to be in.
tbl = table(X', Y');
% Define the model as Y = a + exp(-b*x)
% Note how this "x" of modelfun is related to big X and big Y.
% x((:, 1) is actually X and x(:, 2) is actually Y - the first and second columns of the table.
modelfun = @(b,x) b(1) * exp(b(2)*x(:, 1)) + b(3);  
beta0 = [11, -.1, 6]; % Guess values to start with.  Just make your best guess.
% Now the next line is where the actual model computation is done.
mdl = fitnlm(tbl, modelfun, beta0);
% Now the model creation is done and the coefficients have been determined.
% YAY!!!!

% Extract the coefficient values from the the model object.
% The actual coefficients are in the "Estimate" column of the "Coefficients" table that's part of the mode.
coefficients = mdl.Coefficients{:, 'Estimate'}
% Create smoothed/regressed data using the model:
yFitted = coefficients(1) * exp(coefficients(2)*X) + coefficients(3);
% Now we're done and we can plot the smooth model as a red line going through the noisy blue markers.
hold on;
plot(X, yFitted, 'r-', 'LineWidth', 2);
grid on;

% Compute the detrended data
detrendedY = Y - yFitted;
plot(X, detrendedY, 'g-', 'LineWidth', 3);

title('Exponential Regression with fitnlm()', 'FontSize', fontSize);
xlabel('X', 'FontSize', fontSize);
ylabel('Y', 'FontSize', fontSize);
legendHandle = legend('Training Y', 'Fitted Y', 'Detrended Y', 'Location', 'northeast');
legendHandle.FontSize = 30;
formulaString = sprintf('Y = %.3f * exp(%.3f * X) + %.3f', coefficients(1), coefficients(2), coefficients(3))
text(7, 11, formulaString, 'FontSize', 25, 'FontWeight', 'bold', 'Color', 'r');

% Set up figure properties:
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);
% Get rid of tool bar and pulldown menus that are along top of figure.
% set(gcf, 'Toolbar', 'none', 'Menu', 'none');
% Give a name to the title bar.
set(gcf, 'Name', 'Demo by ImageAnalyst', 'NumberTitle', 'Off')

[Evenor]

Thank you for your answer. However, it assumes that the trend is Y = a*exp(-b*X)+c. I don't want to assume that. Models like Y = c + a1*exp(-b1*X) + ... + an*exp(-bn*X) are also possible. I want to detrend without assuming that I know the formula for the trend.

[Image Analyst]

Well now we don't know what you don't know. You have to know something or else you'll just have to assume everything is signal and nothing is unwanted offset. For example, if you know that you have a sinusoid then you can just to a band pass filter. But if you know absolutely nothing about your signal I think you just have to assume that it's ALL signal.

[Bruno Luong]

IMA: "For example, if you know that you have a sinusoid then you can just to a band pass filter"

No need such restricted constraints: meaning that the signal is mono frequency (e.g, sine wave).

All he needs to know (set) the maximum frequency of the "trend" he wants to remove.

Answer 2

[Bruno Luong]

"Detrending" just looks like a highpass filter to me. There is no such thing as "smart" xxx (even a lot of people pretent there is out there) especially when a computer takes a decision of your parameters.

If you want to detrend correctly, then you have to design a filter with a cut-off frequency that you think that can separates a trend (lower f) and signal+noise (higher f).

[Evenor]

I wrote a low-pass filter (I don't have the Signal Processing Toolbox). I then analyzed a decaying exponent signal ("the trend"). Its spectrum has a non non-negligible tail which overlaps the oscillations frequency. When I filter out the oscillation using a low-pass filter I also filter out some of the high frequency components of the trend (a decaying exponential in this case) so the filtered signal is a distorted decaying exponent, mainly on the endpoints.

See example below:

The red curve is for the original signal x = exp(-10*t) from 0 to 1. The blue curve is for the signal after filtering out frequencies higher than 40 Hz (see right plot, red - org. spectrum, blue - after low-pass filter). As you can in the left plot, the filtered signal has distortions, mainly at the endpoints (e.g. it starts oscillating near t=1). The distortions are clearer if I take the cutoff frequency as 20 Hz. See plot below.

[Bruno Luong]

It doesn't look like your filter(s) do the job

[Evenor]

If I take cutoff frequency as the maximal frequency appearing in the discrete spectrum so that practically there is no filtering, I recover the decaying exponent almost perfectly.

It also filters out well two-component signal:

Fs = 1023;
t=0:(1/Fs):1;
x = sin(2*pi*t) + sin(2*pi*10*t);
z = LowPassFilter(t,x,Fs,5);

It seems to work on simple input signal.

[Evenor]

I tried to write another filter - DespikeFilter. This filter searches for high-frequency spikes in the power spectrum and replaces the corresponding amplitude y(j) with a local mean (e.g. y(j)=(y(j-1)+y(j+1))/2 ). This works well for delta functions (spikes that are point wide) but sometimes even a pure sinus does not produces a point wide spike as it should have.

[Bruno Luong]

When working with discrete data on finite-length interval, you should never expect a sharp spectrum even for mono-frequency signal.

Answer 3

[Bruno Luong]

% Create the X coordinates as 1000 samples from 0 to 100.
X = linspace(0, 100, 1000);

% Define function that the X values obey.
a = 10; % Arbitrary sample values I picked.
b = -0.1;
period = 6.6;
Y = a * exp(b * X) + sin(2 * pi * X / period);
Y = Y + 0.5*randn(size(Y));

Y2 = sin(2 * pi * X / period);

hpFilt = designfilt('highpassfir','StopbandFrequency',1e-10, ...
         'PassbandFrequency',20/length(X),'PassbandRipple',1, ...
         'StopbandAttenuation',33,'DesignMethod','kaiserwin');
detrend = @(y) filtfilt(hpFilt, y);

close all

subplot(2,1,1);
plot(X, Y, 'b-', X, detrend(Y), 'r');
grid on;
legend('input signal', 'high pass signal');

subplot(2,1,2);
plot(X, Y2, 'b-', X, detrend(Y2), 'r');
grid on;
legend('input signal', 'high pass signal');