function [X,opix] = ola(Y,H)
% X = ola(Y,H)
%   Overlap-add columns of Y, offset by H.
%   There's no windowing (i.e. scaling by a tapered window) in here. 
% 2010-11-14 Dan Ellis dpwe@ee.columbia.edu

[W,nw] = size(Y);

% How long should X be in the end?
% one complete window, plus another hop's worth of points for each
% additional window
lenx = (W + (nw-1)*H);

% The approach is to build an index matrix, where each row pulls 
% in successive, nonoverlapping windows from the original data.
% This matrix then has as many rows as needed to get all the
% overlapping windows.  We can then index the input with it, and 
% simply sum up.

% Upper bound on the the number of windows that overlap for each
% point = the number of rows we'll need in our index matrix
ovf = ceil(W/H);

% How many windows will there be in each row?
nwonf = ceil(nw/ovf);

% Unravel Y, then add a zero on the end.  Any indices we don't need
% in our index will point here, so they have no contribution to the
% final sum
Ywz = [Y(:)',0];
% (index for the final zero value)
leny = length(Ywz);

% Here's the indexing magic.  Create one row of the final index
% matrix that pulls out every ovf'th window from the original
% matrix in its unravelled form.  Pad any gaps with the zero index.
opix = reshape(repmat([[1:W],leny*ones(1,ovf*H - W)]', 1, nwonf) ...
               + repmat(ovf*W*[0:(nwonf-1)], ovf*H, 1), ...
               1, ovf*H*nwonf);

% need to pad out opix to cover hangover of final window
opix = [opix,leny*ones(1,lenx-length(opix))];

lopix = length(opix);
% Add in the following rows, which are offset by H points, and with
% indices that are shifted to pick up the interleaving windows
for i = 2:ovf
  opix = [opix;[leny*ones(1,(i-1)*H),W+opix(1,1:(lopix-(i-1)*H))]];
end

% Any points off the end get set to the zero value
opix(opix > leny) = leny;

% Now, the big payoff - index by the rows, including pulling in the
% zeros, in such a way that it's all ready to sum up to get OLA
X = sum(Ywz(opix))';