Remove parabolic or curved trend
7 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
Konvictus177
il 31 Ago 2022
Commentato: Image Analyst
il 31 Ago 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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1112960/image.png)
0 Commenti
Risposta accettata
Star Strider
il 31 Ago 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 Commenti
Star Strider
il 31 Ago 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.
.
Più risposte (1)
Image Analyst
il 31 Ago 2022
Modificato: Image Analyst
il 31 Ago 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);
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1113065/image.png)
2 Commenti
Image Analyst
il 31 Ago 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)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1113320/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1113325/image.png)
Vedere anche
Categorie
Scopri di più su Bartlett in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!