# Sine curve fitting in MATLAB

38 views (last 30 days)

Show older comments

y = [ 0.0060 ;0.0077 ;0.0058];

x = 0:0.1:3;

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

##### 1 Comment

Cris LaPierre
on 24 Mar 2023

### Accepted Answer

Cris LaPierre
on 24 Mar 2023

Edited: Cris LaPierre
on 24 Mar 2023

For fminsearch to work correctly, your function must return a scalar. From the fminsearch documentation:

"fun is a function that accepts a vector or array x and returns a real scalar f (the objective function evaluated at x)."

Your y vector is not the same length as x, so the result of (fit(b,x) - y).^2 is a 3x189 array. Since this is a 2D array, sum will sum the values along the first dimension (columns) resulting in a 1x189 vector.

Instead, define an objective function that captures the overall goodness of fit in a single scalar value. The traditional way to do this would be to define the actual Y value for each x position.

### More Answers (1)

Star Strider
on 24 Mar 2023

@M — I was in the process of answering this when Win 11 crashed. Again. For the fifth time in two days, and three times when I was in the middle of doing something on Answers. My hatred of Windows and its infernal stability issues knows no bounds at this point! Break up Micro$oft!

For what it’s worth (which isn’t much at this point) —

A = readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1334739/test4.txt');

A = fillmissing(A,'linear');

L = size(A,1);

Fs = 1000;

t = linspace(0, L-1, L).'/Fs;

x = t;

for k = 1:size(A,2)

y = A(:,k);

yu = max(y);

yl = min(y);

yr = (yu-yl); % Range of ‘y’

yz = y-yu+(yr/2);

zci = @(v) find(diff(sign(v-mean(v)))); % Returns Approximate Zero-Crossing Indices Of Argument Vector

zt = x(zci(y));

per = 2*mean(diff(zt)); % Estimate period

ym = mean(y); % Estimate offset

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

fcn = @(b) norm(fit(b,x) - y); % Least-Squares cost function

[s,nmrs] = fminsearch(fcn, [yr; -10; per; -1; ym]) % Minimise Least-Squares

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

figure

plot(x,y,':b', 'LineWidth',1.5)

hold on

plot(xp,fit(s,xp), '-r')

hold off

grid

xlabel('Time')

ylabel('Amplitude')

legend('Original Data', 'Fitted Curve')

text(0.3*max(xlim),min(ylim)+0.05*diff(ylim), sprintf('$y = %.3f\\cdot e^{%.3f\\cdot x}\\cdot sin(2\\pi\\cdot x\\cdot %.0f%.3f) %+.2f$', [s(1:2); 1./s(3:4); s(5)]), 'Interpreter','latex')

title("Column "+k)

end

The fitted functions are not perfect, however the data themselves do not appear to have single-exponential decay characteristics that my code assumes.

.

##### 2 Comments

Star Strider
on 24 Mar 2023

Fitting four parameters with three data pairs is actually not possible, or only minimally possible if you estimate one less parameter (for example eliminating ‘b(4)’), otherwise, there are likely an infinity of curves that could fit those four parameters, since they would not be unique. Estimating three parameters with three data pairs is actually not appropriate. Notice that the curve intersects those three points precisely, however the sine curve itself does not approximate the actual data.

For example, these two illustrations —

y = [ 0.0060 0.0077 0.0058];

x = 0:1:2;

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,rn] = 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

per = 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,rn] = 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

per = 4*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,rn] = 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

All of these waveforms fit the data, with approximately the same residual norm (‘rn’) values.

.

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!