function [Z,E] = new_stxcorr(R,T,W,H,N)
% [Z,E] = new_stxcorr(R,T,W,H,N)  short-time cross-correlation
%    R (ref) is segmented into W-point stretches, hopped by H
%    samples. Each is slid over T over the range [-N:N] relative
%    points (or [N(1):N(2)] for two-element N), and the
%    cross-correlation is returned as columns of Z.
%    Result is timing of T relative to R, so that if T is equal to
%    some padding samples followed by R (i.e., R delayed) then the
%    cross correlation has peak values beyond its midpoint (in each
%    column of Z).
%    E returns theoretical maximum (harmonic mean energy) for each
%    col of Z.
% 2004-04-28 dpwe@ee.columbia.edu

% This version takes W-point chunks of R and slides them over T for
% the corresponding lags (-N:N) only.  Thus, the chunks of T are
% sized W+2N points, and the convolution is padded to W + W+2N =
% 2(W+N) points

if nargin < 2;   T = []; end
if nargin < 3;   W = 2000;  end
if nargin < 4;   H = round(W/2); end
if nargin < 5;   N = W; end

if length(N) == 1
  N = [-N,N];
end

%if N > W; N = W; end

npts = N(2)-N(1)+1;
LR = length(R);
nfr = 1+floor((LR-W)/H);

Z = zeros(npts, nfr);
E = zeros(1, nfr);

% Default to xcorr
if isempty(T)
  T = R;
end

% Prepad T with N zeros, so xcorr is centered
prepad = max(0, -N(1));
postpad = max(0, W+N(2)+1+H*(nfr-1)-length(T));
PT = [zeros(prepad, 1); T; zeros(postpad, 1)];
N = N + prepad;

% break it into blocks
maxarraysize = 10000000;
frblk = floor(maxarraysize/(2*W+npts));

nfrblk = ceil(nfr/frblk);

basefr = 0;
for blk = 1:nfrblk
  basefr = (blk-1)*frblk;
  nfrthisblk = min(frblk, nfr - basefr);
  
  RR = fft(frame(R(H*basefr + [1:(W+H*(nfrthisblk-1))]), W, H), 2*W+npts);
  TT = fft(frame(PT(H*basefr + [1:(W+npts+H*(nfrthisblk-1))]), W+npts, H), 2*W+npts);
  
  xx = ifft(TT.*conj(RR));
  % if the last block only has one frame, it returns result as a row
  if nfrthisblk == 1; xx = xx'; end

  xdest = [xx(1+[N(1):N(2)],:)];
  Z(:,basefr + [1:nfrthisblk]) = xdest; 
  E(basefr + [1:nfrthisblk]) = sqrt(sum(abs(RR).^2).*sum(abs(TT).^2))/(2*W+npts);
end