function F=bfmethod(f,L,h,eta,x)
%function F=bfmethod(f,L,h,eta,x)
N = max(size(h));
Hx = x(N) + [0:h(N):L-h(N)]; %coset M_N = x(N) + H_N
fH = f(1+Hx);%sampled values on M_N = x(N) + H_N;
SMf = SM(fH,L,h(N),x(N));
if N==1
  F = SMf;
else
  tmp = 1-exp(2*pi*i*([0:L-1]-x(N))*eta(N-1));
  tmp1=tmp; tmp1(find(abs(tmp < 1.e-14)))=1;%avoid zero divisions
  f1 = (f-SMf)./tmp1;
  fN1 = bfmethod(f1,L,h(1:N-1),eta(1:N-2),x(1:N-1));
  F = fN1.*tmp + SMf;
end
%----------------------------------------