;+ ; NAME: ; PHUNWRAP ; ; AUTHOR: ; Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770 ; craigm@lheamail.gsfc.nasa.gov ; UPDATED VERSIONs can be found on my WEB PAGE: ; http://cow.physics.wisc.edu/~craigm/idl/idl.html ; ; PURPOSE: ; Unwrap phase jumps to recover cycle counts ; ; MAJOR TOPICS: ; Mathematics ; ; CALLING SEQUENCE: ; CYCLES = PHUNWRAP(PHASE, TOLERANCE=, MAXVAL=) ; ; DESCRIPTION: ; ; PHUNWRAP unwraps a sequence of phases to produce a new series of ; cycle counts. Phase jumps due to crossing over the PHASE=0 ; boundary are removed by adding an integral number of cycles. The ; algorithm is based on the MATLAB "unwrap" function. ; ; NOTE: the unwrapping process can be ambiguous if there is a phase ; jump of more than a half cycle in the series. For example, if the ; phase changes by ~0.5 cycles, it is not possible to distinguish ; whether there wasa +0.5 cycle or -0.5 cycle jump. The most ; accurate unwrapping can be performed if the PHASE series is nearly ; continuous and does not have rapid phase changes. ; ; Users can select the tolerance used to determine the phase jump. ; Users can also select the definition of "1 cycle" by changing ; MAXVAL. By default, MAXVAL is 2*!DPI, which correspondes to 1 ; cycle = 2*!DPI radians, but other values of 1 (cycle), or 360 ; (degrees) are possible. ; ; INPUTS: ; ; PHASE - phase series to be unwrapped. Values should range from 0 ; to MAXVAL. The ordering of the series is important. ; ; RETURNS: ; ; A new series, expressed in cycles, with cycle jumps larger than ; TOLERANCE removed. ; ; OPTIONAL KEYWORDS: ; ; TOLERANCE - phase jump tolerance. If the phase from one sample to ; the next changes by more than TOLERANCE, then a single ; cycle jump is assumed to have occurred. ; DEFAULT: 0.5*MAXVAL ; ; MAXVAL - Maximum value for phase. Common values are: 2*!DPI ; (radians; DEFAULT); 1 (cycle); 360 (degrees), but any ; positive value may be used. ; ; EXAMPLE: ; ; ;; Set up some fake data ; x = dindgen(100)/10d ; y = x/2 ; ph = y MOD 1.0 ;; Mock phases ; ; cycles = phunwrap(ph, maxval=1) ; ; MODIFICATION HISTORY: ; Written and documented, CM, July 2003 ; Handle the case of unsigned integer input, CM, Feb 2006 ; ; $Id: phunwrap.pro,v 1.3 2006/03/28 14:19:53 craigm Exp $ ; ;- ; Copyright (C) 2003, 2006, Craig Markwardt ; This software is provided as is without any warranty whatsoever. ; Permission to use, copy, modify, and distribute modified or ; unmodified copies is granted, provided this copyright and disclaimer ; are included unchanged. ;- function phunwrap, ph, tolerance=tol0, maxval=maxval0 common phunwrap_common, idlver if n_elements(idlver) EQ 0 then begin idlver = !version.release endif if n_elements(maxval0) EQ 0 then maxval = 2d*!dpi else maxval = maxval0(0) if n_elements(tol0) EQ 0 then tol = 0.5*maxval else tol = tol0(0)*maxval if n_elements(ph) LT 2 then return, ph sz = size(ph) tp = sz(sz(0)+1) ;; First order difference case tp of 12: dph = [0, long(ph)-long(ph(1:*))] 13: dph = [0, long64(ph)-long64(ph(1:*))] 15: dph = [0, long64(ph)-long64(ph(1:*))] else: dph = [0, ph - ph(1:*)] endcase p = maxval * (fix((dph GT tol) EQ 1) - fix((dph LT (-tol)) EQ 1)) if idlver GT 5.25 then begin ;; Use built-in version if available r = total(p, /cumulative) endif else begin ;; .. if not, then use the lame FOR loop r = p for i = 1L, n_elements(r)-1 do $ r(i) = r(i) + r(i-1) endelse return, ph+r end