File Exchange

gaussian curve fit

version 1.0 (1.71 KB) by

gaussian curve fit

3.65517
26 Ratings

Updated

[sigma,mu,A]=mygaussfit(x,y)
[sigma,mu,A]=mygaussfit(x,y,h)

this function is doing fit to the function
y=A * exp( -(x-mu)^2 / (2*sigma^2) )

the fitting is been done by a polyfit
the lan of the data.

h is the threshold which is the fraction
from the maximum y height that the data
is been taken from.
h should be a number between 0-1.
if h have not been taken it is set to be 0.2
as default.

Carlos Macias

Ethan Montag

Ethan Montag (view profile)

I found a paper that describes how to get a higher resolution estimate of a peak frequency in a 1D FFT.using three points centered on a nominal peak. The paper described how to estimate the new peak frequency but didn't say how to get the new amplitude.

This function can do the fit for you. It gives you the frequency as mu and the amplitude as A when you input the three x values (frequencies, evenly spaced) and the y values (amplitudes, with the middle as the nominal peak). The polyfit fn generates a warning but you can turn it off with warning('off','MATLAB:polyfit:RepeatedPointsOrRescale'). Thanks!

michael scheinfeild

michael scheinfeild (view profile)

i want to thank you great fit
if you want c++ implement just solve the fit with qr

Mat src1(len,3,CV_32FC1);

Mat src2(len, 1, CV_32FC1);

for (int i = 0; i < len;i++)
{
src2.at<float>(i) = float(ylog[i]);
}

for (int i = 0; i < len; i++)
{
src1.at<float>(i, 0) = float(xlog[i] * xlog[i]);
src1.at<float>(i, 1) = float(xlog[i]);

src1.at<float>(i, 2) = float(1);
}

bool ret = cv::solve ( src1,
src2,
dst,
flags ) ;

Yu-jiao Shen

Yu-jiao Shen (view profile)

I am badly in need of this file.Thank you very much for your help.

Dan

Dan (view profile)

Felipe Santibañez

Achut

Georges

Georges (view profile)

This file poorly fits a gaussian (and has a much higher average rating) as compared to others here...
The test main file is wrong: the "rand" noise on the gaussian biases the measurements by 0.5, randn should be used instead.

Same main file used for comparison with another gaussian fit (gaussfit.m):

%% data
x=1:100;
sigma=15; mu=40; A=3;
plot(normpdf(x,mu,sigma))
y=A*exp(-(x-mu).^2/(2*sigma^2))+randn(size(x))*0.5;
hold all;
plot(x,y/(sum(y)),'.');

%% fitting
[sigmaNew,muNew,Anew]=mygaussfit(x,y);
y2=Anew*exp(-(x-muNew).^2/(2*sigmaNew^2));
plot(x,y2/(sum(y2)));

%% bis
[sigmaNew2,muNew2]=gaussfit(x,y/sum(y));
plot(normpdf(x,muNew2,sigmaNew2));
legend('truth','measurements','mygaussfit','gaussianfit')

Qing

chile1987

chile1987 (view profile)

Great idea,

though I wonder when you do the retransformation from second order polynomial (a0 + a1*x + a2*x^2)to logarithmic gaussian (log(A) - x^2/2*sigma^2 + x*mu/2*sigma^2 - mu^2/2*sigma^2)(line 46 in mygaussfit) shouldn't you calculate the mean from p(2) via mu = A1/-A2 or equivalently mu = A1*2*sqrt(-1/2*A2) instead of mu = A1*2*sigma^2 ?

I might be totally wrong or missed sth, just a quick idea.

Manuel A. Diaz

Omer

Namra Aftab

Namra Aftab (view profile)

hi...i have a sequence of 40 frames and i want to plot gaussian distribution of, say, first pixel of all 40 frames.the result does not look like a gaussian at all.what should i correct?

hello every one any one can tell me about gussian curve fitting back groung why we use this instead of other curve fitting method what is the benift of this from other curve fitting method.If any one have some good data regarding gussian curve fitting kindly inform me.

Avi

GREAT

Kiran

Kiran (view profile)

hello...could someone please explain to me how the approximation from polynomial back to gaussian is done in the code? why are the coefficients equated in the way that they are below :

sigma=sqrt(-1/(2*A2));
mu=A1*sigma^2;
A=exp(A0+mu^2/(2*sigma^2));

Thank you very much

Jan

Jan (view profile)

Where do you get he formulas from? and where do you use h for? Can it also without h?

David Holz

David Holz (view profile)

Bob Marley is absolutely correct. The use of a for loop is unnecessary and slows the function down by ~20%.

If you are confident in your inputs you can also remove the warning lines #74-85 of polyfit to knock out 87% of the total previous computation time.

Ryan

Ryan (view profile)

I disagree with Matheca. The function is intended to fit a general gaussian, not necessarily a probability distribution function. The equation is correct.

However, the user should be aware that removing data points in a deterministic manner (i.e. by thresholding) definitely skews the resulting fit.

Rather than fitting to the whole series with negatives removed, try finding the largest contiguous positive subset of the original data series and fitting to that. This method won't work when the noise amplitude is greater than the distribution amplitude, but in most cases it will give you a better fit.

Even better yet: if accuracy is more important than computation speed, use fmincon with a least-squares difference cost function:

p = fmincon(@(p) sum((y-(p(1)*exp(-(x-p(2)).^2/2/p(3)^2))).^2),...
[max(y) mean(x) 1],[],[],[],[],[0 -inf 0],[2*max(y) inf inf])
A=p(1);
mu=p(2);
sigma=p(3);

Ryan

Michael Jordan

an ‰Ê

an ‰Ê (view profile)

thanks a lot

Matheca ProbStock

The formula used for a Gaussian pdf is wrong. pdf(x)=(A/sqrt(2*sigma)) * exp( -(x-mu)^2 / (2*sigma^2) )...should be used.

Changlong Jin

The value of parameter h severely influence the result, in last comment, I use the default value, the fitting result is not correct, it looks more better when h = 0.1, how to solve this automatically?

Changlong Jin

there is a problem,
When I fit a data, for example,
y = [0.0651 0.0548 0.0461 0.0686 0.1268 0.2266 0.2292 0.1187 0.0299 0.0146 0.0092 0.0048 0.0032 0.0024];
it gives out a result like this:
yout = [0.0470 0.0594 0.0743 0.0918 0.1120 0.1352 0.1611 0.1897 0.2208 0.2538 0.2884 0.3238 0.3592 0.3937];
it is not a good fitting, how to solve this problem?

David last name

great, thanks!

Marie-Eve Gagne

Thank you very much!
This is exactly what I was looking for!

Dinh Vo

Bob Marley

The code could be written in a more efficient manner -- i.e., using matlab syntax instead of the 'for' loop. Something like:

%% threshold
if nargin==2, h=0.2; end

%% cutting
indx = (y > max(y)*h);

%% fitting
p=polyfit(x(indx), log(y(indx)), 2);
sigma=sqrt(-1/(2*p(1)));
mu=p(2)*sigma^2;
A=exp(p(3)+mu^2/(2*sigma^2));

Marinna M

Thanks, this was a great help.

ll ss

very helpful but I don't think you should use logx for calculation and you should just do the polyfit with ylog and x.

Gustavo P

Very smart idea

m m