Re: Unwrap a function

[lash]

> 
> There is under Matlab a very useful function which is called unwrap 
> and ... which unwraps (:->) a set of data (for example a phase in
> signal processing).
> 
> Octave does not seem to hold this function (:-<)
> 
> Could-you help me ?
> 
> Surf-you-soon...
> 
> Emmanuel
> 
> address@hidden
> 
> 
> 

Oops, I meant to send this to the list, but only sent it to Emmanuel.

I wrote a version over lunch, and then tidied it up a bit later.  I
don't have access to Matlab, so I don't know if it uses the same
calling conventions or not.  It should be pretty easy to modify.


Any suggestions, improvements, or criticisms of the code are welcome.

Bill Lash
address@hidden


## Usage: b = unwrap(a{,range{,threshold}})
## 
## Unwraps a set of numbers a, (e.g. phase data)
## 
## If only one argument is given, the range of the data is
## determined by finding the max and min values of the entire
## data set.  The wrapping function will then occur when a change
## between any two adjacent samples is greater than half this range.
## when a wrap is detected, the direction is determined, and the
## full range is added or subtracted as appropriate.  If a range
## is given, a threshold can also be specified.  The threshold
## is a number between 0 and 1 that indicates how large (as a 
## proportion of the range) a jump between samples must be.
##
## It is recommended that if the range of the data is known,
## that it be supplied, especially if the the elements of a are
## thinly sampled.  For phase data, setting a range of pi or
## 2pi as appropriate may be needed.
##
## example for testing that the function works:
##
##      m = [0:200,200:-1:0]/10;   # make some data that slopes up, 
##                                 # then down
##      r = tan(m);                # take the tangent 
##      n = unwrap(atan(r));       # then unwrap the arctan letting 
##                                 # unwrap find the range
##      diff = max(abs(n-m))       # Not bad, max err is about 0.1
##      n_pi = unwrap(atan(r),pi); # now specify that the range is
##                                 # pi
##      diff_pi = max(abs(n_pi-m)) # Now the error is really small
##      
##
## Bill Lash
## address@hidden

function retval = unwrap(a,range,thresh)
        
        if ((nargin<1) || (nargin>3))
                usage("unwrap(a,[range,[threshold]])")
        endif
        if (nargin < 3)
                thresh = 0.5;
        endif
        if (nargin < 2)
                range = max(a) - min(a);
        endif

        #
        # want a row vector to work on
        #

        ra = reshape(a,1,columns(a)*rows(a));

        #
        # take first order differnce to see so that wraps will show up
        # as large values, and the sign will show direction
        #
        diff = filter([-1,1],1,ra);

        #
        # Find only the peaks, and multiply them by the range so that there
        # are kronecker deltas at each wrap point multiplied by the range
        # value
        #
        p = range * (((diff > thresh*range)>0) - ((diff < -thresh*range)>0));

        #
        # Now need to "integrate" this so that the deltas become steps
        #
        r = filter(1,[1,-1],p);

        #
        # Now add the "steps" to the original data and reshape to the original
        # shape
        #
        retval = reshape(ra+r,rows(a),columns(a));

endfunction