How to Remove the Periodic Oscillations from a Signal
Their are some problems in the code and same algorithmic issues:
Code:


You are estimating the wrong frequency: f2=1/((365*24*60)/2); should be f2=1/((365*24*60))/2;


Algorithm:
Basically you are using a fourier series (http://en.wikipedia.org/wiki/Fourier_series) with those terms corresponding to the known frequencies. But you should be aware of the limitations of that approach:


You have to estimate the dc part (frequency 0), because your signals' dc value is not zero!
In theory this approach works only for full periods of each frequency. That's why f1 is estimated quite good, but f2 is not. To work arround that you need either more data, because the more periods are in the data the smaller will be the estimation error. And you should use a windowing function, to get rid of some side effects by cutting the signal at the beginnig and the end.


The improved code could look like this:

clear
doWindow = true;
scaleFreq = 1;
%% Creating the data

g=2.5;
t=linspace(0,g,g*365*24*60);              %n-years of 1min data
o=rand(size(t));                         %random noise
%we create two sinusoids  
p1=2*sin(2*pi*t*scaleFreq);
p2=1.5*cos(1*pi*t*scaleFreq);
p=p1+p2;
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=scaleFreq*1/(365*24*60);       %frequency of 1 yr oscillations
f2=scaleFreq*1/((365*24*60))/2;   %frequency of 1/2 year oscillations

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

if(doWindow)
    window = hanning(numel(x)).';
else
    window = ones(size(x));
end
xWin = x.*window;%windowing
xMean = mean(xWin);%calculate dc
xWin = xWin-xMean;%subtract dc

buildFunkCos1 = cos((2*pi*a1*n)/N);
buildFunkCos2 = cos((2*pi*a2*n)/N);
buildFunkSin1 = sin((2*pi*a1*n)/N);
buildFunkSin2 = sin((2*pi*a2*n)/N);

A1=(1/N)*sum(xWin.*buildFunkCos1)
A2=(1/N)*sum(xWin.*buildFunkCos2)
B1=(1/N)*sum(xWin.*buildFunkSin1)
B2=(1/N)*sum(xWin.*buildFunkSin2)

if(doWindow)%do correction due to windowing
    A1 = 2*A1;
    A2 = 2*A2;
    B1 = 2*B1;
    B2 = 2*B2;
end

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

fi1=atan(B1/A1)%phase shift
if(A1&lt;0)
    fi1 = fi1 + pi;
end
fi2=atan(B2/A2)
if(A2&lt;0)
    fi2 = fi2 + pi;
end

%the reconstructed set of data for these two frequencies
v1=C1*cos(2*pi*t*scaleFreq-fi1);
v2=(C2*cos(1*pi*t*scaleFreq-fi2));
v=v1+v2+xMean;%sum estimates (freq 1,freq 2, dc)
r=x-v;    %I wasn't sure what was meant by removing, so i just subtracted

%% Visualization
figure(1)                                  %now lets plot the data
plot(t,[p1;p2;p]);
legend('p1','p2','p');
title('single sine signals and sum');

figure(2)                           %and, original signal and noise separately
hold off
plot(t,[v1;v2;v]);
legend('v1','v2','v');
title('estimated sine signals and sum');

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

figure(4)
hold off
plot(t,[buildFunkCos1;buildFunkCos2;buildFunkSin1;buildFunkSin2]),
legend('buildFunkCos1','buildFunkCos2','buildFunkSin1','buildFunkSin2');


Here you can switch the windowing and also test what happens if more periods are in the signal. Try, scaleFrequ=10 und doWindow=true and the result should be what you are looking for.