function rewrite_aligned_audio(inputname, outputname, offset, skew, maxdur)
% rewrite_audio(inputname, outputname, offset, skew, maxdur)
%    Read in a waveform; write out a new version such that the
%    output is the original *advanced* by <offset> (so -ve offset
%    introduces silence at the start of the output), and sped up by
%    factor <skew> (so skew > 1 leads to shorter output)
%     t_output = (t_input - offset) / skew
%    Optional maxdur will truncate the output at that duration.
% 2013-08-30 Dan Ellis dpwe@ee.columbia.edu
% adapted from deskew.m via skewview.m

if nargin < 3; offset = 0; end % then why are you calling me?
if nargin < 4; skew = 1; end
if nargin < 5; maxdur = 0; end

disp(['rewrite_aligned_audio(', inputname,' ',outputname,' ', ...
      num2str(offset), ' ', num2str(skew), ' ', num2str(maxdur), ')']);

% Find p/q s.t. p/q ~= a -- approx 1 part in 1/(a-1)
p0 = floor(1/abs(skew-1));
%disp(['p0=',num2str(p0)]);
if p0 < 2^15
  % exhaustive search for pair of integers closest to desired rate
  p = p0:(2^15);
  q = round(p./skew);
  er = (skew - p./q); 
  [ee,xx] = min(abs(er));
  P = p(xx);
  Q = q(xx);
else
  % too close to unity
  P = 1;
  Q = 1;
end
disp(['Rewrite: resampling ratio: ',sprintf('%d/%d=%.6f',P,Q,P/Q)]);

[pi,ni,ei] = fileparts(inputname);
[po,no,eo] = fileparts(outputname);
if strcmp(ei,'.mp3') && strcmp(eo,'.mp3')
  % can use streaming for both input and output

  % open input and output streams
  [YI,SR,CH] = mp3streamin_open(inputname);
  YO = mp3streamout_open(outputname, SR, CH);

  if offset < 0
    % prepad silence
    padsil = round(-offset*SR);
    % is silence before or after time scaling?
  else
    % -ve offset is time to trim from input file
    % discard some of the input stream?
    mp3streamin_read(YI, round(offset*SR));
    padsil = 0;
  end
  
  % algorithm from resample_in_parts
  targetblock = 1e7;  % say

  % How many blocks of P can we do at once?
  numP = round(targetblock/P);
%  numloops = ceil(length(X)/(numP*P));

  done = 0;
  Ydone = 0;
  %todo = length(X);
  alldone = 0;
  nextX = zeros(0,CH);
  block = 0;
  
  while ~alldone
    if done == 0; pad = 0; else pad = 1; end
    thistime = numP*P;
    %nextX = X(done - pad*P + [1:(thistime+pad*P)]);
    % Keep padding from last frame
    nextX = nextX(end-pad*2*P+1:end,:);
    % Read in from padsil if nonzero
    padsilnow = min(thistime+pad*P, padsil);
    nextX = [nextX; zeros(padsilnow, CH)];
    padsil = padsil - padsilnow;
    % Pull rest from file
    nextX = [nextX; ...
             mp3streamin_read(YI, thistime + (1-pad)*P - padsilnow)];
    if size(nextX,1)-pad*2*P < thistime
      % Reduce thistime if we couldn't get the whole amount from file
      thistime = min(size(nextX,1)-2*pad*P, thistime);
      % and note that we've hit EOF
      alldone = 1;
    end
    Ypart = [];
    for ch = 1:CH
      Ypart(:,ch) = resample(nextX(:,ch), Q, P);
    end
    % Discard prepad
    Ypart = Ypart((pad*Q+1):end, :);
    % Only keep the valid parts
    lYpart = min(size(Ypart,1), numP*Q);
    %Y(Ydone + [1:lYpart]) = Ypart(1:lYpart);
    mp3streamout_write(YO, Ypart(1:lYpart,:));
    Ydone = Ydone + lYpart;
    % forget postpad when winding forward
    thistime = min(thistime, numP*P);
    done = done + thistime;
    %todo = todo - thistime;
    block = block + 1;
  end
  
  % close streams
  mp3streamin_read(YI,-1);
  mp3streamout_write(YO,[]);
  
else
  % not mp3 to mp3 - resort to in-memory
  % re-read targ
  [dr,sr] = audioread(inputname);
  if offset < 0
    % prepad silence
    dr = [zeros(round(-offset*sr),size(dr,2));dr];
    % is silence before or after time scaling?
  else
    % -ve offset is time to trim from input file
    dr = dr(round(offset*sr)+1:end,:);
  end

  % Apply time scaling via resampling: 
  for i = 1:size(dr,2)
    dmr(:,i) = resample_in_parts(dr(:,i),Q,P);
    % resample_in_parts gives identical results but doesn't save time
  end
  dr = dmr;
  % Maybe limit duration
  if maxdur > 0
    dr = dr(1:min(round(maxdur*sr), size(dr,1)), :);
  end
  % Write it out
  audiowrite(dr,sr,outputname);
  CH = size(dr,2);
end

disp(['Wrote skewed audio to ',outputname,' from ',inputname, ...
      ' with ',sprintf('offset=%.3fs skew=%.6f maxdur=%.1fs', ...
                       offset,skew,maxdur)]);

disp('Equivalent sox command:')
offsetstretch = offset/skew;
if offset < 0
  cmd = ['sox ', inputname, ' ', outputname, ' speed ', sprintf('%.6f', skew), ...
         ' delay ', num2str(-offsetstretch)];
  if CH == 2
    % need to specify delay for both channels
    cmd = [cmd, ' ', num2str(-offsetstretch)];
  end
  if maxdur > 0
    cmd = [cmd, ' trim 0 ', num2str(maxdur)];
  end
else
  cmd = ['sox ', inputname, ' ', outputname, ' speed ', sprintf('%.6f', skew), ...
         ' trim ', num2str(offsetstretch)];
  if maxdur > 0
    cmd = [cmd, ' ', num2str(offsetstretch+maxdur)];
  end
end
disp(cmd);

