function [ ynew ] = sqbandpass( x,y,long,short )
%returns bandpass filtered version of data y
%   retains portion of y that gives rise to the portion of the signal
%   within long and short, which define the window of periods or
%   frequencies  this uses a simple rectangular frequency window
%   
f1=1/long;
f2=1/short;
n=length(x);
fc=fft(detrend(y));
fNyq=.5/abs(x(2)-x(1));

fre=linspace(0,2*fNyq,n)';

[m,k1]=min(abs(fre-f1)); % find k1 corresponding to f1
[m,k2]=min(abs(fre-f2)); % find k2 corresponding to f2
k3=n-k2+2; % conjugate points
k4=n-k1+2;
% zero all but the desired band
fc(1:k1)=zeros(k1,1);
fc(k2:k3)=zeros(k3-k2+1,1);
fc(k4:n)=zeros(k1-1,1);
%take inverse FFT; remove the imaginary part
ynew=real(ifft(fc));

end

