% usage: [Amplitude, Phase] = inter(amp_tracks, frequency_tracks, phase_tracks, alpha, 
% beta, winspace)
%
% This function interpolates the amplitude and phase of the signal. 
% It takes the matrices calculated by freqmat2.m and uses the along with the
% alpha and beta parameters to interpolate linearly for the amplitude and
% cubically for the phase. 

function [Amps, Phase]=inter(data, frek, faze, alpha, beta, winspace);  

[row, col] = size(data);

% an array of zeros to fill into the output arrays as appropriate
temp = zeros(winspace+1,1);

% T, the time between windows = # samples between windows * original sampling rate
T = winspace * 0.000125;

% n is a discrete time measurement starting at 0 and going up to the amount of time
% between consecutive frames
n = 0 : 0.000125 : T;
N = n';
for i = 1:col,
   % the row index into the new array  
   r_index = 1;
   for j = 1:row-1,
      % check to see if the two endpoints are 0, in which case insert more 0's into output arrays
      if data(j,i) == 0 & data(j+1,i) == 0
         Amps(r_index:r_index+winspace,i) = temp;
         Phase(r_index:r_index+winspace,i) = temp;   
    
      % check to see if neither endpoint is 0, in which case interpolate normally   
      % the interpolation equation is being implemented as described in the paper referenced
      % earlier in this report.
      elseif data(j,i)~=0 & data(j+1,i)~=0
         t = data(j,i) : (data(j+1,i)-data(j,i))/(winspace) : data(j+1,i);
         Amps(r_index : r_index+winspace, i) = t(:);
         fz = zeros(winspace+1,1)+faze(j,i);
         Phase(r_index : r_index+winspace, i) = (fz + N.*frek(j,i) + N.^2*alpha(j,i) + N.^3*beta(j,i));
   
      % check for the death of a track, in which case maintain the old frequency and reduce 
      % amplitude to 0 
      elseif data(j,i)~=0 & data(j+1,i)==0
	 t = data(j,i) : (data(j+1,i)-data(j,i))/(winspace) : data(j+1,i);
         Amps(r_index : r_index+winspace, i) = t(:);
	 fz = zeros(winspace+1,1)+faze(j,i);	
         Phase(r_index : r_index+winspace, i) = fz;
   
      % check for the birth of a track, in which case start with the next frequency (the first
      % measured for the track) and increase the amplitude up to the first measured value 
      elseif data(j,i)==0 & data(j+1,i)~=0
	 t = data(j,i) : (data(j+1,i)-data(j,i))/(winspace) : data(j+1,i);
         Amps(r_index : r_index+winspace, i) = t(:);
	 fz = zeros(winspace+1,1)+faze(j+1,i);	
      	 Phase(r_index : r_index+winspace, i) = fz;
    
      end % of if data() statements
   r_index = r_index + winspace;
   end % of for j loop
end % of for i loop