MATLAB: Sine fit in matlab vs closed-form expressions (linear least squares)

curve fittingleast squaresmldividesin

Given the arrays:
t = 0:0.25:2;
y = [3.0800    3.4900    3.1500    2.2300    1.2300    0.6900    0.8900    1.7300    2.7600];
The y array was created based on the expression 2+sin(3*t)+cos(3*t), added some noise and truncated to 2 decimal places.
We want to obtain the parameters of the model
y = A0 + A1*cos(w0*t) + B1*sin(w0*t)
This is a Fourier analysis, of course, but we want to study it in a least squares sense.
If we have N observations equispaced at intervals dt and a total record length T=(N-1)*dt, the coefficients can be determined as:
(Chapra, Applied Numerical Methods with Matlab, 4th ed)
So I code:
N = length(t)
A0 = sum(y)/N
A1 = 2*sum(y.*cos(3*t))/N
B1 = 2*sum(y.*sin(3*t))/N
And get the results
N =
9
A0 =
2.138888888888889
A1 =
1.337796008190744
B1 =
0.868485945765937
We can, however, use MATLAB mldivide operator, provided that, according to mldvide help:
If A is a rectangular m-by-n matrix with m ~= n, and B is a matrix with m rows, then A\B returns a least-squares solution to the system of equations Ax= B.*
I can code, then:
X = [ones(length(t),1) (cos(w0*t)).' (sin(w0*t)).']
a = X\(y.')
Which gives me the results:
a =
2.079400007449057
0.999055551110904
1.000625226465150
Those results are not the same as the ones obtained before. What am I missing here? Why those results are not the same?
I suspect I can not use the closed form expressions, but why exactly?

Best Answer

Hi Thales,
Before doing the least squares calculation it makes sense to try the less ambitious result of finding the right amplitudes without any added noise. Your time array has N = 9 points, and an array spacing of delt = 1/4 sec. The formulas only work if the fourier components are periodic with period T = N*delt. So the frequencies are multiples of w0 = 2*pi/T = 8*pi/9 ~~ 2.79. You are using w = 3 which isn't commensurate with w0, so you can't expect a fourier contribution at just one frequency, or that the fourier coefficient = 1.
N = 9;
x = (0:N-1)*(1/4);
w0 = 8*pi/9;
y = cos(x*w0);          % good results
(1/N)*sum(y)
(2/N)*sum(cos(x*w0).*y)
(2/N)*sum(cos(2*x*w0).*y)
(2/N)*sum(sin(x*w0).*y)
(2/N)*sum(sin(2*x*w0).*y)
% etc

ans =  -8.6351e-17
ans =  1.0000
ans =  -3.7007e-17
ans =   1.4803e-16
ans =  -2.4672e-17
y = cos(x*3);          % bad results
(1/N)*sum(y)
(2/N)*sum(cos(x*w0).*y)
(2/N)*sum(cos(2*x*w0).*y)
(2/N)*sum(sin(x*w0).*y)
(2/N)*sum(sin(2*x*w0).*y)
% etc
ans =  0.0695
ans =  1.0040
ans =  -0.0492
ans =  -0.2224
ans =   0.0210
Related Question