function Y=BPF(X,fl,fh,fs,m)
  wh=2*pi*fh/fs;
  wl=2*pi*fl/fs;
  n=(-(m-1)/2:(m-1)/2);
  nZeros = find(n == 0);
  n(nZeros) = eps;
  Hh=sin(wh*n)./(pi*n);
  Hl=sin(wl*n)./(pi*n);
  H=Hh-Hl;
  Y=conv(X, H, 'same');
end

