Fitting a sinusoidal curve to a set of data points

I have a few sets of data points that I need to fit a sine curve of the form A*sin(B*X+C)+D on each data set. I don't know A, B, C and D (which defer for each data set) but I have access to cftool. I wrote a function (by generating the code from cftool) which doesn't work well for some data sets and I believe it's because of the wrong "start point". I wonder how to determine the startpoint for each fitting?! and if this is not possible, is there any other advice to fit the sine curve and get the coefficients (or generated data)?!

Thanking in advance,

Maryam

[Richard Garner]

I had this same problem awhile ago and I came up with a way to attack it which works extremely well. It involves a linear fit, and it's fast. The thing to do is first find the frequency of you data using an fft. Make sure you use a window function (e.g. Hann) otherwise frequency will be off. Thus, fft(data.*window), then get the frequency by centroiding around the place where it peaks up. Here you need a little discretion to decide on the window around the peak to choose.

Having found frequency, call it B, to use your notation, you now write you sinusoid (which will be a fit function as follows:

f(x;A,C,D)=A*sin(BX)*cos(C) + A*cos(BX)*sin(C)+D.

where A,C,and D are the free parameters to be determined. In this case you know B (very well, I think), so it's not a parameter. You now do a LINEAR fit of your data to this function to find parameters A*cos(C), A*sin(C), and D. This is done by forming the 2xN data matrix: D=[cos(Bx1),sin(Bx1),1 ; cos(Bx2),sin(Bx2),1 ; ... ; cos(BxN),sin(BxN),1] and the parameter matrix: p=[A*cos(C); B*sin(C); D]. Note that D*p=y, where y=[y1;y2;...;yN], and also the [x1;x2;...; xN] are the corresponding x's. Thus, to be specific, if you consider: p_fit=D\y, the p_fit will contain the best fit values of Acos(C), Asin(C), and D. From these you can get A and C. By the way, this produces much better results than doing a nonlinear fit by allowing B to be unknown parameter. This technique is much better with regard to noise in data. Furthermore, the noise in the parameters is Gaussian like if the noise in the data is also Gaussian. In the nonlinear fit this is not necessarily the case.

[John D'Errico]

As it turns out, Richard has posed what is ALMOST the complete and correct solution, but for one more step.

The trig identity almost makes it seem more complicated, but in reality, it reduces the problem in complexity, by effectively allowing you to solve for all but one of the parameters using a linear least squares.

The idea of what is called separable nonlinear least squares (sometimes also called a partitioned nonlinear least squares) is that some of these parameters are intrinsically linear. Thus given a value of the parameter B, we can solve for the other parameters using a LINEAR least squares solver. Since linear least squares does not need a starting value, it is not an iterative scheme that you need to deal with, linear least squares is a very nice thing.

The parameter B in this problem is the only one that is intrinsically nonlinear. So a separable nonlinear least squares tool (like my own fminspleas, as it is found on the file exchange) uses a nonlinear solver to iterate over those truly nonlinear parameters, and a linear solver to generate estimates of the other parameters.

So as I said, Richard has almost the right solution. Get a decent starting value for B using an fft. Then use a tool like fminbnd to solve for the value of B that minimizes the sum of squares, where the parameters A,C,D are all solved using the transformation and a linear least squares.

The virtue of the above scheme is it is very well behaved. It is FAST. It is robust to starting values for B, since the parameter search space is only a one dimensional search space.

[Maryam]

Thanks for your helps. The methods you've suggested all look sound but I am very new to Matlab and it took me a few hours to figure out what you were suggesting. Now, I have some questions that your further help would be greatly appreciated:

1- For what 'Star Strider' has suggested: I'm not sure how to proceed without any estimate of 'C'. Can I just put the first guess for C as, say 1?! Could you please give me a short example?!

2- I might need to be more accurate later, so regarding to 'John' and "Richard' methods: what would be L for Hann function?! and I also don't know how to do a linear fit and get those parameters. I think a short example could be very helpful. FYI_ I have different data sets, each has a 15101-velocity data measured at 25Hz. These are velocities from a wave flow. Some of those data sets are perfectly sinusoidal but some of them (from different run) are extremely noisy. I'm trying to replace those noises with a fitted sine curve. And I need to do so for all different sets of data in one code.

[John D'Errico]

There is a tool on the FEX, called fminspleas . It does all of the hard work for you.

http://www.mathworks.com/matlabcentral/fileexchange/10093-fminspleas

All you will need to do is define the function to be used. Here they would be something like...

funslist = {1, @(coef,xdata) cos(coef*xdata), @(coef,xdata) sin(coef*xdata)};

Fminspleas will return 4 coefficients, in two different return variables. Here, Bstart is the starting value for B. It is the ONLY one you need to provide.

[B,ILP] = fminspleas(funslist,Bstart,x,y);

D = ILP(1);
C = atan2(ILP(2),ILP(3));
A = ILP(2)/sin( C );

(I think I got those back-transformations right, but I wrote them down quickly, so worth checking.)

[Image Analyst]

Maryam's "Flag" moved here to be a comment:

Thanks for the detailed answer!

Answer 1

[Star Strider]

The easiest way to get ‘B*X’ is to use the fft function. You know ‘X’, so calculate to estimate ‘B’. (I leave the details of that to you.)

Find the mean of the curve to estimate ‘D’, and find the maximum of the absolute value of the function after subtracting the mean to estimate ‘A’.

This leaves ‘C’ without an initial estimate, but with three relatively close initial estimates for the others, ‘C’ should be relatively straightforward to estimate in your regression.

Good luck!

[Star Strider]

@Maryam —

With respect to ‘C’, since it is a phase term and likely between 0 and 2*pi, an initial parameter estimate with any number in that range, for instance 1 as you suggested, would probably work.

As for an example, an initial parameter estimates vector, ‘B0’ here, would be (in ABCD order) and denoting your data as ‘y’:

B0 = [max(abs(y-mean(y))), |B|, 1, mean(y)];

I left ‘B’ as is because you would get that estimate from the fft of your data.