
function [XF] = lanczos_fltr(X,FILTER,HFc,LFc,N)

% LANCZOS_FLTR High, low, band pass filters based on Lanczos filter
%
% Y = LANCZOS_FLTR(X,FILTER_TYPE,HFc,LFc,N)
%
% Filter 1D signal X through a Lanczos filter.
%
% Inputs:
%	X: Input signal (1 dimension of length n)
%
%	FILTER_TYPE: Type of filtering you want:
%			1 : Low-pass  (cutoff frequency is LFc)
%			2 : High-pass (cutoff frequency is HFc)
%			3 : Band-pass (retain signal frequency within HFc<->LFc)
%
%	HFc: High frequency cut off (any k/n between 1/n and 1)
%	LFc: Low frequency cut off (any k/n between 1/n and 1)
%
%	N: (optional) Number of points in the filter. Note that the transition 
%     width of the transfer function is 2/(N-1)
%
% Outputs:
%	Y: Filtered signal
%
% Example:
%  Filter some time series x with a cutoff period of 10 years given annual
%  input data into x
%  Y=lanczos_fltr(x,2,1/10,1/10);
%
% Created by Nathan Steiger, UW Jan 2015
% Based on functions by some French guys...
%

if FILTER < 0 || FILTER > 3
   error('FILTER must be 1, 2 or 3 !');
end

X  = X(:);
if nargin < 5
   N  = length(X);
end
NJ = fix((N-1)/2);


switch FILTER
   case 1 % LOW PASS
      XF = lanczos(X,LFc,NJ);
      
   case 2 % HIGH PASS
      XF = X - lanczos(X,HFc,NJ);
      
   case 3 % BAND PASS
      XF = lanczos(X,LFc,NJ);
      XF = XF - lanczos(XF,HFc,NJ);
      
end %switch


end % prep function



% FILTER SUBFUNCTION

function [y] = lanczos(x,fn,nj)

% lanczos Filtre de lanczos
%
% [y] = lanczos(x,fn,nj)
%
% y  : output filtered data
% x  : input data
% fn : the cut off frequency (any k/n between 1/n , 0.5)
% nj : (number of points in the filter-1)/2 ; transition width of the
%      transfer function is 1/nj.
%
% The filter uses lanczos smoothing of a square box transfert function.

if nargin ~= 3
   error('Il faut 3 parametres en entree de la fonction lanczos')
end

nk = nj+1;

n=length(x);
%
% Pre-allocation memoire.
%
y=zeros(1,n);
a=zeros(1,nk);
b=zeros(1,nk);
ex=zeros(1,nk);
ey=zeros(1,nk);
pds=zeros(1,nk);
%
% Inverse vecteur en entree (on travaille sur vecteur colonne).
%
x=x';
%
% Calcul des poids du filtre.
%
pds(1)= 2*fn;
i=2:nk;
a(i)=pi*(i-1);
b(i)=2*a(i)*fn;
ex(i)=(sin(b(i))./a(i));
ey(i)=a(i)/nj;
pds(i)=ex(i).*sin(ey(i))./ey(i);

%
% Somme des poids.
%
ss=pds(1)+sum(2*pds(2:nk));

%
% Filtrage des nj premiers points.
%
for i=1:nj
   k2=2*i-1;
   j=1:k2;
   val=abs(i*ones(size(j))-j)+1;
   y(i)=sum(pds(val).*x(j));
   den=sum(pds(val));
   y(i)=y(i)./den;
end

%
% Filtrage du centre de la serie.
%
for i=nk:n-nj
   k1=i-nj;
   k2=i+nj;
   j=k1:k2;
   val=abs(i*ones(size(j))-j)+1;
   y(i)=sum(pds(val).*x(j));
   y(i)=y(i)/ss;
end

%
% Filtrage des nj derniers points de la serie.
%
for i=n-nj+1:n
   k1=2*i-n;
   j=k1:n;
   val=abs(i*ones(size(j))-j)+1;
   y(i)=sum(pds(val).*x(j));
   den=sum(pds(val));
   y(i)=y(i)./den;
end

%
% Inverse vecteur en sortie (pour obtenir vecteur ligne).
%
y=y';

end

