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