function y = fdfilter(x,D,window,N)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% output = fdfilter(input, delay, window, length) 
%%
%% Implements fractional delay interpolation filter.
%%
%% INPUTS:
%%   input    - data series to be interpolated
%%   delay    - # of samples to delay (-0.5 <= delay <= 0.5)
%%   window   - window type for filter
%%               supported windows:
%%                  'blackman' - Blackman Window
%%                  'blackman3' - Blackman Window^3
%%                  'lagrange' - Lagrange Filter
%%                  'none' - No window
%%   length    - length of filter kernel, must be odd.
%%
%% OUTPUTS:
%%   output    - interpolated data set.  The output begins
%%             at (N-1)/2 samples and goes to L-(N-1)/2 
%%             where N is the filter length and L is the
%%             number of samples.  For an input x(n) the
%%             output will then be x(n-(N-1)/2 + D), where
%%             D is the fractional delay.
%%
%% NOTES:
%%   For reference see "Post-processed time-delay interferometry
%%    for LISA" by D.A. Shaddock, et. al.
%%
%%   Author: Ira Thorpe
%%   Created: 8-29-06
%%   Modified: 9-21-06
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if mod(N,2) == 0
    disp('Error: Filter length must be odd');
    y = zeros(size(x));
end

k = (-(N-1)/2):1:((N-1)/2);

% select window function
switch lower(window)
    case {'rectangular','boxcar','none'}
        w = ones(1,N);
    case 'blackman'
        w = 0.42+0.5*cos((2*pi*k)/(N-1))+0.08*cos((4*pi*k)/(N-1));
    case 'blackman3'
        w = (0.42+0.5*cos((2*pi*k)/(N-1))+0.08*cos((4*pi*k)/(N-1))).^3;
    case 'lagrange'
        td = (N-1)/2+D;
        b1 = gamma(td+1)/(gamma(N+1)*gamma(td-N+1));
        b2 = gamma(N)./(gamma(k+(N-1)/2+1).*gamma(N-k-(N-1)/2));
        w = ((pi*N)/(sin(pi*td)))*b1*b2;
    otherwise
        disp('Error: unknown window');
        output = zeros(size(x));
end

% compute filter kernel
h = sinc(D-k).*w;
L = length(x);

%Use This section if the filter function is installed (much faster)
[junk, states] = filter(h,1,x(1:N));  % Compute Intial filter states
clear junk;
y = filter(h,1,x(N:L),states);        % filter subsequent data


% Use this section for systems w/o filter function (slow)
%for n = 1:L-N+1
%        y(n) = sum(h.*x(n:n+N-1));
%end