function [X,ph,pp] = synthphtrax(F, M, P, SR, WIN, HOP, DUR)
% X = synthphtrax(F, M, P, SR, WIN, HOP, DUR) Reconstruct from tracks incl phase.
% Each row of F, M and P contains a series of frequency, magnitude
% and phase samples for a particular track. These will be reconstructed
% and summed into the output sound X which will run at sample rate SR,
% although the columns in F and M are derived from an STFT with window
% WIN and advance HOP.
% If DUR is nonzero, X will be padded or truncated
% to correspond to just this much time.
% dpwe@icsi.berkeley.edu 1994aug20, 1996aug22
if(nargin<7)
DUR = 0;
end
rows = size(F,1);
cols = size(F,2);
cols = size(F,2);
opsamps = round(DUR*SR);
if(DUR == 0)
opsamps = WIN+((cols-1)*HOP);
end
X = zeros(1, opsamps);
% Time step
dt = HOP/SR;
% quadratic or cubic interpolation
iorder = 3;
for row = 1:rows
mm = M(row,:);
ff = F(row,:);
pp = P(row,:);
% Where mm = 0, ff is undefined. But interp will care, so set them.
% First, set all nan values of mm to zero
mm(find(isnan(mm))) = zeros(1, sum(isnan(mm)));
ff(find(isnan(ff))) = zeros(1, sum(isnan(ff)));
pp(find(isnan(pp))) = zeros(1, sum(isnan(pp)));
% For speed, chop off regions of initial and final zero magnitude -
% but want to include one zero from each end if they are there
nzv = find(ff);
firstcol = min(nzv);
lastcol = max(nzv);
zz = [max(1, firstcol-1):min(cols,lastcol+1)];
mm = mm(zz);
ff = ff(zz);
pp = pp(zz);
nzcols = length(zz);
% Find onsets - points where mm goes from zero (or NaN) to nzero
mz = (ff==0);
mask = mz & (0==[mz(2:nzcols),1]);
% Set frequencies at bins before onsets to match bin at onset
ff = ff.*(1-mask) + mask.*[ff(2:nzcols),0];
% Wind phase back to exactly match that freq
pp = pp.*(1-mask) + mask.*[pp(2:nzcols)-ff(2:nzcols)*dt*2*pi,0];
% Do offsets too
mask = mz & (0==[1,mz(1:(nzcols-1))]);
ff = ff.*(1-mask) + mask.*[0,ff(1:(nzcols-1))];
pp = pp.*(1-mask) + mask.*[0,pp(1:(nzcols-1))+ff(1:(nzcols-1))*dt*2*pi];
% Ok. Can interpolate now
% This is actually the slow part
ph = zeros(1,1+(nzcols-1)*HOP);
for col = 1:(nzcols-1)
% time index
tt = [0:(HOP-1)]/HOP*dt;
dp = pp(col+1) - pp(col);
ww = 2*pi*ff(col+[0 1]);
% How many cycles at given frq?
ncycu = ((ww(1) + ww(2)) * dt/2 - dp)/(2*pi);
ncyc = round( ncycu );
dp = dp + 2*pi*ncyc;
if iorder == 2
% Re-estimate final frq to match final phase
ww(2) = 2/dt * dp - ww(1);
% .. and recover quadratic phase coefficient
a = (ww(2) - ww(1))/(2*dt);
% Hence phase is...
ph((col-1)*HOP + [1:HOP]) = pp(col) + ww(1)*tt + a*(tt.^2);
else
% iorder = 3
a2 = 3/(dt*dt)*(dp - dt/3*(2*ww(1)+ww(2)));
a3 = 1/(3*dt*dt)*(ww(2)-ww(1)-2*a2*dt);
ph((col-1)*HOP + [1:HOP]) = pp(col) + ww(1)*tt + a2*(tt.^2) + a3*(tt.^3);
end
end
% Last value in phase function directly from input
ph((nzcols-1)*HOP + 1) = pp(nzcols);
% Simple linear interpolation of magnitudes
mmi = slinterp(mm, HOP);
% run the oscillator and apply the magnitude envelope
xx = mmi.*cos(ph);
% overlap-add it in to the correct place in the array
base = 1+WIN/2-HOP+HOP*(zz(1)-1);
sizex = prod(size(xx))-1;
xx = xx(1:sizex);
ww = (base-1)+[1:sizex];
X(ww) = X(ww) + xx;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Helper function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Y = slinterp(X,F)
% Y = slinterp(X,F) Simple linear-interpolate X by a factor F
% Y will have ((size(X)-1)*F)+1 points i.e. no extrapolation
% dpwe@icsi.berkeley.edu fast, narrow version for SWS
% Do it by rows
sx = prod(size(X));
% Ravel X to a row
X = X(1:sx);
X1 = [X(2:sx),0];
XX = zeros(F, sx);
for i=0:(F-1)
XX((i+1),:) = ((F-i)/F)*X + (i/F)*X1;
end
% Ravel columns of X for output, discard extrapolation at end
Y = XX(1:((sx-1)*F+1));