function [RawZeroedTS, RawMean, InterpMean] = ...
do_interp( interpFactor, rawTS, Y, PST, pst )

global GabLab;
if ( ~isfield(GabLab, 'DumpPlotDataToFile') ) 	GabLab.DumpPlotDataToFile = 0; end;

%	Text files that are written if GabLab.DumpPlotDataToFile = 1
%
FILE1	= 'interp_means.txt';
FILE2	= 'raw_means.txt';

%  Not sure what means to subtract here...
%
rawzTS 	= rawTS - mean(rawTS) + mean(Y);
HRFpts 	= 1:length(PST);
uniqPst	= unique(sort(pst));
uniqPst	= uniqPst(HRFpts);

%	set up time vectors for interpolation
%
t1 = uniqPst';                                        
t2 = linspace( t1(1), t1(end), interpFactor*length(t1))';    

for i = HRFpts		%% compute means of raw data
    lastOne	 = length(rawzTS) - length(PST) + i;
    goodVals	 = i:length(PST):lastOne;
    rawMean(i)   = mean( rawzTS( goodVals ) );
end

%	Sinc interpolation is from the examples in the Matlab Signal
%	Processing Toolbox manual
%
interpMean1 = interp( rawMean, interpFactor );
interpMean2 = sinc( t2(:,ones(size(t1))) - t1(:,ones(size(t2)))')*rawMean';

%	pst2 is a new pst vector rotated by one time point
%
pst2 = [pst(end) pst(1:end-1)];

% 	we want these to be long column vectors
%
RawZeroedTS 	= [ pst2' rawzTS ];
RawMean 	= [ uniqPst' rawMean' ];
InterpMean 	= [ t2 interpMean1' interpMean2 ];

%===========================================================================

%	Dump data to ASCII file if desired
%
if ( GabLab.DumpPlotDataToFile )

    %	Interpolated mean values
    %
    A = InterpMean( :, [1 2] );		% time & non-sinc interp mean
    %A = InterpMean( :, [1 3] );	% time & sinc interp mean
    fid = fopen( FILE1, 'w' );
    fprintf( fid, '%f\t%f\n', A' );
    fclose( fid );

    %	Interpolated mean values
    %
    fid = fopen( FILE2, 'w' );
    fprintf( fid, '%f\t%f\n', RawMean' );
    fclose( fid );
end
