Matlab for Unwrapping Spectral Phase

When spectral phase is processed, it is often necessary to unwrap the phase to make it a ``continuous'' function. Below is a simple matlab function for this purpose. It is based on the assumption that phase jumps by more than $ \pi$ radians must have been ``wrapped''. That is, multiples of $ 2\pi$ are added or subtracted so that the phase changes by no more than $ \pm\pi$ from one spectral bin to the next. Bin 0 (corresponding to dc) is arbitrarily chosen as ``unwrapped'' and used as a starting point for the unwrapping algorithm.

Matlab listing: unwrap.m

function up = unwrap(p)
%UNWRAP unwrap phase

  N = length(p);
  up = zeros(size(p));
  pm1 = p(1);
  up(1) = pm1;
  po = 0;
  thr = pi - eps;
  pi2 = 2*pi;
  for i=2:N
    cp = p(i) + po;
    dp = cp-pm1;
    pm1 = cp;
    if dp>thr
      while dp>thr
        po = po - pi2
        dp = dp - pi2;
    if dp<-thr
      while dp<-thr
        po = po + pi2
        dp = dp + pi2;
    cp = p(i) + po;
    pm1 = cp;
    up(i) = cp;

Next Section:
Non-Parametric Frequency Warping
Previous Section:
Matlab for Computing Spectrograms