function q = unwrap(p,cutoff,dim)
%UNWRAP Unwrap phase angle.
%   UNWRAP(P) unwraps radian phases P by changing absolute 
%   jumps greater than pi to their 2*pi complement. It unwraps 
%   along the first non-singleton dimension of P. P can be a 
%   scalar, vector, matrix, or N-D array.
%
%   UNWRAP(P,TOL) uses a jump tolerance of TOL rather
%   than the default TOL = pi.
%
%   UNWRAP(P,[],DIM) unwraps along dimension DIM using the
%   default tolerance. UNWRAP(P,TOL,DIM) uses a jump tolerance
%   of TOL.
%
%   See also ANGLE, ABS.

%   Original: J.N. Little, 4-1-87.
%   Revised:  C R. Denham, 4-29-90. 
%             D.L. Chen, 3-22-95.
%   Copyright (c) 1984-98 by The MathWorks, Inc.
%   $Revision: 5.8 $  $Date: 1997/11/21 23:28:40 $

% Overview of the algorithm:
%    Reshape p to be a matrix of column vectors. Perform the 
%    unwrap calculation column-wise on this matrix. (Note that this is
%    equivalent to performing the calculation on dimension one.) 
%    Then reshape the output back.

% Treat row vector as a column vector.
rflag = 0;
if (ndims(p)==2) & (size(p,1)==1), 
   rflag = 1; 
   p = p.';
end;

% Initialize parameters.
nshifts = 0;
perm = 1:ndims(p);
if nargin == 1
   [p,nshifts] = shiftdim(p);
   cutoff = pi;     % Original UNWRAP used pi*170/180.
elseif nargin == 2
   [p,nshifts] = shiftdim(p);
else    % nargin == 3
   perm = [dim:max(ndims(p),dim) 1:dim-1];
   p = permute(p,perm);
   if isempty(cutoff), cutoff = pi; end
end;
   
% Reshape p to a matrix.
siz = size(p);
p = reshape(p, [siz(1) prod(siz(2:end))]);

% Unwrap phase angles.
[m, n] = size(p); 

pmin = min(p); pmin = pmin(ones(m, 1), :);   % To force REM to behave.
p = rem(p - pmin, 2 .* pi) + pmin;           % Phases modulo 2*pi.

if m > 1
   b = [p(1, :); diff(p)];                   % Differentiate phases.
else
   b = p;
end
c = -(b > cutoff); d = (b < -cutoff);        % Locations of jumps.
e = (c + d) .* 2 .* pi;                      % Array of 2*pi jumps.
f = cumsum(e);                               % Integrate to get corrections.

q = p + f;                                   % Phases + corrections.

% Reshape output
q = reshape(q,siz);
q = ipermute(q,perm);
q = shiftdim(q,-nshifts);
if rflag, q = q.'; end
