[SciPy-user] phase wrapping algorithm

[Bryan Cole]

I do this quite a lot in the context of FFTs.
try:

def unwrap(thetaArray):
    """takes an array of theta values in degrees
    returns the data in unwrapped form"""
    diff = zeros_like(thetaArray)
    diff[1:] = thetaArray[1:] - thetaArray[:-1]
    upSteps = diff > 180
    downSteps = diff < -180
    shift = cumsum(upSteps) - cumsum(downSteps)
    return thetaArray - 360*shift

Obviously, a Weave version will be a lot faster, but this is at least 
fully vectorised, so it's not bad.

BC

On Tue, 2006-10-10 at 12:32 -0500, Ryan Krauss wrote:
> I am rethinking an algorithm for unwrapping phase data from arctan2
> for a vector of complex numbers (coming from lti models of feedback
> control systems).  I had written something a while ago using weave,
> but my old code causes seg faults with a recent scipy update.  So, I
> am rethinking the algorithm.
> 
> The problem occurs when either because of noise or just the definition
> of the inverse tangent, their are discontinuities in the phase (which
> is coming from arctan2(imag(vector), real(vector)).  For example, in
> the attached plot and data file, there is a point where the phase
> passes through -180 degrees.  arctan2 jumps to returning +180, +175
> etc, when it makes a pretty plot and makes more sense to me to convert
> +175 to -185.  What I had done previously is write a for loop that
> compares the current value of the phase to the value at the previous
> index (i.e. if phase[i]-phase[-1]>300 or if phase[i]-phase[i-1]<-300)
> to find the discontinuities.  Is there a cleaner/faster way to do
> this?  (I had done it in weave and will probably now do it in f2py for
> comparison on speeds).  I have attached/included a first attempt at an
> algorithm that is more vectorized.  But is there a fast, easy way to
> compare a vector to the same vector shift one element up or down?
> Right now I am doing this
> 
> phase2=r_[phase[0:1],phase[0:-1]]
> 
> to shift the vector (and reuse the first element so there isn't a jump
> from phase[-1] to phase[0] which may legitimately be far apart).
> 
> Has anyone else already done this?  Are there obvious bottlenecks with
> the algorithm?
> 
> Thanks for your thoughts,
> 
> Ryan
> 
> =====================================
> from scipy import arange, io, r_
> from pylab import show, figure, subplot, clf, cla, semilogx, xlabel,
> ylabel, legend, xlim
> 
> import copy
> 
> data=io.read_array('phase_data.txt')
> f=data[:,0]
> phase=data[:,1]
> figure(1)
> clf()
> startphase=copy.copy(phase)
> semilogx(f,startphase)
> 
> keepgoing=1
> 
> inds=arange(0,len(phase))
> 
> while keepgoing:
>     phase2=r_[phase[0:1],phase[0:-1]]#artificially repeat the first element
>     diff=phase-phase2
>     jumps=inds[abs(diff)>250]
>     keepgoing=jumps.any()
>     if not keepgoing:
>         break
>     ind=jumps.min()
>     if diff[ind]>0:
>         phase[ind:]-=360
>     else:
>         phase[ind:]+=360
> 
> 
> semilogx(f,phase,'r--')
> 
> legend(['Before unwrapping','after unwrapping'],2)
> 
> xlim([1,30])
> 
> show()
> _______________________________________________
> SciPy-user mailing list
> SciPy-user at scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-user