Fit data to sawtooth shape

18 次查看（过去 30 天）
Angela
Angela 2018-6-21
评论： Daniel M 2019-10-18
I have a data sample and i would like to fit it to a sawtooth shape but it does not seem to be one of the functions that Matlab supports in the curve fitting. Does anyone know how to do that? I would like to do the fit inside a program (not interactively).
  4 个评论
显示 2更早的评论
Angela
Angela 2019-10-17
Unfortunately I was never able to do that through Matlab, although the shape I was trying had uniform stride and slope but they were both unknown.
xi
xi 2019-10-17
write a function of your sawtooth model with any given stride and slope(s) as parameters. Optimize these unknow parameters such that the difference between the real curve and your model curve is minimum such as least squre error.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（2 个）

John D'Errico
John D'Errico 2019-10-17
编辑：John D'Errico 2019-10-17
Too late for an answer for Angela, since I never saw the question a year ago. But a simple model for such a curve will not be too difficult. There would be at least 4 parameters to deal with.
  1. (A) The value of y at the bottom of one of the gullets. This is a simple translation in Y.
  2. (B) The slope of the teeth.
  3. (C) The location in x of ONE of the breaks between teeth. (Which one is irrelevant.)
  4. (D) The stride.
As I said, the curve fitting toolbox might fail. It is certainly true that any estimates of parameter confidence intervals would be garbage. But other tools that are more robust, perhaps fminsearch or better yet, a tool like GA would succeed, if the model truly fits the data. If the stride is not perfectly uniform however, expect problems in the fit. And you will need good starting values on any such problem.
The model might look like this:
y(x) = A + B*mod((x - C)/D,1)
Here is some fake data. I won't show you how I made it, as that would tell you the parameters.
plot(x,y,'o')
grid on
For kicks, I'll try out the CFTB on it.
ft = fittype('A + B*mod((x - C)/D,1)','indep','x','coef',{'A','B','C','D'});
Starting values can be difficult.
The minium of the gullets is easy to guess, just min(y).
If we have an estimate of the stride, then the slope can be estimated as max(y) - min(y))/stride.
You could use findpeaks to find the peak of each tooth. That would give you an estimat of one of the tooth breaks. You could also use it to find the stride as just the average difference between peak locations. Or you could just use the plot to graphically guess at a break location. LEts say I guess 25 as a break between teeth.
It looks like there are roughly 5 teeth in the curve. So an estimate of the stride would be simply (max(x) - min(x))/5.
strideguess = (max(x) - min(x))/5;
breakguess = 25;
abcdstart = [min(y),(max(y) - min(y))/strideguess,breakguess,strideguess];
mdl = fit(x,y,ft,'start',abcdstart)
mdl =
General model:
mdl(x) = A + B*mod((x - C)/D,1)
Coefficients (with 95% confidence bounds):
A = 2.613 (-1.168e+07, 1.168e+07)
B = 3.536 (2.861, 4.21)
C = 25.37 (-6.265e+07, 6.265e+07)
D = 18.96 (18.14, 19.79)
>> plot(mdl)
>> hold on
>> plot(x,y,'o')
untitled.jpgAs expected, it fits like crap. The confidence limits on the parameters are complete crap. Far better estimates for the starting parameters would be necessary. And even then, if your stride were not absolutely constant, forget it.
And that data is not very bad data. Pretty good, actually. Actually, perfectly splendid data I would claim.
  1 个评论
Fabio Freschi
Fabio Freschi 2019-10-17
编辑：Fabio Freschi 2019-10-17
My answer arrived later, I had my window open while testing the code and I didn't see your reply. In any case, results are similar even if I used lsqcurvefit instead of fit. I guess they share the same kernel.

请先登录，再进行评论。

Fabio Freschi
Fabio Freschi 2019-10-17
As usual, John is right. I tried to implement a strategy to fit sampled data with noise with a sawtooth waveform. I am pretty proud of the choice of the inital guess, but the fitting is anything but simple. Try yourself
clear variables, close all
%% params
% period
T = 1;
% peak
a = 2;
% delay
t0 = .3;
% noise
k = 0.3;
%% wave
% analytical expression of a sawtooth wave
y = @(x,t)2*x(2)*((t-x(3))/x(1)-floor(.5+(t-x(3))/x(1)));
%% samples
% params
x0 = [T a t0];
% sample data at random points in [0,2*T];
N = 50;
t0 = sort(2*T*rand(N,1));
% values with some random noise
y0 = y(x0,t0)+k*rand(size(t0));
%% the selection of initial guess is critical
% estimated peak
aStart = max(y0); % the max value available
% estimated period
idx = find(-diff(y0) > aStart); % I estimate the interval between discontinuities
TStart = mean(diff(t0(idx)));
% estimated delay
t0Start = t0(find(y0(2:end) > 0 & y0(1:end-1) < 0,1,'first')); % I estimate the firs zero crossing
% put all together
xStart = [TStart aStart t0Start];
% best fit
% optimization options (here you can play with thresholds, but it is not easy)
options = optimoptions('lsqcurvefit','FunctionTolerance',1e-20,'StepTolerance',1e-20);
% bounds
lb = [0 0 0]; % all params are poritive
ub = [];
[x,resnorm,residual,exitflag,output] = lsqcurvefit(y,xStart,t0,y0,lb,ub,options);
%% final plot
figure, hold on
plot(t0,y0,'o');
t = linspace(0,2*T,1000);
plot(t,y(x,t));
As you can see, different runs, with different random sampled data, provide different results. Maybe your data are more "regular" and the fitting becomes accetable. Let us know!
  1 个评论
Daniel M
Daniel M 2019-10-18
I had the idea that you could try to fit the parameters of a sine wave, and take the instantaneous phase of it to get a sawtooth. Here is some code I added to the end of your script. Sometimes it works well, sometimes it returns garbage.
% inst phase of a sine wave is a sawtooth
zfunc = @(z,t) z(1)*angle(hilbert(sin(2*pi*z(2)*t + z(3)))) +z(4);
% amplitude of phase, frequency, phase shift, vertical shift
amp = ((max(y0)-min(y0))/2)/pi; % angle() always between -pi and pi
z0 = [amp, 1./TStart, 0, 0];
lb = [0, 0, -pi -1];
ub = [4, 10, pi 1];
[z,resnormz,residualz,exitflagz,outputz] = lsqcurvefit(zfunc,z0,t0,y0,lb,ub,options);
figure, hold on
plot(t0,y0,'o');
plot(t,zfunc(z,t));

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Data Preprocessing 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by