Remove parabolic or curved trend

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.

Answer 1

[Star Strider]

Use the detrend function —

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)

L = 799

x = linspace(0, L-1, L);

y_data_dt = detrend(y_data, 5)

y_data_dt = 1×799

0.0028 0.0027 0.0025 0.0023 0.0022 0.0020 0.0018 0.0016 0.0014 0.0013 0.0011 0.0010 0.0008 0.0006 0.0004 0.0003 0.0001 -0.0001 -0.0002 -0.0002 -0.0001 -0.0001 -0.0000 -0.0000 0.0000 0.0001 0.0000 -0.0001 -0.0001 -0.0001

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.

.

[Konvictus177]

Thanks for the quick reply.

I want to measure the distance to both elevations from my base profile which should be around 10 microns.

Simple detrending seems not to be the right thing to do in my opinion since detrending also removes the trend caused by the elevations and higher detrending orders will detrend more then they should.

[Star Strider]

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)

L = 799

x = linspace(0, L-1, L);

figure
plot(x,y_data)
grid

Fs = 1/(x(2)-x(1))

Fs = 1

Fn = Fs/2;
NFFT = 2^nextpow2(L)

NFFT = 1024

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.

.

Answer 2

[Image Analyst]

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);

[Konvictus177]

Looks promising. Just tried it out quickly but my results looks a bit different than yours. Used exactly the code in your answer.

Actually each time I run the code I will get a different result.

[Image Analyst]

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)