function [Fn,An,Pn] = FFT_for(fs,x)
f=50;
m=3;
N = length(x);
Fn=zeros(m,1);% 谐波频率
An=zeros(m,1);% 幅值
Pn=zeros(m,1);% 初相角
fsN=fs/N;
f0=50;
% 2-Hanning FFT
hw=hann(N,'periodic');
Xh=fft(x.*hw');
Xh=Xh(1:N/2)*2;
Xhabs=abs(Xh);
for i= 1 : m
[Amax,index]=TriFind(Xhabs,floor((i*f0-15)/fsN),ceil((i*f0+15)/fsN));
if(index==-1)
Fn(i) = 0;
An(i) = 0;
Pn(i) = 0;
else
if(Xhabs(index-1) > Xhabs(index+1))
a2 = Xhabs(index-1) / Xhabs(index);
r2 = (2-a2)/(1+a2);
k02 = index -1;
else
a2 = Xhabs(index) / Xhabs(index+1);
r2 = (2-a2)/(1+a2);
k02 = index;
end
Fn(i) = (k02+r2-1)*fs/N;
An(i) = 2*pi*r2*(1-r2*r2)*Xhabs(k02)/(N*sin(r2*pi));
Pn(i) = phase(Xh(k02))-pi*r2;
Pn(i) = mod(Pn(i),pi);
end
end
end