# 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;
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; 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. May 19, 2012 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] May 19, 2012 at 16:11
• forums.wolfram.com/mathgroup/archive/1998/May/msg00107.html May 20, 2012 at 16:26

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] • Nice test function. Aug 10, 2019 at 7:46
• This doesn't work if the phase shift is around the origin, try: argPlot[50/((s + 4) (s + 3) (s - 1)) + 1 /. s -> \[ImaginaryJ] \[Omega] /. \[Omega] -> Tan[x], {x, 0, \[Pi]/2}] May 28, 2021 at 5:26

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

• Just a side note, the following is the compiled version of the first function: Clear[phase]; phase = Compile[{{l, _Complex, 1}}, Module[{args = Arg[l]}, args + Prepend[2 Pi Accumulate@-IntegerPart@Differences[args/Pi], 0]], RuntimeOptions -> "EvaluateSymbolically" -> False];  Jun 11, 2018 at 19:01

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 Here is a simple answer for an ordinary list. The OP asked about:

vec = E^(I 2 \[Pi] #/10) & /@ Range;
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.

• I'd suggest to use Threshold instead of Chop because the former preserves packing of lists. And using 0. instead of 0. Aug 9, 2019 at 18:12
• This code doesn't work on this dataset (which I import by resMMA = Flatten@Import["/tmp/angles.csv"];). Jan 21, 2020 at 10:45