Hi Star Strider, I have further question on this topic. May I know how to contact you for sending you the problem.
Curve fitting to a sinusoidal function
554 次查看（过去 30 天）
显示 更早的评论
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.
采纳的回答
Star Strider
2014-3-15
编辑：Star Strider
2014-3-15
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 个评论
Patrick Fairclough
2015-6-27
编辑：Patrick Fairclough
2015-6-27
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
2015-6-27
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.
Star Strider
2016-3-5
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
2016-8-11
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
2016-8-11
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
2016-11-3
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.
Carly Hauser
2017-6-12
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.
Kunal Kumar Saraf
2017-8-4
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
2017-9-7
Thank you for the code. However, I am struggling what x is in the line where you are finding zero crossings. Thanks, Ahmad
Star Strider
2019-6-17
@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.
2020-2-29
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
2020-2-29
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
2020-10-23
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
2020-10-23
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
2021-8-3
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
2021-8-3
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.
.
zizo hiro
2021-10-24
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
2021-10-24
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
2022-1-16
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
2022-1-16
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
2022-4-25
Thank you for the incredible code.
How do I write such a code when all 'y' values are positive?
Star Strider
2022-4-25
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
2022-6-10
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
2023-3-7
编辑：Ashfaq Ahmed
2023-3-7
@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
2023-3-7
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
2023-3-7
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
2023-3-7
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.
.
Hector Camilo Clavijo Zarate
2023-5-1
@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
2023-5-1
It is essentially just this:
The β values are the elements of the parameter vector to be estimated.
.
更多回答（2 个）
Jos (10584)
2014-3-14
1 个评论
Dejan
2014-3-15
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
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Linear and Nonlinear Regression 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!发生错误
由于页面发生更改，无法完成操作。请重新加载页面以查看其更新后的状态。
Translated by 
您也可以从以下列表中选择网站：
如何获得最佳网站性能
选择中国网站（中文或英文）以获得最佳网站性能。其他 MathWorks 国家/地区网站并未针对您所在位置的访问进行优化。
美洲
- América Latina (Español)
- Canada (English)
- United States (English)
欧洲
- 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)
亚太
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)