function out=synthtrc(sr,file,synth,scal)
% out=synthtrc(sr,file,synth,scal)
%                        Last change: <2000-08-22 00:59:56 roebel>
% 
% additive synthesizer for sdif sinusoid tracks data
% Written by Axel Roebel
% This is Matlab!
% based on IRCAM's sdif mex-functions 
%
% Attention: This code is only a quick hack, bugs are included with 
%            high probablity!
%             
% Output:
%    out    Model Data as piecewise linear amplitude spline model
%           and piecewise quadratic phase trajectory (synth == 0)
%           or
%           synthesized sound (synth == 1 )
% Input:
%    sr    Sample rate of data
%    file  Filename to read data
%    synth (0 = return model data/1= return synthesiszed sound)
%    scal  time scaling factor (>1 is stretching)

% this function is restricted to read 1TRC data

if nargin < 4
  scal = 1;
end

if nargin < 3
  synth =1;
end

filetrc=[file,':::1TRC/1TRC'];
[data, header, framesig, matsig] =loadsdiffile(filetrc);
fprintf('Frames read %d \n',size(data,2))

% Modelstructure 
md=struct('a0',[],'na',[],'p0',[],'nw',[],'ind',[],'f',[]);

datind =[];
for i=1:length(data)
  if ~isempty(data{i})
    datind(end+1) = i;
  end
end

  
anzframes    = length(datind);
fram=data{datind(1)};  
partials=fram(1,1);
for i=1:anzframes
  fram=data{datind(i)};  
  for ii=1:size(fram,1)  
    if isempty(find(fram(ii,1)==partials))
      partials(end+1) = fram(ii,1);
    end
  end    
end

numparts = length(partials);

md(1:numparts)=md;

for i=1:numparts
  md(i).ind = partials(i);
end

for ifram=1:anzframes

  % Store Indices  
  frammat = data{datind(ifram)};
  for ii=1:size(frammat,1)
    
    ip=find(partials==frammat(ii,1));

    if isempty(md(ip).a0)
      fprintf('Init Partial %d \n',partials(ip))      
      md(ip).p0      = frammat(ii,4);
      md(ip).a0(1)   = frammat(ii,3);
      md(ip).f(1)    = frammat(ii,2);
      md(ip).na(1:2) = round(header(ifram,1)*sr);
      md(ip).nw(1:2) = round(header(ifram,1)*sr);      
    else
      md(ip).a0(end+1) = frammat(ii,3);
      md(ip).f(end+1)  = frammat(ii,2);
      md(ip).na(end+1) = round(header(ifram,1)*sr);      
      md(ip).nw(end+1) = round(header(ifram,1)*sr);      
    end
  end
end

% Time scaling
for i=1:numparts
  % position the nodes and make  double nodes at the end  
  md(i).nw=round(md(i).nw([1:end end])*scal);
  md(i).na=round(md(i).na([1:end end])*scal);
end  

% Transform frequency to phase
for i=1:numparts
  % skip all partials that have zero length
  % due to missing last frame bug in some versions of sdif reader     
  if length(unique(md(i).nw))>1 
    % B-Spline representation of frequency trajectory (radiant)
    sp=spmak(md(i).nw,md(i).f/sr*2*pi);
    
    % phase is integral of f with given initial phase
    sp=fnint(sp,md(i).p0(1));
    md(i).p0 = spbrk(sp,'c');
    md(i).nw = spbrk(sp,'t');  
  end
end

if synth == 1
  fprintf('Synthesizing \n');
  out=zeros(1,max([md.na]+1));

  for ip = 1:numparts
    % skip all partials that have zero length
    % due to missing last frame bug in some versions of sdif reader     
    if length(unique(md(ip).nw))>1 & max(md(ip).f) <sr/2 
      fprintf('Partial %d \n',ip)      
      
      % partial samples 
      ind = md(ip).na(1):md(ip).na(end);
      
      % B-spline for amplitude 
      sp=spmak(md(ip).na,md(ip).a0);
      % Eval amplitude 
      aval=spval(sp,ind);
      % B-spline for phase
      sp=spmak(md(ip).nw,md(ip).p0);
      % Eval phase
      pval=spval(sp,ind);
      
      % Add partial    
      out(ind+1)= out(ind+1)  + aval.*cos(pval);
    end
  end    
else
  ind = 1;  
  for ip=1:numparts  
    if max(md(ip).f )< sr/2
      out(ind) = md(ip);  
      ind = ind+1;
    end
  end
  
end

return

% NAME
%     loadsdiffile  -  Load a complete SDIF file
%
% USAGE 
%     [ data, header, frame, matrix ] = loadsdiffile (name)
%
%     where data{i} is the matrix at time header(i, 1) in stream header(i, 2) 
%     with signatures char(frame(i)), char(matrix(i)).
%
% SEE ALSO
%     loadsdifflat, loadsdif, sdifexist
% 
% AUTHOR
%     Diemo Schwarz (schwarz@ircam.fr), 31. January 2000
%
% CVS REVISION
%     $Id: loadsdiffile.m,v 1.3 2000/07/27 18:22:33 schwarz Exp $

% $Log: loadsdiffile.m,v $
% Revision 1.3  2000/07/27  18:22:33  schwarz
% man-style comments, better makedefs.
%
% Revision 1.2  2000/05/11  12:37:17  schwarz
% Init of output variables to avoid 'output argument not assigned' warning.
% Optional types file argument.
%
% Revision 1.1  2000/05/04  13:24:06  schwarz
% Matlab mex extension and support functions to load SDIF files.

function [ data, header, frame, matrix ] = loadsdiffile (name, types)

    data   = [];			% init to avoid 'non assigned' warning
    header = [];
    frame  = [];
    matrix = [];

    sdifexist (name);			% quit with error if file's not there

    if nargin == 1,			% open file with/without types
	NVTinfo = loadsdif (name);
    else
	NVTinfo = loadsdif (name, types);
    end
	
    n = 1;

    while (1)				% read frame by frame
	[ d, t, s, f, m ] = loadsdif;

% changed end condition to t=[]
% otherwise empty frames force termination of read
%       if isempty (d),  break;  end
	if isempty (t),  break;  end
    
	data{n}		= d;
	header(n, 1)	= t;
	header(n, 2)	= s;
	frame (n, 1:4)  = f;
	matrix(n, 1:4)  = m;

	n = n + 1;
    end
    
    loadsdif ('close');			% close files
return

% Usage:   [ data, time, stream, frame, matrix ] = loadsdiffile (name)
% where data{i} is the matrix at time(i) in stream(i) with signatures 
% frame{i}, matrix{i}.
%
%	frame{n}  = f;
%	matrix{n} = m;




%% Local Variables: 
%% time-stamp-start: "Last change:[ \t]+\\\\?[\"<]+"
%% End: 


