**You are now following this question**

- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.

# Curve fitting to a sinusoidal function

1,041 views (last 30 days)

Show older comments

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.

##### 1 Comment

Douglas Lim
on 29 Feb 2016

### Accepted Answer

Star Strider
on 15 Mar 2014

Edited: Star Strider
on 15 Mar 2014

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.

##### 41 Comments

Patrick Fairclough
on 27 Jun 2015

Edited: Patrick Fairclough
on 27 Jun 2015

Star Strider
on 27 Jun 2015

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
on 5 Mar 2016

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

Star Strider
on 5 Mar 2016

Star Strider
on 11 Aug 2016

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
on 3 Nov 2016

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
on 3 Nov 2016

Thank you!

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

Carly Hauser
on 12 Jun 2017

Star Strider
on 12 Jun 2017

As always, my pleasure!

Kunal Kumar Saraf
on 4 Aug 2017

Ahmad Suliman
on 7 Sep 2017

Edward Blocker
on 17 Jun 2019

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

Star Strider
on 17 Jun 2019

@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.
on 29 Feb 2020

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
on 29 Feb 2020

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
on 23 Oct 2020

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
on 23 Oct 2020

Quentin COSSON
on 3 Aug 2021

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
on 3 Aug 2021

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
on 12 Oct 2021

Cool answer! It helps me a lot! Thanks!

Star Strider
on 12 Oct 2021

Youngju Kim — My pleasure!

A vote is always appreciated!

.

zizo hiro
on 24 Oct 2021

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
on 24 Oct 2021

Maya Eyal
on 16 Jan 2022

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
on 16 Jan 2022

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
on 25 Apr 2022

Thank you for the incredible code.

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

Star Strider
on 25 Apr 2022

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)’.

Ashfaq Ahmed
on 7 Mar 2023

Edited: Ashfaq Ahmed
on 7 Mar 2023

@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
on 7 Mar 2023

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
on 7 Mar 2023

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
on 7 Mar 2023

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
on 7 Mar 2023

It worked. Thank you.

Star Strider
on 7 Mar 2023

My pleasure!

Hector Camilo Clavijo Zarate
on 1 May 2023

Star Strider
on 1 May 2023

It is essentially just this:

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

.

Hector Camilo Clavijo Zarate
on 2 May 2023

@Star StriderWonderful! thank you for all your help

Star Strider
on 2 May 2023

My pleasure!

A Vote would be appreciated!

### More Answers (2)

Jos (10584)
on 14 Mar 2014

##### 1 Comment

Dejan
on 15 Mar 2014

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

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!**An Error Occurred**

Unable to complete the action because of changes made to the page. Reload the page to see its updated state.

Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list

How to Get Best Site Performance

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

Americas

- América Latina (Español)
- Canada (English)
- United States (English)

Europe

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)

Asia Pacific

- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)