%	open the file
clear all

dname = 'C:\Documents and Settings\Bruce Land\My Documents\WebSites\ee476site\Speech\Haskins\';
filename = 'atmelmega32parameters3.swi';
[fid, msg] = fopen([dname,filename]);
if fid == -1							
	error([msg, ' (', filename, ')']);
end;

%read SWI data
nOscs = fscanf(fid, '%f', 1);

theend = 0;
thistime = 1;

while (~theend)
    %need the try/catch because there may be bogus lines at the end of the
    %file
    try
        %the time stamp
        timeslice(thistime) = fscanf(fid, '%f', 1);
    catch
        break
    end
    % the freq and amplitude for each sinewave component
    for thisosc=1:nOscs
	 freq(thisosc,thistime) = fscanf(fid, '%f', 1);
	 amp(thisosc,thistime) = fscanf(fid, '%f', 1);
    end
    thistime = thistime + 1 ;
    theend = feof(fid);
end
fclose(fid);

%at this point we have three arrays: timeslice, freq, amp
%now interpolate and produce the sine waves

endTime = timeslice(end-3) * 0.001 ; %sceonds
%at 8 kHz get the sample times
Fs = 8000 ;
dt = 1/Fs ;
stepsize = timeslice(2)*8; % number of samples in timeslice(2) msec
t = 0:dt:endTime ;
y = zeros(size(t));

sliceindex = 1;
timeindex = 1;
slicetimeindex = 0;

F(1:nOscs) = freq(:,1);
dF = (freq(:,2)-freq(:,1))/stepsize;
phase(1:nOscs) = 0;
A(1:nOscs) = amp(:,1);
dA = (amp(:,2)-amp(:,1))/stepsize;

for time=t
    for i=1:nOscs
        %compute output signal
        y(timeindex) = y(timeindex) + ...
            A(i)*sin(phase(i));
        %compute new instantaneous frequency
        F(i) = F(i) + dF(i);
        A(i) = A(i) + dA(i);
        %and integrate the new phase
        phase(i) = phase(i) + 2*pi*F(i)*dt ;
    end
    y(timeindex) = y(timeindex)/nOscs ; %normalize
    
    timeindex = timeindex + 1; %running time
    slicetimeindex = slicetimeindex + 1; %time within a time slice
    
    if (slicetimeindex == stepsize)
        sliceindex = sliceindex + 1
        slicetimeindex = 0;
        F = freq(:,sliceindex);
        Fnext = freq(:,sliceindex+1);
        %check for fade out
         if Fnext==0
             Fnext = F;
         end
         A = amp(:,sliceindex);
         Anext = amp(:,sliceindex+1);
         dF = (freq(:,sliceindex+1)-freq(:,sliceindex))/stepsize;
         dA = (amp(:,sliceindex+1)-amp(:,sliceindex))/stepsize;
    end
end
sound(y,Fs)



