Asked by Thales
on 16 May 2019

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

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?

Answer by David Goodmanson
on 16 May 2019

Edited by David Goodmanson
on 17 May 2019

Accepted 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

David Goodmanson
on 16 May 2019

John D'Errico
on 17 May 2019

Another way of looking at this is to note that

omega = 3;

t = 0:0.25:2;

dot(cos(omega*t),sin(omega*t))

ans =

-0.09224

So over that interval, cos(omega*t) and sin(omega*t) are not orthogonal.

However, now change t.

t = linspace(0,2*pi,10);

dot(cos(omega*t),sin(omega*t))

ans =

7.2164e-16

they are now orthogonal. Those formulas are only valid if we have orthogonality.

David Goodmanson
on 17 May 2019

John makes a good point, but to demonstrate orthogonality in general, one has to be careful about the length of the array.

The extent of the time array here is 0 to 2*pi, so the fundamental frequency is omega0 = 1. Linspace(0,2*pi,10) does not represent a periodic array, because for the trig functions the last value is identical to the first value, meaning that one of the waveform points is listed twice. That works ok to show orthogonality if at least one of the basis functions is a sine function, since sin() = 0 at the end point, but not for two cos functions.

To create a periodic array, the last point in t needs to go away, leaving nine points.

t = linspace(0,2*pi,10);

omega0 = 1;

dot(cos(3*omega0*t),sin(3*omega0*t))

dot(cos(3*omega0*t),cos( omega0*t))

% ans = 7.0850e-16

% ans = 1.0000 % not 0

t(end) = []; % periodic array

dot(cos(3*omega0*t),sin(3*omega0*t))

dot(cos(3*omega0*t),cos( omega0*t))

% ans = 1.4433e-15

% ans = -1.7764e-15

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.