function [wl_l Hs] = Band_filter_V2(T,wl,f)

%% input
tolerance = 0.001;
fmax= 1;
%% Adjust data
% average water level
wlm = mean(wl);
% zero centered signal
wl = wl-wlm;
% detrend
wl = detrend(wl);
%% Fourier analyse
%sampling freq
fs = 1/(T(2) - T(1));
% nyquist freq
fn = fs/2;
% maximum freq
fn = min(fmax,fn);
n       = length(wl);
vy      = fft(wl,n);
df      = 1/(T(end)-T(1));

%vy = abs(vy);
fx = df*[0:1:round(length(T)/2) -1*floor(length(T)/2)+1:1:-1]';%'+0.5*df;

%% afjust frequentie to the fx spacing
indices = find(abs(fx-f)<df);
indices = indices(1);
f = fx(indices);

%% band filter
ind = ~rem(fx,f);
ind_t = ind';
for i=1:ceil(tolerance/df)
    indplus = [0 ind'];
    indplus = indplus(1:end-1);
    indmin = [ind' 0];
    indmin = indmin(2:end);
    % select surounding points
    ind_t = ind_t + indplus + indmin;
end
% make sure never get a value higher 1
ind_t(ind_t>1) =1;
% make a logical
ind_t = logical(ind_t);
% select frequencies
vy(~ind_t) = 0;
vy_l = vy;
% remove f > f nyquist
vy_l(fx>fn) = 0;

%% Compute water level
wl_l = real(ifft(vy_l,n));
wl_l = wl_l + wlm;
%% compute Hs low
%E         = 0.5 * vy_l.^2* T(end)/n;
m0        = trapz(fx,vy_l(:));
Hs        = 4.*sqrt(m0);

end
