function y = blockresample(x,p,q,varargin)
% BLOCKRESAMPLE - block method of re-sampling
%   y = blockresample( x, p, q, N, beta )
%
%   (N, beta optional)
%
%   For x, p, q, N, beta see RESAMPLE, *except* that blockresample only
%   operates on x a vector
%
%   BLOCKRESAMPLE reduces memory usage for re-sampling by doing it on fixed
%   length blocks and then stiching the time-series back together.  It uses
%   'resample'. The goal is to provide mathematically the same result as a
%   straight call to re-sample 
%
% $Id: blockresample.m 2280 2005-07-04 00:06:38Z lsf $

% Goal:
%    resample input in blocks such that input blocksize + output blocksize
%    is kept close to a memory goal

% Parse input
if(nargin < 4)
    N = 10;
else
    N = varargin{1};
end
if(nargin < 5)
    beta = 5;
else
    beta = varargin{2};
end

% only block resample if x is a vector
if (numel(x) ~= max(size(x)))
    y = resample(x,p,q,N,beta);
    return
end

% reduce p/q to relatively prime fraction
[p,q] = rat(p/q);

% Mathematically the resample operation involves three steps
% 1. Zero insertion: x is expanded to length p*length(x) by zero insertion
% 2. low-pass filtering
% 3. decimation by keeping only every q sample
% To make the resample zero-delay the filtering is run forwards and
% backwards
% Correspondingly, to block resample using resample() we must
% A. Resample x in blocks that are a multiple of q in length
% B. Pad each block, ahead and behind, with data corresponding to the
%    length of the filter in samples of x (to set the filter state
%    correctly)
% C. Extract from each resampled, padded block, just the appropriate
%    resampled data (i.e., toss the parts that are associated with the
%    filter ringing up or ringing down)

% nFilt: length of fir filter. This is the # x samples that need to be
%        reserved to "restart" resample
nFilt = 2*q*N;
% nToss: how much of y to "throw-away"
nToss = nFilt*p/q;

% Memory goal: 
nMemMB = 2^6;
nMemB = nMemMB*2^20;

% x must be resampled in blocks that are multiples of q in length
% Arrange so that we don't use more than nMemMB MB data for input + output
%   x memory per min block: q*8 (for real);
%   y memory per min block: p*8 (for real);
%   total memory per min block: (p+q)*8 (for real)
%   block size (# x samples) (based on memory goal):
%       q*round(nMemB/((p+q)*8)); 
nXSamp = q*round(nMemB/((p+q)*8));
nYSamp = nXSamp*p/q;

% if nXSamp == 0 then we simply can't make a blocks small enough to meet
% memory goal 
if (0 == nXSamp)
  msgId = 'blockresample:qTooLarge';
  warning(msgId,'%s: q = %d',msgId,q);
%   NOTE we also use 'lmresample'  a low-memory resample which avoid a
%   memory hog at the end
  y = lmresample(x,p,q,N,beta);
  return
end

% get indices of block starts
% nFullBlocks: number of full blocks
% nSampRem: number of remainder samples
nFullBlocks = floor(length(x)/nXSamp);
nSampRem = length(x)-nFullBlocks*nXSamp;

% if nFullBlocks == 0 then we meet memory goal without blocking
if (0 == nFullBlocks)
%   NOTE we also use 'lmresample'  a low-memory resample which avoid a
%   memory hog at the end
    y = lmresample(x,p,q,N,beta);
    return
end

% if nSampRem < nFilt then we must increase nSampRem at expense of
% nFullBlocks
if (nSampRem < nFilt)
    nFullBlocks = nFullBlocks - 1;
    nSampRem = nSampRem + nXSamp;
end

% Resample x in blocks
%
% Pre-allocate result
% process full blocks
%   pre-allocate state (for first full block)
%   foreach full block
%     set-up block
%     resample block
%     extract result
%     place result
%     set-up state
%   ENDFOR
%   set-up remainder
%   resample remainder
%   extract result
%   place result

% shape y like x
[r,c] = size(x);
ny = ceil(length(x)*p/q);
if (r == 1)
    y = zeros(ceil([1,ny]));
else
    y = zeros(ceil([ny,1]));
end
preState = zeros([1,nFilt]);
xSub = zeros([1,nXSamp+2*nFilt]);
for iBlock = 1:nFullBlocks
    ndxX = (iBlock-1)*nXSamp;
    xSub(1:nFilt) = preState;
    xSub(nFilt+(1:(end-nFilt))) = x(ndxX + (1:(nXSamp+nFilt)));

    % Use lmresample - a low-memory resample that avoids memory
    % hogging steps at the end of resample

    ySub = lmresample(xSub,p,q,N,beta);

    ndxY = (iBlock-1)*nYSamp;
    y(ndxY + (1:nYSamp)) = ySub(nToss + (1:(end-2*nToss)));

    preState = x((1:nFilt) + (iBlock*nXSamp-nFilt));
end
ndxX = nFullBlocks*nXSamp;
% xSub = [preState(:)', x((ndxX+1):end)'];
xSub = zeros([1,nFilt+nSampRem]);
xSub(1:nFilt) = preState;
xSub(nFilt+(1:end-nFilt)) = x((ndxX+1):end);
% Use lmresample - a low-memory resample that avoids memory
% hogging steps at the end of resample
ySub = lmresample(xSub,p,q,N,beta);
ndxY = nFullBlocks*nYSamp;
y((end+1-length(ySub)+nToss):end) = ySub(nToss + (1:(end-nToss)));

return