matlab 1 yr oscillations

doWindow = true;

scaleFreq = 1;

%% Creating the data

g=2;

t=linspace(0,g,g*365*24);

x=tiyingbian;

x=detrend(x);

%% The method

N=length(x);            %no of data

n=(1:N);

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

f1=scaleFreq*1/(365*24);       %frequency of 1 yr oscillations

f2=scaleFreq*1/((365*24))/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

x=x';

xWin = x.*window;%windowing

xMean = mean(xWin);�lculate 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<0)

    fi1 = fi1 + pi;

end

fi2=atan(B2/A2);

if(A2<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,x);

% 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');