function out=gfilter(s,dt,n)
% GFILTER filter with cosine taper 
% gfilter(s,dt)
% Megnin's filter for COSY
% t1=31.4;
% t2=37;
% t3=250;
% t1=400;

 t1=31.4;
 t2=37;
 t3=250;
 t4=400;


% pad up to n^2
not=length(s);

imax=15;

for i=1:imax,
	if not >= 2^i,
	   if not < 2^(i+1),
	   ii=i+1;
	   end
	   if not == 2^i,
	   ii=i;
	   end
	end
end

s(not:2^ii)=(1:2^ii-not+1)*0.;
s=taper(s,-200);
s=taper(s,200);

disp(sprintf(' Padded %g points up to %g ',not,2^ii));


fs=fft(s);

% build filter

nt=length(fs);
fny=1./(2*dt);
finc=fny/(nt/2);

w1=1./t4;
w2=1./t3;
w3=1./t2;
w4=1./t1;

for i=1:nt/2,

w=i*finc;

if w<w1, filter(i)=0.;end

if w>=w1,
   if w<=w2,
     eta=pi*(w-w1)/(w2-w1);
     filter(i)=.5*(1-cos(eta));
   end
end

if w >=w2,
   if w <=w3,
     filter(i)=1.;
   end
end

if w >= w3,
   if w <= w4,
   eta=pi*(w-w3)/(w4-w3);
   filter(i)=.5*(1+cos(eta));
   end
end

if w > w4, filter(i)=0.;end

filter(nt+1-i)=filter(i);

end


filter=filter+j*filter;

% filter 

fs=fs.*filter';
l=length(fs);

ns=2^n*l;
rfs=j*(1:ns);
rfs=rfs*0;
rfs(1:l/2)=fs(1:l/2);
rfs(ns-l/2+1:ns)=fs(l/2+1:l);

%if n==1,
%rfs(1:l/2)=fs(1:l/2);
%rfs(l/2+1:l/2+l)=(l/2+1:l/2+l)*0;
%rfs(l/2+l+1:2*l)=fs(l/2+1:l);
%end

%if n==2,
%rfs(1:l/2)=fs(1:l/2);
%rfs(l/2+1:l/2+3*l)=(l/2+1:l/2+3*l)*0;
%rfs(l/2+3*l+1:4*l)=fs(l/2+1:l);
%end

% backtransform

sfilt=ifft(rfs);
out=real(sfilt);
out=out(1:2^n*not);


% end


