processTimeStretch(data,stAlgo)

function [stAlgo, status, data] = processTimeStretch(data,stAlgo)
%     Time-scale a signal to r times faster with or without phase vocoder
% 2000-12-05, 2002-02-13 [e-mail].  Uses pvsample, stft, istft
% $Header: /homes/dpwe/public_html/resources/matlab/RCS/pvoc.m,v 1.2 2002/02/13 16:14:54 dpwe Exp $
status = 0;

if stAlgo.pitchcorrection == 0
    data = resample(data,100,(stAlgo.stretchfactor*100));
else

n = 1024;

% With hann windowing on both input and output, 
% we need 25% window overlap for smooth reconstruction
hop = n/4;
% Effect of hanns at both ends is a cumulated cos^2 window (for
% r = 1 anyway); need to scale magnitudes by 2/3 for
% identity input/output
scf = 2/3;

% Calculate the basic STFT, magnitude scaled
X = scf * stft(data', n, n, hop);

% Calculate the new timebase samples
[rows, cols] = size(X);
t = 0:stAlgo.stretchfactor:(cols-2);
% Have to stay two cols off end because (a) counting from zero, and 
% (b) need col n AND col n+1 to interpolate

N = 2*(rows-1);

if hop == 0
  % default value
  hop = N/2;
end

% Empty output array
c = zeros(rows, length(t));

% Expected phase advance in each bin
dphi = zeros(1,N/2+1);
dphi(2:(1 + N/2)) = (2*pi*hop)./(N./(1:(N/2)));

% Phase accumulator
% Preset to phase of first frame for perfect reconstruction
% in case of 1:1 time scaling
ph = angle(X(:,1));

% Append a 'safety' column on to the end of b to avoid problems 
% taking *exactly* the last frame (i.e. 1*b(:,cols)+0*b(:,cols+1))
b = [X,zeros(rows,1)];

ocol = 1;
for tt = t
  % Grab the two columns of b
  bcols = b(:,floor(tt)+[1 2]);
  tf = tt - floor(tt);
  bmag = (1-tf)*abs(bcols(:,1)) + tf*(abs(bcols(:,2)));
  % calculate phase advance
  dp = angle(bcols(:,2)) - angle(bcols(:,1)) - dphi';
  % Reduce to -pi:pi range
  dp = dp - 2 * pi * round(dp/(2*pi));
  % Save the column
  c(:,ocol) = bmag .* exp(j*ph);
  % Cumulate phase, ready for next frame
  ph = ph + dphi' + dp;
  ocol = ocol+1;
end

data = istft(c, n, n, hop)';
end