fit a curve to data without using curve fitting toolbox

I have a set of data=[x1 x2] which looks periodical. I want to fit them into this Fourier transform equation:

x2 = A1 + A2.*sin(x1) + A3.*cos(x1) + A4.*sin(2*x1) + A5.*cos(2*x1)

by using least squares optimisation to know the optimum A = [A1;A2;A3;A4;A5]

I don't have curve fitting toolbox. please let me know how to do it without using toolbox.

[Andreas Goser]

I see you already have answers, but I wonder why you do not have the Curve Fitting Toolbox. Is the information about the university you work at current (your profile)?

Answer 1

[Mohammad Abouali]

Assuming your X1 and X2 are vector, i.e. size(x1)=size(x2)=[N 1] and N>=5 (since there are 5 coefficients, A1 to A5); then

C=[ones(numel(x1),1) sin(x1(:)) cos(x1(:)) sin(2*x1(:)) cos(2*x1(:))];

A=(C'*C)\(C'*x2(:));

or even:

A=C\x2(:);

[Mohammad Abouali]

By the way have a look at Harmonic ANalysis of Time Series (HANTS) code. It does pretty much the same thing except that it also looks for missing values and some smoothing features. So it does a bit more.

But if you say you want 2 frequencies pretty much that's what you get.

[John D'Errico]

The latter form using \ is of course always the best, although it is unlikely that the resulting equations would be ill-conditioned.

[Mohammad Abouali]

If you decided to increase the number of sine and cosine terms; it can be easily done. Let's say you want up to sin(n*x) and cos(n*x) then just construct C matrix as follow:

C=zeros(numel(x1),2*n+1); % n is the coefficient in sin and cos
C(:,1)=ones(numel(x1),1);
for i=1:n
 C(:,2*i:2*i+1)=[sin(i*x1(:)) cos(i*x1(:))];
end

the rest is the same to get A

A=C\x2(:);

in this Case A would have 2*n+1 coefficient though.

In the problem you defined above n=2;

[Anqi Li]

Thank you so much for your fast and simple answer. I am a beginner in programming. Could you kindly tell me what's the meaning by using "A=C\x2(:);" ?It works for me so far but I want to learn how to use "\". Thanks a lot!

[Mohammad Abouali]

so x2(:) makes sure that your data is a column vector data, i.e. size(x2(:))=[N 1]. I just wanted to be sure that it is not a row data but a vertical column vector. If it is already a column vector you can just replace it with X2 and drop (:).

A=C\x2 conceptually is pretty much solving C*A=x2; Your system of equations can be written as matrix C (the one with sine and cosine functions) multiplied by the unknown column vector, and x2 is your known variable vectors. once you do A=C\x2 pretty much you are solving the system of linear-equations. This operation is known as mldivide. If you want to know more about it go to mldivide help section, there you can find more information on how to use this.

[Anqi Li]

It's very useful! thank you very much!