diff options
Diffstat (limited to 'macros/unwrap2.sci')
| -rw-r--r-- | macros/unwrap2.sci | 244 |
1 files changed, 210 insertions, 34 deletions
diff --git a/macros/unwrap2.sci b/macros/unwrap2.sci index b8ea9de..7844e38 100644 --- a/macros/unwrap2.sci +++ b/macros/unwrap2.sci @@ -1,35 +1,211 @@ -function Y = unwrap2 (X, TOL, DIM) -//Unwrap radian phases by adding or subtracting multiples of 2*pi. -//Calling Sequence -//B = unwrap(X) -//B = unwrap(X, TOL) -//B = unwrap(X, TOL, DIM) -//Parameters -//Description -//This function unwraps radian phases by adding or subtracting multiples of 2*pi as appropriate to remove jumps greater than TOL. -// -// TOL defaults to pi. -// -//Unwrap will work along the dimension DIM. If DIM is unspecified it defaults to the first non-singleton dimension. -//Examples -//unwrap2([1,2,3]) -//ans = -// 1. 2. 3. -funcprot(0); -lhs = argn(1); -rhs = argn(2); -if (rhs < 1 | rhs > 3) -error("Wrong number of input arguments."); -end - -select(rhs) - - case 1 then - Y = callOctave("unwrap",X); - case 2 then - Y = callOctave("unwrap",X,TOL); - case 3 then - Y = callOctave("unwrap",X,TOL,DIM); - end - +/*Description: + The unwrap function adjusts radian phases in the input array x by adding or subtracting multiples of + 2π as necessary to remove phase jumps that exceed the specified tolerance tol. If tol is not provided, it defaults to 𝜋 + Radian Phases: These are typically angles or phases expressed in radians, commonly encountered in signal processing and communication systems. + Tolerance (tol): Determines the maximum allowable discontinuity in the phases. + If the difference between consecutive phases exceeds tol, unwrap adjusts the phase by adding or subtracting 2π. + Dimension (dim): Specifies the dimension along which the unwrapping operation is applied. + By default, unwrap operates along the first non-singleton dimension of the input array x. +Calling Sequence: + b = unwrap(x) + b = unwrap(x, tol) + b = unwrap(x, tol, dim) +Parameters: + x: Input array containing radian phases to be unwrapped. + tol (optional): Tolerance parameter specifying the maximum jump allowed between consecutive phases before adding or subtracting 2π. Defaults to 𝜋 + dim (optional): Dimension along which to unwrap the phases. If unspecified, dim defaults to the first non-singleton dimension of the array x. +Dependencies : ipermute*/ +function retval = unwrap2 (x, tol, dim) + nargin = argn(2) + if (nargin < 1) + error("invalid number of inputs"); + end + if (~ (type(x) == [ 1 5 8]) || or(type(x)==[4,6])) + error ("unwrap2: X must be numeric"); + end + if (nargin < 2 || isempty (tol)) + tol = %pi; + end + // Don't let anyone use a negative value for TOL. + tol = abs (tol); + nd = ndims (x); + sz = size (x); + if (nargin == 3) + if (~(or(type(dim)==[1 5 8])&& isscalar (dim) && ... + dim == fix (dim)) || ~(1 <= dim)) + error ("unwrap2: DIM must be an integer and a valid dimension"); + end + else + // Find the first non-singleton dimension. + dim = find (sz > 1, 1) + if isempty(dim) + dim = 1; + end + end + rng = 2*%pi; + // Handle case where we are trying to unwrap a scalar, or only have + // one sample in the specified dimension (a given when dim > nd). + if ((dim > nd) || ( sz(dim) == 1)) + retval = x; + return; + end + if (and(abs(x(:))<%inf ) ) + // Take first order difference so that wraps will show up as large values + // and the sign will show direction. + if length(sz) < 3 + sz(3) = 1 ; + end + sz(dim) = 1; + zero_padding = zeros (sz(1),sz(2),sz(3)); + d = cat (dim, zero_padding, -diff (x, 1, dim)); + // Find only the peaks and multiply them by the appropriate amount + // of ranges so that there are kronecker deltas at each wrap point + // multiplied by the appropriate amount of range values. + p = round (abs (d)./rng) .* rng .* (double(double(d > tol) > 0) - double(double(d < -tol) > 0)); + // Integrate this so that the deltas become cumulative steps to shift + // the original data. + retval = cumsum (p, dim) + x; + else + // Unwrap needs to skip over NaN, NA, Inf in wrapping calculations. + if (isvector (x)) + // Simlpified path for vector inputs. + retval = x; + xfin_idx = abs(x)<%inf ; + xfin = x(xfin_idx); + d = cat (dim, 0, -diff(xfin, 1, dim)); + p = round (abs (d)./rng) .* rng .* (double(double(d > tol) > 0) - double(double(d < -tol) > 0)); + retval(xfin_idx) = xfin + cumsum (p, dim); + else + // For n-dimensional arrays with a possibly unequal number of non-finite + // values, mask entries with values that do not impact calcualation. + // Locate nonfinite values. + nf_idx = ~ abs(x)<%inf; + if (and(nf_idx(:))) + // Trivial case, all non-finite values + retval = x; + return; + end + // Permute all operations to occur along dim 1. Inverse permute at end. + permuteflag = dim ~= 1; + if (permuteflag) + perm_idx = [1 : nd]; + perm_idx([1, dim]) = [dim, 1]; + x = permute (x, perm_idx); + nf_idx = permute (nf_idx, perm_idx); + sz([1, dim]) = sz([dim, 1]); + dim = 1; + end + // Substitute next value in dim direction for nonfinite values(ignoring + // any at trailing end) to prevent calculation impact. + x_nf = x(nf_idx); // Store nonfinite values. + zero_padding = zeros ([1, sz(2:$)]); + x = __fill_nonfinite_columnwise__ (x, nf_idx, zero_padding, sz, nd); + d = [zero_padding; -diff(x, 1, 1)]; + p = round (abs (d)./rng) .* rng .* ... + (((d > tol) > 0) - ((d < -tol) > 0)); + retval = x + cumsum (p, 1); + // Restore nonfinite values. + retval(nf_idx) = x_nf; + // Invert permutation. + if (permuteflag) + retval = ipermute (retval, perm_idx); + end + end + end endfunction +function y = repelems(x,r) + y = []; + for i = 1:size(r,2) + y = [y, x(r(1,i)*ones(1, r(2,i)))]; + end +endfunction +function x = __fill_nonfinite_columnwise__ (x, nonfinite_loc, zero_padding, szx, ndx) + // Replace non-finite values of x, as indicated by logical index + // nonfinite_loc, with next values. + flip_idx(1:ndx) = {':'}; + flip_idx(1) = {szx(1):-1:1}; + // Isolate nf values by location: + nf_front = cumprod (nonfinite_loc, 1); + nf_back = cumprod (nonfinite_loc(flip_idx{:}), 1)(flip_idx{:}); + nf_middle = nonfinite_loc & ~ (nf_back | nf_front); + // Process bound/middle elements + locs_before = [diff(nf_middle, 1, 1); zero_padding] == 1; + locs_after = diff ([zero_padding; nf_middle], 1, 1) == -1; + mid_gap_sizes = find (locs_after) - find (locs_before) - 1; + x(nf_middle) = repelems (x(locs_after), ... + [1 : numel(mid_gap_sizes); mid_gap_sizes'])'; + // Process front edge elements + nf_front = nf_front & ~ and (nonfinite_loc, 1); // Remove all nf columns. + locs_after = diff ([zero_padding; nf_front], 1, 1) == -1; + front_gap_sizes = (sum (nf_front, 1))(any (nf_front, 1))(:); + x(nf_front) = repelems (x(locs_after), ... + [1:length(front_gap_sizes); front_gap_sizes'])'; +endfunction +/* +//Test cases +i = 0; +t = []; +r = [0:100]; // original vector +w = r - 2*%pi*floor ((r+%pi)/(2*%pi)); // wrapped into [-pi,pi] +tol = 1e3*%eps; +assert_checkalmostequal (r, unwrap2 (w), tol) +assert_checkalmostequal (r', unwrap2 (w'), tol) +assert_checkalmostequal ([r',r'], unwrap2 ([w',w']), tol) +assert_checkalmostequal ([r; r ], unwrap2 ([w; w ], [], 2), tol) +assert_checkalmostequal(r + 10, unwrap2 (10 + w), tol) +assert_checkequal (w', unwrap2 (w', [], 2)) +assert_checkequal(w, unwrap2 (w, [], 1)) +assert_checkequal([w; w], unwrap2 ([w; w])) +// Test that small values of tol have the same effect as tol = pi +assert_checkalmostequal(r, unwrap2 (w, 0.1), tol) +assert_checkalmostequal(r, unwrap2 (w, %eps), tol) +// Test that phase changes larger than 2*pi unwrap properly +assert_checkalmostequal([0; 1], unwrap2([0; 1])) +assert_checkalmostequal([0; 4 - 2*%pi], unwrap2 ([0; 4])) +assert_checkalmostequal([0; 7 - 2*%pi], unwrap2 ([0; 7])) +assert_checkalmostequal([0; 10 - 4*%pi], unwrap2 ([0; 10])) +assert_checkalmostequal([0; 13 - 4*%pi], unwrap2 ([0; 13])) +assert_checkalmostequal([0; 16 - 6*%pi], unwrap2 ([0; 16])) +assert_checkalmostequal([0; 19 - 6*%pi], unwrap2 ([0; 19])) +//test +A = [%pi*(-4), %pi*(-2+1/6), %pi/4, %pi*(2+1/3), %pi*(4+1/2), %pi*(8+2/3), %pi*(16+1), %pi*(32+3/2), %pi*64]; +assert_checkalmostequal (unwrap2 (A), unwrap2 (A, %pi)); +assert_checkalmostequal (unwrap2 (A, %pi), unwrap2 (A, %pi, 2)); +assert_checkalmostequal (unwrap2 (A', %pi), unwrap2 (A', %pi, 1)); +//test +A = [%pi*(-4); %pi*(2+1/3); %pi*(16+1)]; +B = [%pi*(-2+1/6); %pi*(4+1/2); %pi*(32+3/2)]; +C = [%pi/4; %pi*(8+2/3); %pi*64]; +D = [%pi*(-2+1/6); %pi*(2+1/3); %pi*(8+2/3)]; +E(:, :, 1) = [A, B, C, D]; +E(:, :, 2) = [A+B, B+C, C+D, D+A]; +F(:, :, 1) = [unwrap2(A), unwrap2(B), unwrap2(C), unwrap2(D)]; +F(:, :, 2) = [unwrap2(A+B), unwrap2(B+C), unwrap2(C+D), unwrap2(D+A)]; +assert_checkalmostequal (unwrap2 (E), F); +// Test trivial return for m = 1 and dim > nd +assert_checkalmostequal (unwrap2 (ones(4,1), [], 1), ones(4,1)) +assert_checkalmostequal (unwrap2 (ones(4,1), [], 2), ones(4,1)) +assert_checkalmostequal (unwrap2 (ones(4,1), [], 3), ones(4,1)) +assert_checkalmostequal (unwrap2 (ones(4,3,2), [], 99), ones(4,3,2)) +// Test empty input return +assert_checkalmostequal (unwrap2 ([]), []) +assert_checkalmostequal (unwrap2 (ones (1,0)), ones (1,0)) +assert_checkalmostequal (unwrap2 (ones (1,0), [], 1), ones (1,0)) +assert_checkalmostequal (unwrap2 (ones (1,0), [], 2), ones (1,0)) +assert_checkalmostequal (unwrap2 (ones (1,0), [], 3), ones (1,0)) +// Test trivial return for m = 1 and dim > nd +assert_checkalmostequal (unwrap2 (ones(4,1), [], 1), ones(4,1)) +assert_checkalmostequal (unwrap2 (ones(4,1), [], 2), ones(4,1)) +assert_checkalmostequal (unwrap2 (ones(4,1), [], 3), ones(4,1)) +assert_checkalmostequal (unwrap2 (ones(4,3,2), [], 99), ones(4,3,2)) +// Test empty input return +assert_checkalmostequal (unwrap2 ([]), []) +assert_checkalmostequal (unwrap2 (ones (1,0)), ones (1,0)) +assert_checkalmostequal (unwrap2 (ones (1,0), [], 1), ones (1,0)) +assert_checkalmostequal (unwrap2 (ones (1,0), [], 2), ones (1,0)) +assert_checkalmostequal (unwrap2 (ones (1,0), [], 3), ones (1,0)) +// Test handling of non-finite values +x = %pi * [-%inf, 0.5, -1, %nan, %inf, -0.5, 1]; +assert_checkalmostequal (unwrap2 (x), %pi * [-%inf, 0.5, 1, %nan, %inf, 1.5, 1], %eps) +assert_checkalmostequal (unwrap2 (x.'), %pi * [-%inf, 0.5, 1, %nan, %inf, 1.5, 1].', %eps) +*/
\ No newline at end of file |