//
//  timetrack.sci
//
//  ------------------------------------------------------------------------
//  copyright 2007  Alfred Steffens Jr.
//  -----------------------------------------------------------------------
//  -----------------------------------------------------------------------
//  Copying Permission:
//
//    This is free software; you can redistribute it and/or modify
//    it under the terms of the GNU General Public License as published by
//    the Free Software Foundation.
//
//    This software is distributed in the hope that it will be useful,
//    but WITHOUT ANY WARRANTY; without even the implied warranty of
//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//    GNU General Public License for more details.
//
//    You should have received a copy of the GNU General Public License
//    along with this software (see file called "COPYING"); if not, write to
//    the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
//    MA 02111-1307  USA
//  ------------------------------------------------------------------------
//
//  2nd step in building a frequency track for additive synthesis
//  (step 1 was ffttrack.sci)
//
//  1.  Assumes that ffttrack.sci was just run, so that a matrix called
//      dataF[] exists with the frequencies and amplitudes.
//
//  2.  Does a spline fit on each frequency track.
//
//  3.  Creates 3 matrices: freqmat, realmat, and imagmat in which each row
//      contains the data for a particular frequency harmonic as it changes
//      over time.
//
//  Set the values of "row1" and "rowN", which are the beginning and ending
//  rows in the data matrix, "dataF".
//
getf('fit_the_track.sci', 'c');

Fs = 44100;
PcMx = 0.1;             // fraction of max amp to keep
dt = 1 / Fs;
[Rows, Cols] = size(dataF);
row1 = 1;
rowN = Rows;
FRQ = 0;
//
//  the time increments
//
tm = dataF(row1:rowN,1);
Tdel = max(tm) - min(tm);
M = int(Fs * Tdel);              // number of samples to create
N = 1:M;                    // array of sample numbers
xd = linspace(min(tm), max(tm), M);   // time points
//
//
//
freqmat = [];
realmat = [];
imagmat = [];
//
//
//  frequencies 1
//
FRQ = FRQ+1;
col1 = (FRQ-1)*3 + 2;
col2 = (FRQ-1)*3 + 3;
col3 = (FRQ-1)*3 + 4;
f1 = dataF(row1:rowN,col1);
r1 = dataF(row1:rowN,col2);
i1 = dataF(row1:rowN,col3);
//
// spline fit of the data
//
ydf = interpln([tm' ; f1'], xd);        // 1st cut at curve
[yl, dl] = lsq_splin(xd, ydf, tm');
f1s = interp(xd, tm', yl, dl);         // frequencies
ydr = interpln([tm' ; r1'], xd);        // 1st cut at curve
[yl, dl] = lsq_splin(xd, ydr, tm');
r1s = interp(xd, tm', yl, dl);         // real amplitudes
ydi = interpln([tm' ; i1'], xd);        // 1st cut at curve
[yl, dl] = lsq_splin(xd, ydi, tm');
i1s = interp(xd, tm', yl, dl);         // imag amplitudes
//
//
//
freqmat = [freqmat ; f1s];
realmat = [realmat ; r1s];
imagmat = [imagmat ; i1s];
//
//
//  the rest of the frequencies
//
for  k = 3:10,

    FRQ = FRQ+1;
    col1 = (FRQ-1)*3 + 2;
    col2 = (FRQ-1)*3 + 3;
    col3 = (FRQ-1)*3 + 4;
    f2 = dataF(row1:rowN,col1);
    r2 = dataF(row1:rowN,col2);
    i2 = dataF(row1:rowN,col3);

    [f22s, r22s, i22s] = fit_the_track(tm, xd, f2, r2, i2);

    //  zero any negative frequencies in the freq row
    fmx = max(f22s) * PcMx;
    indfr = find(f22s < fmx);
    f22s(indfr) = zeros(indfr);

    freqmat = [freqmat ; f22s];


    r22s(indfr) = zeros(indfr);
    realmat = [realmat ; r22s];


    i22s(indfr) = zeros(indfr);
    imagmat = [imagmat ; i22s];
end


fprintfMat('freqmat.dat', freqmat, '%f');
fprintfMat('realmat.dat', realmat, '%f');
fprintfMat('imagmat.dat', imagmat, '%f');