% usage: [alpha, beta] = unwarp(frequency matrix, phase matrix, amplitude matrix)
% 
% This function calculates a matrix of values for M, alpha, and beta given the
% 3 indicated matrices, all of which should have the same dimensions.  M is a
% smoothing parameter for the phase unwrapping, and alpha and beta depend upon it.
% Alpha and beta themselves are used in inter.m as parameters of the cubic inter-
% polation function. 

function [alpha,beta] = unwarp(frek, fase, amp, winspace)

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

[m,q] = size(frek);

% The equations for M, alpha, and beta are given in the paper referenced at the
% beginning of this report.
for a=1:m-1
   for b=1:q
      if(amp(a,b)~=0)
	 M(a,b) = round( (1/(2*pi))*( (fase(a,b) + frek(a,b)*T - fase(a+1,b)) + (T/2)*(frek(a+1,b)-frek(a,b)) ) );
         alpha(a,b) = (3/T^2)*(fase(a+1,b) - fase(a,b) - frek(a,b)*T + 2*pi*M(a,b)) -1/T*(frek(a+1,b) - frek(a,b));
         beta(a,b) = -(2/T^3)*(fase(a+1,b) - fase(a,b) - frek(a,b)*T + 2*pi*M(a,b)) + (1/T^2)*(frek(a+1,b) - frek(a,b));
      end
   end
end