MATLAB Answers

How to fit this signal to x axis at y=0 using polyfit/polyval?

5 views (last 30 days)
Ebrahim Abdullah Ahmed Albabakri
Commented: John D'Errico on 30 Nov 2020
I want to fit this signal to the x axis because it looks more shifted upwards, any help will be much appreciated.
x = Graph1(:, 1);
y = Graph1(:, 2);
plot(x, y);


Sign in to comment.

Answers (1)

John D'Errico
John D'Errico on 29 Nov 2020
Edited: John D'Errico on 29 Nov 2020
This is just another case of someone not understanding how polynomials work or behave.
Seriously, to fit a polynomial to that data would be an unjust thing to do to a poor, unsuspecting polynomial. Why? Well, yes, you can force fit a polynomial through the origin. The data is massively noisy but you could still do so.
You would need to discard the data below x == 0 though, since it appears your function is essentially constant below zero.
As well, above 2 or 3 for x, the function again appears to settle down to a constant value. This is impossible for any polynomial function (except for the trivial case of a constant itself.) All polynomials, as x grows large in magnitude, ALWAYS MUST grow to either +inf or -inf.
As such, your data shows strongly non-polynomial behavior.
If you were willing to discard all data below zero, as well as all data above say x=2.5, could you fit what remains with a polynomial that is forced to pass through zero? Well, yes.
Polyfit does not have this capability. And you would NOT want to use a remotely high order polynomial. Your data is wildly too noisy to support anything of any significant order. (Perhaps a cubic is as high as I would go.) And, since I lack your data, I can't do too much. A picture of numbers is nice, but not terribly useful. So I;ll make up some garbage data. Actually, what I have here is a bit better than what I see in the picture.
x = rand(100,1)*2.5;
y = sin(x) + randn(size(x))/3;
So, pretty noisy data. You would agree the function should pass through the origin, in this case, only because we know the underlying function does exactly that. Now, fit that with a cubic polynomial, passing through the origin. This takes only one line of code.
coef = (x.^[3:-1:1])\y;
I could also have trivially used the curve fitting toolbox, or my own polyfitn from the file exchange. But one line of code that requires no other tools is hard to beat.
The constant term in that model is implicitly zero. So that I can now use polyval to evaluate it, I'll append the constant term to the coefficients.
coef = [coef;0];
hold on
fplot(@(x) polyval(coef,x),[0,2.5])
legend('data','Cubic fit through the origin')
A higher order model would be wholly unjustified for data this noisy, but the cubic seems reasonable. Could I have added a up to maybe 5th order term? Sigh, perhaps. At some point the extra terms will just start chasing noise in the data, and given this much noise, that will happen quite quickly.
Do NOT try to use that polynomial to predict outside the support of the data. It will do obscene things, as that is just the nature of the polynomic beast.


Show 1 older comment
John D'Errico
John D'Errico on 30 Nov 2020
I'm at a loss here, because what you ask seems to make no sense.
  1. Away from polyfit, to me, means you don't want to use polyfit. Yet, I showed how to fit the data without using polyfit. So what is the problem? Yes, I used a polynomial model, but that is what you asked to fit. No?
  2. You just want to use the data greater than zero, so then use it. Create a new variable, where you have excluded all data where x<0. Are you asking how to do that?
ind = T>=0;
That = T(ind);
Ahat = A(ind);
You won't be happy, since as I said, a polynomial will NOT fit that data. You need to read my answer. A polynomial does NOT have the behavior your data shows.
So there is now no data below x == 0. A polynomial is still totally inappropriate to fit this data. Must I repeat what I said?
A polynomial has the property that as the variable grows large, it will ALWAYS approach either +inf or -inf. As long as the highest order non-zero coefficient is non-zero, this will be true. And if the highest order coeffiicent is zero, then the polynomial is the trivial one - a constant. So you have two simple choices in the world of polynomials. Either discard all of your data above roughly a time of 2e-6 (as I did for the data below zero) or do not fit a polynomial to your data.
In my eyes, they seem to be simple choices.
Perhaps you mean to ask how to just fit the data. But that is a meaningless question. You cannot just "fit" data. You must fit some model to your data. And until you pose some model for this data, you have no ability to just fit it. To fit implies a model. You have no model, at least, you have posed no such thing. And your data seems to lack sufficient signal to noise ratio to justify more than a LOW order polynomial model. What I see below a time of 2e-6 is a highly oscillatory signal.
I could fit your data using a least squares spline model, at least, if I use a model that I force to be monontonic increasing, that goes through zero, and has a zero slope as x grows large. These constraints are sufficiently constraining that the model looks pretty smooth. It must, because there is insufficient flexilbility for anything else to happen. The fact remains though, I will have posed a model to fit the data. So, here, using my SLM toolbox (free download from the file exchange.)
slm = slmengine(That,Ahat,'knots',linspace(0,max(That),10),'leftvalue',0,'rightslope',0,'increasing','on');
The data overwhelms the curve, so here I have plotted only the curve:
Again, you cannot just fit some data, without specifying a model to use for the fit. In the above, I used a regression spline that allows constraints on the curve, which serve to regularize the model.
Or, I might have chosen a nonlinear model. Again, to fit without a model is a meaningless operation. I can try using the curve fitting toolbox instead.
In this model, the curve is a simple exponential rise to an asymptote at y==a, with a rate constant of b. Here, the slope of the curve at T==0 is a/b.
mdl = fittype('a*(1-exp(-t/b))','indep','t')
mdl =
General model:
mdl(a,b,t) = a*(1-exp(-t/b))
fittedmdl = fit(That,Ahat,mdl,'start',[.06 1e-6])
fittedmdl =
General model:
fittedmdl(t) = a*(1-exp(-t/b))
Coefficients (with 95% confidence bounds):
a = 0.05844 (0.05757, 0.05931)
b = 4.646e-07 (4.174e-07, 5.117e-07)
I might have chosen a subtly different sigmoidal shape, to predict an exponential rise to an asymptote, where the curve is forced to pass through the origin at time T=0. In this model, c is a variation of a rate parameter, and a is the value of the asymptote. In this model, the parameter b alows the model to infer the slope at T==0.
mdl = fittype('a.*(1./(1 + b*exp(-x/c))*2 - 1)','indep','x')
mdl =
General model:
mdl(a,b,c,x) = a.*(1./(1 + b*exp(-x/c))*2 - 1)
>> fittedmdl = fit(That,Ahat,mdl,'start',[.06 1 1e-6])
fittedmdl =
General model:
fittedmdl(x) = a.*(1./(1 + b*exp(-x/c))*2 - 1)
Coefficients (with 95% confidence bounds):
a = 0.0587 (0.05784, 0.05957)
b = 1.002 (0.8621, 1.141)
c = 3.617e-07 (3.139e-07, 4.095e-07)
In that model, a is the asymptote at infinite time, and c is a rate parameter that measures how quickly the model approaches the asymoptote. So the units on c are time. b is a shape parameter, that controls the slope of the rise at T==0.
In any case, you cannot just fit something without a model. And you need to understand the model you will use.
Ebrahim Abdullah Ahmed Albabakri
I don't want any change of values in the signal, I just want to draw the x-axis at y==0 instead of at the buttom of the graph. I heared it could be done by polyfit but if you have any other way it is okay.
John D'Errico
John D'Errico on 30 Nov 2020
SigH. Is that all you want? How much time have I wasted trying to answer the question you asked, because you explicitly asked how to fit your signal, using polyfit. No, polyfit or polyval would seem to have nothing to do with that. I have no clue why you even thought they would.
Your alternatives there would seem to be one of: 'bottom', 'top', 'origin', with the default as 'bottom'.

Sign in to comment.




Community Treasure Hunt

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

Start Hunting!

Translated by