# Implementing continuous phase/Arg function

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?

-
Indeed, very useful. I wonder how to implement similar thing for Plot3D. It'd be needed for example when solving a Schrodinger equation and plotting Arg of wave-function versus time and space. – Vitaliy Kaurov May 19 '12 at 9:16
For completeness, here is a straightforward Mathematica translation of the code used within MATLAB: phase[vec_?VectorQ] := Module[{n = Length[vec], ph = Arg[vec], df}, df = Differences[ph]; Do[ph -= 2 Pi Sign[df[[k]]] PadRight[ConstantArray[0, k], n, 1], {k, Flatten[Position[Abs[df], _?(# > Pi &)]]}]; ph] – J. M. May 19 '12 at 16:11

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]


-
+1 Very nice ;-) – Vitaliy Kaurov May 19 '12 at 9:13
+1 Greeeeeeeeat – Rojo May 19 '12 at 16:35

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]]

-

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

-