function ry=FouFilter(y,samplingtime,centerfrequency,frequencywidth,shape,mode)
% Fourier filter function for time-series signal vector y; 'samplingtime'
% is the total duration of sampled signal in sec, millisec, or microsec;
% 'centerfrequency' and 'frequencywidth' are the center frequency and width 
% of the filter in Hz, KHz, or MHz, respectively; 'Shape' determines the 
% sharpness of the cut-off. If shape = 1, the filter is Gaussian; as  
% shape increases the filter shape becomes more and more rectangular. 
% Set mode = 0 for band-pass filter, mode = 1 for band-reject (notch) filter.  
% FouFilter returns the filtered signal.
% Example: plot(FouFilter(tan(1:1000),1,19,2,2,0))
%  T. C. O'Haver (toh@umd.edu),  version 1.5, May, 2007

center=centerfrequency*samplingtime; %  center harmonic (fourier component)
width=frequencywidth*samplingtime; %  width of filter (in harmonics)

fy=fft(y); % Fourier transform of signal
lft1=[1:(length(fy)/2)];
lft2=[(length(fy)/2+1):length(fy)];
% Compute filter shape
ffilter1=ngaussian(lft1,center+1,width,shape);
ffilter2=ngaussian(lft2,length(fy)-center+1,width,shape);
ffilter=[ffilter1,ffilter2];
if mode==1, 
    ffilter=1-ffilter; 
end
if length(fy)>length(ffilter), ffilter=[ffilter ffilter(1)];,end
ffy=fy.*ffilter;  % Multiply filter by Fourier transform of signal
ry=real(ifft(ffy)); % Recover filter signal from Fourier transform

function g = ngaussian(x,pos,wid,n)
%  ngaussian(x,pos,wid) = peak centered on x=pos, half-width=wid
%  x may be scalar, vector, or matrix, pos and wid both scalar
%  Shape is Gaussian when n=1, becomes more rectangular as n increases.
% Example: ngaussian([1 2 3],1,2,1) gives result [1.0000    0.5000    0.0625]
g = exp(-((x-pos)./(0.6006.*wid)) .^(2*round(n)));