它们是代码中的一些问题和相同的算法问题：代码：

• 您正在估计错误的频率：f2=1/((365*24*60)/2);应该是f2=1/((365*24*60))/2;

算法：基本上，您使用的是傅立叶级数（http://en.wikipedia.org/wiki/Fourier_series），这些术语对应于已知频率。但是您应该意识到这种方法的局限性：

• 您必须估计直流部分（频率 0），因为您的信号的直流值不为零！

• 从理论上讲，这种方法仅适用于每个频率的完整周期。这就是为什么 f1 估计很好，但f2不是。要解决您需要更多数据的问题，因为数据中的周期越多，估计误差就越小。并且您应该使用窗口功能，通过在开始和结束处切断信号来消除一些副作用。

改进后的代码可能如下所示：

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

在这里，您可以切换窗口并测试如果信号中有更多周期会发生什么。试一试，scaleFrequ=10结果doWindow=true应该是你想要的。

您可以使用傅里叶变换，将要消除的频率设置为零，然后进行傅里叶逆变换。

说你有功能$f(x)$和你 FT 获得功能$\hat{f}(w)$，然后定义$\hat{g}(w) = \hat{f}(w)$对全部$w \neq 365/2\pi$和$\hat{g}(365/2\pi)=0$. 现在旅游学院$\hat{g}$要得到$g(x)$这将是一个没有年波动的函数。

您可以使用作为递归滤波器实现的陷波滤波器。在计算上，这比两个傅里叶变换要有效得多。