Implementing continuous phase/Arg function

Question

In MATLAB there's a function called phase that is just like our Listable friend Arg, but that, when evaluated with a list, it tries to keep the result "continuous", allowing itself to return values outside the [-Pi, Pi] range. So, for example, sampling a complex exponential

vec = E^(I 2 π #/10) & /@ Range[100];
resMMA = Arg[vec] // N;

The red line is Mathematica's result with Arg, and the blue line is MATLAB's with phase

So, I want to implement a function that does a good job at allowing you to see phase properties such as "linearity", more easily, like phase, in Mathematica. Ideally, it should also work on some symbolic functions and not only lists, but what's more important to me is for it to work as good as possible on lists.

My attempt at the "list version" is very simple, but seems to be working reasonably, and fast. I also tested it with a function that's positive and negative, like vec2 = E^(I 2 Pi/16. #) Sinc[2 Pi #/20.] & /@ Range[100]; Here's the code

phase = Compile[{{l, _Complex, 1}}, 
   FoldList[
    Function[{prev, new}, # + Round[prev - #, 2 Pi] &@Arg@new], 
    Arg@First@l, Rest@l], 
   CompilationTarget -> "C", RuntimeOptions -> "Speed"];

Any improvements, would be appreciated. For example, it would probably be better if it could look at more than just the nearest neighbouring point, to be more robust to an occasional noise. I'd also appreciate ideas on how to achieve something analogous for a continuous well-behaved, not-created-by-mathematicians-just-to-prove-people-wrong symbolic function?

Answer 1

For a continuous function you could do something like this:

SetAttributes[argPlot, HoldAll];
Options[argPlot] = Options[Plot];

argPlot[exp_, {x_, x0_, x1_}, opt : OptionsPattern[argPlot]] :=
 Module[{pts, pl},
  pl = Plot[Arg[exp], {x, x0, x1}, PlotRange -> All, 
    PlotPoints -> OptionValue[PlotPoints]];
  pts = SortBy[Cases[pl, Line[pts_] :> pts, Infinity], #[[1, 1]] &];
  pts = Reap[Fold[Module[{ptsn},
         ptsn = #2;
         ptsn[[All, 2]] -= Round[ptsn[[1, 2]] - #1, 2 Pi];
         Sow[ptsn];
         ptsn[[-1, 2]]] &, 0, pts];][[2, 1]];      
  ListLinePlot[Flatten[pts, 1], opt]]

argPlot[3 + 2 Exp[3 I a] + Exp[(1 - I a^2)], {a, -8, 8}]

Compared to an ordinary plot of the Arg[3 + 2 Exp[3 I a] + Exp[(1 - I a^2)]]

Plot[Arg[3 + 2 Exp[3 I a] + Exp[(1 - I a^2)]], {a, -8, 8}, PlotRange -> All]

Answer 2

This solution only considers direct differences, but is a bit shorter than the original (it doesn't contain compilation, though, and I didn't test the speed):

phase[l: {_?NumericQ ..}] :=
  Module[{args = Arg[l]},
         args+Prepend[2Pi Accumulate@-IntegerPart@Differences[args/Pi],0]]

The following solution definitely isn't the fastest, but generally should give the best possible result:

phase[l: {_?NumericQ ..}]:=
  With[{vars = Table[Unique[], {Length@l}]},
       Arg@l+2 Pi Accumulate@vars /.
         Last@NMinimize[{Total@Abs@Differences[Arg@l+2\[Pi] Accumulate@vars],
                         (Alternatives@@vars) \[Element] Integers &&
                         First@vars==0},
                        vars]]

Answer 3

We can transfer the phase function from Matlab to Mathematica:

phase[vec_List] := 
  Module[{phi, df, len, i},
    phi = Arg @ vec;
    df = Differences @ phi;
    len = Length @ phi;
    i = Flatten @ Position[df, x_ /; Abs[x] > 3.5];
    Do[phi = phi - (2 Pi*Sign[df[[j]]]*UnitStep[# - (j + 1)] & /@ Range[len]),
      {j, i}]; 
    phi]

data = Table[1/((3/2 Exp[I a])^4 - 1), {a, -8, 8, 0.01}];
data // phase // ListLinePlot

The output is

while

data // Arg // ListLinePlot

gives

Answer 4

Here is a simple answer for an ordinary list. The OP asked about:

vec = E^(I 2 \[Pi] #/10) & /@ Range[100];
resMMA = Arg[vec] // N;

Just a simple few functions is enough:

unwraped =   resMMA + Accumulate[2*Pi*Unitize[Threshold[Prepend[Differences[resMMA], 0.], 1]]];
ListLinePlot[{resMMA, unwraped}]

Sometimes you might need to subtract instead of add, if your saw-tooth is upside down.