How to remove the periodic oscillations from a signal?

The task that I have is to remove the annual and semiannual oscillation from a set of temperature measurements, taken over several years, by means of least squares method.

I found the method described here and since I have to wait for my colleague to first make some analysis on the data and pass them on to me, I decided to make a test run in Matlab, using a signal with two simple sine and cosine functions and random noise.

%%Creating the data

g=2.5;
t=linspace(0,g,g*356*24*60);              %n-years of 1min data
o=rand(size(t));                         %random noise
p=2*sin(2*pi*t)+1.5*cos(pi*t);           %we create two sinusoids  
x=p+o;                                     %and add noise to them

%%The method

N=length(x);            %no of data
n=(1:N);
T=g*365*24*60;          %time that the data covers
f1=1/(365*24*60);       %frequency of 1 yr oscillations
f2=1/((365*24*60)/2);   %frequency of 1/2 year oscillations

a1=f1*T;                
a2=f2*T;

A1=(2/N)*sum(x.*cos((2*pi*a1*n)/N));
A2=(2/N)*sum(x.*cos((2*pi*a2*n)/N));
B1=(2/N)*sum(x.*sin((2*pi*a1*n)/N));
B2=(2/N)*sum(x.*sin((2*pi*a2*n)/N));

C1=sqrt(A1^2+B1^2);                     %amplitude 
C2=sqrt(A2^2+B2^2);

fi1=atan(B1/A1);                        %phase shift
fi2=atan(B2/A2);

v=(C1*cos(2*pi*t-fi1))+(C2*cos(2*pi*t-fi2)); %the reconstructed set of data
                                             for these two frequencies

r=x-v;    %I wasn't sure what was meant by removing, so i just subtracted

%%Visualization

figure                                  %now lets plot the data
plot(t,x)
xlabel('t[min]')
ylabel('Arb. units :)')
title('Data')

figure                           %and, original signal and noise separately
plot(t,p,'b')
hold on
plot(t,o,'-r')
xlabel('t[min]')
ylabel('Arb. units :)')
title('Data separated')

figure                    %and once again data (blue), reconstructed signal 
plot(t,x,'b-')            %(green), and the final result afer subtraction
hold on                   %(red)
plot(t,v,'g--')
plot(t,r,'r:')
xlabel('t[min]')
ylabel('w/e')

The problem now is, that the reconstructed signal (v) does not resemble the original signal so much. Even when I remove the noise and work only with the pure sinusoid signal, what I get is not really a good representation of it. So, what I'm most interested is:

-Is this at all the correct use of this method? Am I using the correct equations?

-Is there an error somewhere in the code? By error I mostly consider the wrong definition, or calculation of some of the variables.

-In this "ideal case" (sinusoid without noise) should the method be able to reproduce the original curve (almost) exactly?

I hope I've laid out my problem clearly and understandably, if anyone has additional questions, please ask. I know this is rather basic stuff, but my knowledge and experience in this area is very limited (up to non-exsistant), so please bear with me.

Thank you in advance for any help.

[Tsehaye Negash]

What is g represent?

[Image Analyst]

@Tsehaye Negash see how annoying and confusing it is when people don't use descriptive variable names? By looking at

t=linspace(0,g,g*356*24*60);              %n-years of 1min data

we can deduce that "g" is what he also calls "n" in the comment, which is the number of years. So in this case he's setting it up to look at 2.5 years.

If he had simply called it numYears or numberOfYears instead of g, you would know instantly what it was. My two biggest peeves in coding are short nondescriptive variable names, and too few (or no) comments.

[Tsehaye Negash]

brother, 356 is not one year,ok

[Image Analyst]

@Tsehaye Negash you are correct. They must have made a typo. They did get it correct later though:

T=g*365*24*60;          %time that the data covers

Answer 1

[Image Analyst]

Try the Curve Fitting Toolbox. I don't have it so I can't help you much more. Otherwise, try taking the FFT and zeroing out low frequencies until you get a signal without much seasonal oscillations.

You could also try John D'Errico's curve fitting app: SLM