function [thetaCurrent,n,N]=FSFSingle(u,y)
%using frequency sampling filter to estimate single frequency
%based on a pair of periodic signals 
% If it is relay testing data, the output is from the negative error signal 
u=u-mean(u);
y=y-mean(y);
n=9;
L=length(u);
Ujw=fft(u);
[n1,m1]=max(abs(Ujw));
 P1=m1;
 N=round(L/(P1-1));
%Generate FSF filter output using Fourier transform
alpha=exp(j*2*pi/N); %basic element
beta=1;  %basic row
for i=1:N-1
beta=[beta alpha^i];
end
zeta=beta.^0;
for i=1:(n-1)/2
zeta=[zeta; beta.^i]; %coefficient matrix
end
%use the first N data points to estimate the initial recursive parameters
U_input=zeros(N,1);

for kk=1:N
Pt=zeta*U_input/N;
Pn=Pt(1);
 for i=2:(n+1)/2
Pn=[Pn;Pt(i);conj(Pt(i))];
end
Ps(:,kk)=Pn;
U_input=[u(kk);U_input(1:N-1,1)];
end
%least squares solution to obtain the initial estimate
M=inv(Ps(:,1:N)*Ps(:,1:N)');
thetaPast=M*Ps(:,1:N)*y(1:N);
P_ls=M;
lam=1;  %forgetting factor of one for time invariant system
%recursive least squares to estimate the parameters
for k=N+1:length(y)-1;
%updating FSF filters in real time    
U_input=[u(k);U_input(1:N-1,1)];
Pt=zeta*U_input/N;
Pn=Pt(1);
 for i=2:(n+1)/2
Pn=[Pn;Pt(i);conj(Pt(i))];
 end
%recursive least squares
  eCu=(y(k)-Pn'*thetaPast);
  thetaCurrent=thetaPast+P_ls*Pn*eCu/(lam+Pn'*P_ls*Pn);
  P_ls=(1/lam)*(P_ls-P_ls*Pn*Pn'*P_ls/(lam+Pn'*P_ls*Pn));
 thetaPast=thetaCurrent;
end

