Remove parabolic or curved trend
6 views (last 30 days)
Show older comments
Konvictus177
on 31 Aug 2022
Commented: Image Analyst
on 31 Aug 2022
Hi,
I collected data from a laser sensor that measures a flat plate with two elevations. Due to the curved lense of the sensor the data has a parabolic trend.
How do I get rid of that parabolic nature of the graph so I can analyze the elevations based on a flat plate/signal? Is there a way to filter it out and make the signal "flat" with the two elevations?
I tried detrending(y_data,1) but it does not filter out a polynomial trend and it does not leave out the real elevations.
The data is attached.
0 Comments
Accepted Answer
Star Strider
on 31 Aug 2022
LD = load(websave('y_data','https://www.mathworks.com/matlabcentral/answers/uploaded_files/1112945/y_data.mat'));
y_data = LD.y_data;
L = numel(y_data)
x = linspace(0, L-1, L);
y_data_dt = detrend(y_data, 5)
figure
yyaxis left
plot(x, y_data)
ylabel('Original')
yyaxis right
plot(x, y_data_dt)
ylabel('Detrended')
grid
Experiment with different polynomial orders to get the result you want. Another option might be to use the bandpass function (ideally with 'ImpulseResponse','iir') depending on the result you want.
.
2 Comments
Star Strider
on 31 Aug 2022
A highpass filter seems to work to eliminate the low-frequency trends without changing the rest of the data. Choose the cutoff frequency to get the result you want.
This uses the first ‘valley’ freqeuency —
LD = load(websave('y_data','https://www.mathworks.com/matlabcentral/answers/uploaded_files/1112945/y_data.mat'));
y_data = LD.y_data.';
L = numel(y_data)
x = linspace(0, L-1, L);
figure
plot(x,y_data)
grid
Fs = 1/(x(2)-x(1))
Fn = Fs/2;
NFFT = 2^nextpow2(L)
FTds = fft(y_data-mean(y_data))/L;
Fv = linspace(0, 1, NFFT/2+1)*Fn;
Iv = 1:numel(Fv);
figure
plot(Fv, abs(FTds(Iv))*2, 'DisplayName','Fourier Transform')
grid
xlim([0 0.1])
[pks,plcs] = findpeaks(abs(FTds(Iv))*2, 'MinPeakProminence',0.00025);
[vys,vlcs] = findpeaks(-abs(FTds(Iv))*2, 'MinPeakProminence',0.00025);
hold on
plot(Fv(plcs), pks, '^r', 'DisplayName','Peaks')
plot(Fv(vlcs), -vys, 'vr', 'DisplayName','Valleys')
hold off
xlabel('Frequency')
ylabel('Magnitude')
legend('Location','best')
Fvpk = Fv(plcs); % Peak Frequencies
Fvvy = Fv(vlcs); % Valley Frequencies
% BPv = [Fvvy([1 end])] % Peak Frequencies ± Offset To Use In 'bandpass' Call
N = 50; % Padding Vector Length
% y_data_filt = bandpass([zeros(N,1)+mean(y_data); y_data], BPv, Fs, 'ImpulseResponse','iir'); % Filter Signal
y_data_filt = highpass([zeros(N,1)+y_data(1); y_data], Fvvy(1), Fs, 'ImpulseResponse','iir'); % Filter Signal
y_data_filt = y_data_filt(N+1:end); % Adding Vector Of Mean Values Eliminates Initial Filter Transient, Remove Those Values Later
figure
plot(x,y_data,'-b', 'DisplayName','Original')
hold on
plot(x,y_data_filt, '-r', 'DisplayName','Filtered')
hold off
grid
legend('Location','best')
This is the best I can do with this signal.
.
More Answers (1)
Image Analyst
on 31 Aug 2022
Edited: Image Analyst
on 31 Aug 2022
If you have the Computer Vision Toolbox you can fit a quadratic to the lower part of your curve, without explicitly telling it what elements those are, by using fitPolynomialRANSAC
s = load('y_data.mat')
y = s.y_data;
subplot(2, 1, 1);
plot(y, 'b.-', 'LineWidth', 2);
hold on;
xlabel('Index', 'FontSize',fontSize);
ylabel('y', 'FontSize',fontSize);
grid on;
x = 1 : length(y);
xy = [x(:), y(:)]
N = 2; % second-degree polynomial
maxDistance = 1; % maximum allowed distance for a point to be inlier
coefficients = fitPolynomialRANSAC(xy, N, maxDistance)
yFitted = polyval(coefficients, x);
plot(yFitted, 'r.-', 'LineWidth', 2);
% Subtract
baseLineCorrectedY = y - yFitted
subplot(2, 1, 2);
plot(baseLineCorrectedY, 'b.-', 'LineWidth', 2);
hold on;
grid on;
xlabel('Index', 'FontSize',fontSize);
ylabel('y', 'FontSize',fontSize);
2 Comments
Image Analyst
on 31 Aug 2022
Yeah. Hmmm...I guess that's just because of the way RANSAC works by sampling random points. You could put it in a loop and call it, say 50 times and take the one set of y values where the mean absolute corrected y value is closest to 0.
% Demo by Image Analyst
% 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 = 18;
s = load('y_data.mat')
y = s.y_data;
subplot(2, 1, 1);
plot(y, 'b.-', 'LineWidth', 2);
hold on;
xlabel('Index', 'FontSize',fontSize);
ylabel('y', 'FontSize',fontSize);
grid on;
x = 1 : length(y);
xy = [x(:), y(:)];
N = 2; % second-degree polynomial
maxDistance = 1; % maximum allowed distance for a point to be inlier
numTrials = 50;
meanAbsValue = zeros(1, numTrials);
% Try RANSAC 50 times and we'll choose the best.
subplot(2, 1, 2);
for k = 1 : numTrials
coefficients{k} = fitPolynomialRANSAC(xy, N, maxDistance);
yFitted = polyval(coefficients{k}, x);
% Subtract
baseLineCorrectedY = y - yFitted;
plot(baseLineCorrectedY, '-', 'LineWidth', 1);
hold on;
grid on;
xlabel('Index', 'FontSize',fontSize);
ylabel('y', 'FontSize',fontSize);
% Compute mean value
meanAbsValue(k) = mean(abs(baseLineCorrectedY));
end
% Find the closest
[~, index] = min(meanAbsValue);
yFitted = polyval(coefficients{index}, x);
% Subtract
baseLineCorrectedY = y - yFitted;
subplot(2, 1, 1);
plot(yFitted, 'r.-', 'LineWidth', 2);
% Plot final baseline corrected signal on a separate figure.
figure;
plot(baseLineCorrectedY, 'r.-', 'LineWidth', 2);
grid on;
xlabel('Index', 'FontSize',fontSize);
ylabel('y', 'FontSize',fontSize);
% Put a line along 0
yline(0, 'LineWidth', 2)
See Also
Categories
Find more on Bartlett in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!