Curve fitting to a sinusoidal function

I have a series of data points that are governed by a sinusoidal function.

I want to fit, plot and generate a sinusoidal function to these data points.

I do not wish to fit an nth degree polynomial to this no matter how close it is to the sinusoidal function.

I understand that there is no standard tool in the toolbox that does this. I have looked at numerous and plenty of old threads and internet posts. None of the answers help me.

Please take into account that I am new to Matlab and can only curve fit very basic data points.

What I therefore need is an exact and step by step guide in how to fit a sine curve to data points.

Please assist.

[Douglas Lim]

Hi Star Strider, I have further question on this topic. May I know how to contact you for sending you the problem.

Answer 1

[Star Strider]

Here’s my suggested solution, using only core MATLAB functions:

yu = max(y);
yl = min(y);
yr = (yu-yl);                               % Range of ‘y’
yz = y-yu+(yr/2);
zx = x(yz .* circshift(yz,[0 1]) <= 0);     % Find zero-crossings
per = 2*mean(diff(zx));                     % Estimate period
ym = mean(y);                               % Estimate offset

fit = @(b,x)  b(1).*(sin(2*pi*x./b(2) + 2*pi/b(3))) + b(4);    % Function to fit
fcn = @(b) sum((fit(b,x) - y).^2);                              % Least-Squares cost function
s = fminsearch(fcn, [yr;  per;  -1;  ym])                       % Minimise Least-Squares

xp = linspace(min(x),max(x));

figure(1)
plot(x,y,'b',  xp,fit(s,xp), 'r')
grid

The elements of output parameter vector, s ( b in the function ) are:

s(1): sine wave amplitude (in units of y)

s(2): period (in units of x)

s(3): phase (phase is s(2)/(2*s(3)) in units of x)

s(4): offset (in units of y)

It provides a good fit.

[Patrick Fairclough]

A very elegant solution, I like how you made allowance for the original signal not crossing zero and then shifting it to find the zero. Would this count data 2, 0, -2 as two zero crossings? (<=0)

[Star Strider]

Thank you, Professor Fairclough. I appreciate your compliment.

Using this construction:

t = 10:12;
x = [2 0 -2]';
zx = t(x .* circshift(x,[0 1]) <= 0)

it returns:

zx =
    11

because:

L = (x .* circshift(x,[0 1]) <= 0)

returns:

L =
     0
     1
     0

So it returns only one zero-crossing.

In this instance (and others in which I’ve used it), it provides an initial estimate the zero-crossings in the context of a regression (such as this), or to define a range of values for a linear interpolation. In that context, close enough is good enough, since it is an initial estimate of a more precise value.

Defining a ‘zero-crossing’ using this method is by definition inexact because it ‘chooses’ (pardon the anthropomorphism) the index on one side or the other of the actual zero-crossing if the function is not exactly zero. (That depends on whether the circshift call uses +1 or -1.) Here it was equal to zero, and the logic detected it correctly.

[Katie Cannon]

How can I then find out what the equation of my sine wave is?

[Star Strider]

You need an independent variable vector (here ‘x’) and a signal vector (your sine wave, here ‘y’) of the same lengths. Subtract the mean from the sine wave if it is not already close to zero so it has zero-crossings, then use my code.

[Karl M]

Hi, I was wondering what '-1' in s = fminsearch(fcn, [yr; per; -1; ym]) stand for? In my case other numbers showed to be a better fit. Thank you!

[Star Strider]

My pleasure!

The ‘-1’ was part of the phase term, and that choice of initial parameter estimates made the function converge. (Nonlinear parameter estimation routines can be extremely sensitive to the initial parameter estimates, so experimenting to see what works is necessary. Here, ‘-1’ for the initial phase period worked best.) My code is designed to derive its initial parameter estimates from the data, to make it as robust as possible with data similar to the data here. It will need tweaking in some situations to produce the best fit.

[Johannes Thewes]

Great solution! For my use case, I included a phase estimation that worked well for me:

phase = mod(-per/zx(1),per)

@Star Strider: If you have no copyright objections, I would like to publish my adapted function on the File Exchange (mentioning this thread) and put a link here.

[Star Strider]

Thank you!

I have no objections at all. I do appreciate the credit.

[Carly Hauser]

Thank you so much for sharing your code! It does exactly what I need it to do for my team's research project. I couldn't find any helpful code anywhere else on the internet.

[Star Strider]

As always, my pleasure!

[Kunal Kumar Saraf]

Hello Sir, I have a small doubt in your solution. As I understood it b3 is the variable for Phase. And when you are solving for fminsearch, why you have taken that variable as (-1). can you plz throw some light on it?

[Ahmad Suliman]

Thank you for the code. However, I am struggling what x is in the line where you are finding zero crossings. Thanks, Ahmad

[Baha411]

This code is not working man, e.g. for, x = 2* pi*rand(100,1); x = sortrows(x); y = 9 + 7*sin(2*x + 4*pi) + randn(100,1);

Am I missing something on this?

[Edward Blocker]

How could I edit this to fit the same function, but with sine^2 ??

[Star Strider]

@Edward Blocker —

Probably:

fit = @(b,x)  b(1).*(sin(2*pi*x./b(2) + 2*pi/b(3))).^2 + b(4);    % Function to fit

Note: UNTESTED with respect to .

[M.A.G.]

How can I interpret yz (line 4)

[M.A.G.]

also, reading circshift documentation:

The default behavior of circshift(A,K) where K is a scalar changed in R2016b. To preserve the behavior of R2016a and previous releases, usecircshift(A,K,1). This syntax specifies 1 as the dimension to operate along.

should this allter the output in this case?

[Star Strider]

Interesting that you brought that up.

I recently updated that in another Answer, adding a reference to the ‘zci’ function that I wrote many years ago:

x = linspace(0,5*pi);
y = 100*sin(x)+700;

yu = max(y);
yl = min(y);
yr = (yu-yl);                                                   % Range of ‘y’
yz = y-yu+(yr/2);
zci = @(v) find(v(:).*circshift(v(:), 1, 1) <= 0);              % Returns Approximate Zero-Crossing Indices Of Argument Vector (>= R2016b)
zx = x(zci(yz));                                                % Find zero-crossings
per = 2*mean(diff(zx));                                         % Estimate period
ym = mean(y);                                                   % Estimate offset
fit = @(b,x)  b(1).*(sin(2*pi*x.*b(2) + 2*pi*b(3))) + b(4);     % Function to fit
fcn = @(b) sum((fit(b,x) - y).^2);                              % Least-Squares cost function
s = fminsearch(fcn, [yr;  1/per;  -1;  ym])                     % Minimise Least-Squares
xp = linspace(min(x),max(x));
figure(1)
plot(x,y,'b',  xp,fit(s,xp), 'r')
grid

This also changed the period representations to frequencies. Both versions will work.

[Steven]

Hi Star

Thank you for this function.

1.It can fairly fit a good sin fucntion. How can we find the frequency from this fitted function? (any quick simple way without fft, as that didn't work)

Here is my data:

y = [142 137 147 131 103 107 145 157 151 139 138];

x = 0:1:length(y)-1;

Thank you

[Star Strider]

In this function, ‘b(2)’ is the frequency (in units of cycles/(independent variable unit), with that defined in the independent variable vector).

[Quentin COSSON]

This code works great.

I've used it for to fit data during my internship.

I've included the link to your code in mine for my supervisor to have the reference.

[Star Strider]

Quentin COSSON —

Thank you!

One update since I originally posted this nearly 7½ years ago:

zxi = find(diff(sign(yz)));
zx = x(zxi);

or, using the anonymous function implementation:

zci = @(v) find(diff(sign(v)));
zx = x(zci(yz));

This is more robust, and does not have problems with ‘end effects’ when there is a ‘false’ zero-crossing at the end of the vector due to the ‘wrap-around’ effect circshift introduces.

.

[Youngju Kim]

Cool answer! It helps me a lot! Thanks!

[Star Strider]

Youngju Kim — My pleasure!

A vote is always appreciated!

.

[zizo hiro]

The function used to fit, I would like to know where did you get the 2*pi/b(3) from? I expected according to the general form of a sin wave to have something like

b(1).*(sin(2*pi*x./b(2) + b(3))) + b(4);

with b(3) being the phase shift

[Star Strider]

zizo hiro — The parameters that are arguments to the sin call are in terms of periods, not frequencies. That’s just the way I defined them because of the way I defined the initial parameter estimates.

[Maya Eyal]

NOTICE: About the sine parameter's units, notice that the phase can sometimes be unitless (or in Radians). For example if Y is describing an angular velocity.

@Star Strider, you may want to edit your (well described) answer...

[Star Strider]

Maya Eyal — I appreciate your compliment!

The arguments to transcendental functions are always unitless (dimensionless), so for example angular frequency is multiplied by time rendering it unitless, and the phase is likewise unitless. The angular frequency fundamental units (radians/second, cycles/fortnight, etc.) and phase units (radians, cycles) must always be the same, and must not change throughout the system being calculated.

.

[Negin Rahmati]

Thank you for the incredible code.

How do I write such a code when all 'y' values are positive?

[Star Strider]

Negin Rahmati —

The ‘y’ offset is initially estimated by:

ym = mean(y); % Estimate offset

and that is used throughout the code to estimate the relevant parameters, including the zero-crossing calculations. It also becomes the initial parameter estimate for the offset parameter ‘b(4)’.

[Bastiaan]

I noticed that the initial estimate of b(1) is off by a factor of 2. It think it should be equal to yr/2 instad of yr:

s = fminsearch(fcn, [yr/2;  1/per;  -1;  ym]) 

[Ashfaq Ahmed]

@Star Strider this is actually a great code that you wrote. Helped me a lot doing a few projects. However, in some cases it has limitations and I am trying to understand why. For example, this -

x = -3*pi:0.1:3*pi;
y = rand(1,size(x,2))*0.5+sin(x);
yu = max(y);
yl = min(y);
yr = (yu-yl);                                                   % Range of ‘y’
yz = y-yu+(yr/2);
zci = @(v) find(v(:).*circshift(v(:), 1, 1) <= 0);              % Returns Approximate Zero-Crossing Indices Of Argument Vector (>= R2016b)
zx = x(zci(yz));                                                % Find zero-crossings
per = 2*mean(diff(zx));                                         % Estimate period
ym = mean(y);                                                   % Estimate offset
fit = @(b,x)  b(1).*(sin(2*pi*x.*b(2) + 2*pi*b(3))) + b(4);     % Function to fit
fcn = @(b) sum((fit(b,x) - y).^2);                              % Least-Squares cost function
s = fminsearch(fcn, [yr;  1/per;  -1;  ym])                     % Minimise Least-Squares

s = 4×1

0.0501 0.4469 -1.5754 0.2469

xp = linspace(min(x),max(x));
figure(1)
plot(x,y,'b',  xp,fit(s,xp), 'r')
grid

[Star Strider]

Thank you!

My code will have problems with noisy, unfiltered data.

x = -3*pi:0.1:3*pi;
y = rand(1,size(x,2))*0.5+sin(x);
yf = lowpass(y, 0.01, 0.1);                                     % Lowpass Filter
yu = max(yf);
yl = min(yf);
yr = (yu-yl);                                                   % Range of ‘y’
yz = yf-yu+(yr/2);
zci = @(v) find(diff(sign(v)));                                 % Returns Approximate Zero-Crossing Indices Of Argument Vector (>= R2016b)
zx = x(zci(yz));                                                % Find zero-crossings
per = 2*mean(diff(zx));                                         % Estimate period
ym = mean(y);                                                   % Estimate offset
fit = @(b,x)  b(1).*(sin(2*pi*x.*b(2) + 2*pi*b(3))) + b(4);     % Function to fit
fcn = @(b) sum((fit(b,x) - yf).^2);                             % Least-Squares cost function
s = fminsearch(fcn, [yr;  1/per;  -1;  ym])                     % Minimise Least-Squares

s = 4×1

1.0102 0.1584 -1.0002 0.2417

xp = linspace(min(x),max(x));
figure(1)
plot(x,y,'b',  xp,fit(s,xp), 'r')
grid

Once the noise is filtered out so that the code can accurately estimate the parameters of the filtered signal, it works.

I defined the fitler parameters empirically by looking at the data. In other circumstances, it might be necessary to do a Fourier transform of ‘y’ first in order to determine the optimal filter cutoff frequency. Any reasonably robust filter, including the smoothdata function and its friends, will likely work.

I also updated the ‘zci’ function here to be more robust.

.

[Ashfaq Ahmed]

Yes, this works MUCH better now! I have one question out of curiosity: how can I save the values for the 's' if I run this code in a for loop? like this -

for i = 1:10
    figure(1),
    x = 1:size(X{i},1);
    y = mean(SST(X{i}(1):X{i}(end),1:end),2,'omitnan')';


    yf = lowpass(y, 0.01, 0.1);                                     % Lowpass Filter
    yu = max(yf);
    yl = min(yf);
    yr = (yu-yl);                                                   % Range of ‘y’
    yz = yf-yu+(yr/2);
    zci = @(v) find(diff(sign(v)));                                 % Returns Approximate Zero-Crossing Indices Of Argument Vector (>= R2016b)
    zx = x(zci(yz));                                                % Find zero-crossings
    per = 2*mean(diff(zx));                                         % Estimate period
    ym = mean(y);                                                   % Estimate offset
    fit = @(b,x)  b(1).*(sin(2*pi*x.*b(2) + 2*pi*b(3))) + b(4);     % Function to fit
    fcn = @(b) sum((fit(b,x) - yf).^2);                             % Least-Squares cost function
    s{i} = fminsearch(fcn, [yr;  1/per;  -1;  ym])                     % Minimise Least-Squares
    xp = linspace(min(x),max(x));
end

It says

Unable to perform assignment because brace indexing is not supported for variables of
this type.

[Star Strider]

I don’t have the data, so i can’t run this to test it.

First optionally preallocate cell array ‘sc’ before the loop:

sc = cell(10,1);

then after determing ‘s’, save it to elements of ‘sc’:

sc{i} = s;

That should work.

.

[Ashfaq Ahmed]

It worked. Thank you.

[Star Strider]

My pleasure!

[Hector Camilo Clavijo Zarate]

@Star Strider Hi! nice to meet you. I would like to ask you if you could please share the sine equation with the parameters that the code provides? I used your code ( I will give you credits on my work) and it's amazing how it approaches the data curve but I have been unable to gather the equation itself, thank you!

[Star Strider]

@Hector Camilo Clavijo Zarate — The equation is the ‘fit’ anonymous function.

It is essentially just this:

The β values are the elements of the parameter vector to be estimated.

.

[Hector Camilo Clavijo Zarate]

@Star StriderWonderful! thank you for all your help

[Star Strider]

My pleasure!

A Vote would be appreciated!

Answer 2

[Chad Greene]

I turned this into a function called sinefit for climatological data with a periodicity of 1 year.

Answer 3

[Jos (10584)]

Take a look here:

http://www.mathworks.com/matlabcentral/answers/36999-how-do-i-regression-fit-a-sinwave-to-a-dataset

[Dejan]

Im relatively new to Matlab and whilst the link provided I kinda get, its not an exact step by step guide on how to fit a Sine wave.

Here are my data points:

x = [ 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 ]

y = [ 16.5 14.32 11.58 10.017 9.629 10.2 12.16 15.08 16.97 16.75 14.331 11.508 10.013 9.617 10.22 12.15 15.304 17.38 16.853 ]

I need a Sine wave to fit that data and its governing equation