Smart detrending for both real noisy signal and pure sinus

12 visualizzazioni (ultimi 30 giorni)
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?
  4 Commenti
Greg Dionne
Greg Dionne il 17 Ott 2018
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
Evenor il 21 Ott 2018
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.

Accedi per commentare.

Risposte (3)

Image Analyst
Image Analyst il 21 Ott 2018
Modificato: Image Analyst il 21 Ott 2018
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);
  5 Commenti
Image Analyst
Image Analyst il 29 Ott 2018
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
Bruno Luong il 29 Ott 2018
Modificato: Bruno Luong il 29 Ott 2018
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.

Accedi per commentare.


Bruno Luong
Bruno Luong il 22 Ott 2018
"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).
  5 Commenti
Evenor
Evenor il 28 Ott 2018
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
Bruno Luong il 28 Ott 2018
When working with discrete data on finite-length interval, you should never expect a sharp spectrum even for mono-frequency signal.

Accedi per commentare.


Bruno Luong
Bruno Luong il 29 Ott 2018
% 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');

Categorie

Scopri di più su Signal Generation, Analysis, and Preprocessing in Help Center e File Exchange

Prodotti


Release

R2014a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by