function [yout,subst] = reconstruct(source, target, s_fs, t_fs, cheat, noise)

% Render the target sound in chunks of the source sound.  Currently
% takes chunks in the time domain.  Cheat is a number between 0 and
% 1 specifying how much of the original song to mix into the
% reconstruction.

if(s_fs ~= t_fs)
  error('Sampling rates need to be the same for now');
end
if(nargin < 5) cheat = 0; end
if(nargin < 6) noise = .1; end

frame = 512;

s_ons = onsets(source, s_fs, frame); zxi(5); pause(.1);
t_ons = onsets(target, t_fs, frame); zxi(5); pause(.1);
s_feat = calcFeatures(source, s_ons, s_fs);
t_feat = calcFeatures(target, t_ons, t_fs);

D = dist(s_feat, t_feat');
D = D + noise*rand(size(D))*diag(mean(D));
n = min(size(D));
D(sub2ind(size(D), [1:n], [1:n])) = inf;  % Prevent self-substitution
imagesc(D), colorbar
[dummy, subst] = min(D);

% Length of all of the source chunks
s_len = diff(s_ons')';

yout = 0 * target;
% $$$ yout = cheat * target;
for i=1:length(subst)
  yout(t_ons(i,1):t_ons(i,1)+s_len(subst(i))) = ...
      yout(t_ons(i,1):t_ons(i,1)+s_len(subst(i))) + ...
      window(source(s_ons(subst(i),1):s_ons(subst(i),2)));
end

if(cheat > 0)
  yout = [yout(:) target(:)];
end


function y = window(x)

%  Do some simple linear windowing to avoid clicks
len = 50;

if(length(x) < 2*len+2)
  w = triang(length(x));
else
  w = [[0:len]/len ones(1,length(x)-2*len-2) [len:-1:0]/len]';
end

y = x.*w;