function [xOut,tOut] = TVFDfilter(x,fs,beta,Tau,N)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% [xOut,tOut] = TVFDfilter(x,fs,beta,Tau,N)
%
% Computes a time-stretched/compressed and time-shifted data-set
% with an N-point Lagrange fraction delay filter as defined in
% "Post-processed time-delay interferometry for LISA" by Shaddock et. al.
%
% INPUTS:
%  x = time series input [amplitude]
%  fs = sampling frequency [Hz]
%  beta = fractional shift rate (in s/s)
%  Tau = absolute time-shifting value in seconds
%  N = interpolation length (number of points to use in
%     Fractional Delay Filter)
%
% OUTPUTS:
%  xOut = output amplitude vector
%  tOut = output time vector
%
%   Author: Shawn Mitryk
%   Created: 9-18-11
%   Updated: 3-22-12
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


	L = length(x);
	M = (N-1)/2;
	k = (-(N-1)/2):1:((N-1)/2);

	offset = floor(Tau*fs);
	Dfrac = -(Tau-(offset/fs))*fs;
	xOutShift(1:L-N+1) = 0;


for i=1:L-N+1
	DTaufrac = beta*i;
	td = (N-1)/2+Dfrac-DTaufrac;
	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;
	h = sinc(Dfrac-DTaufrac-k).*w;
	xOutShift(i) = sum(h.*(x(i:i+N-1)'));
end


xOut = xOutShift(N+M+1:L-N+1)';
tOut = (0:length(xOut)-1)'/fs;

