5
$\begingroup$

The problem of unwrapping 1D phase often comes up in Fourier analysis: someone computes a DFT of some data, the resulting phase values fall between $-\pi$ and $\pi$, and they'd like to "unwrap" these phase values to not be bounded to (-$\pi$,$\pi$]. It's fairly straightforward in the case of low noise and sufficient discrete sampling, as seen in the excellent answers here.

However, by working more with the continuous case and/or no noise, the essential problems typically associated with 1D phase unwrapping are not dealt with. Here's an example of where things can go wrong in a discrete case with some noise and limited sampling.

phaseUnwrap[vec_?VectorQ] := (*translation of Matlab's implementation*)
 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];

data = Import["http://pastebin.com/raw/haAxt0br", "Package"];

ListPlot[phaseUnwrap[data[[;; , 2]]], Joined -> True]

phase plot

Even with noise and low-ish sampling the algorithm does a pretty good job, but there are still some obvious jumps in the phase unwrapping process, notably just after point 300 and just before point 500.

Given how easy it is to pick these out by eye, it seems like there could be a simple way to correct them. Maybe minimizing the derivative somehow? Checking afterward for deviations from a linear fit?

$\endgroup$
  • $\begingroup$ @J. M. I refocused the question a little to attempt to distinguish it from the other phase-unwrapping post. My thinking is that for the discrete case with low sampling and non-zero noise, this becomes a different sort of problem, and one that is not addressed in the other post. If it's still too similar perhaps I can synthesize my question and answer and move it to the other post. $\endgroup$ – nadlr Feb 2 '17 at 19:01
  • $\begingroup$ Perhaps you could make a smooth line through the noisy data in the complex plane. f = BSplineFunction[ReIm[data[[All, 2]]]]; ListLinePlot[phaseUnwrap[Table[{1, I}. f[x], {x, 0, 1, 1/671}]]] $\endgroup$ – Simon Woods Feb 2 '17 at 22:00
  • $\begingroup$ It looks sufficiently different now. Thankfully it's been reopened in the meantime. $\endgroup$ – J. M. will be back soon Feb 3 '17 at 4:25
3
$\begingroup$

Alright, I found a decent way to do this given some assumptions about what the shape of the data "should" be. This is a post-processing function, i.e., it should be applied after a function like phaseUnwrap. Essentially I assume that the data should be linear, and further that in some part of the domain there is already a nice linear region. As input we take the lower and upper bounds of this region, as well as the messy unwrapped data.

magnetize[phaseData_, lower_, upper_] :=
 Module[{fitdata, fit, magnet, magnetizedData, tries},
  fitdata = Transpose[{
 Range[lower, upper],
 phaseData[[lower ;; upper]]}];
 fit = LinearModelFit[fitdata, x, x];
 magnet = First@Last@Reap@
      For[i = 1, i < (Length@phaseData + 1), i++,
   tries = Sort[
     Table[{k, 
       Abs[(fit[i] - phaseData[[i]] - (2 Pi k))]}, {k, -10, 10, 1}],
     #1[[2]] > #2[[2]] &];
   Sow[
    tries[[-1, 1]]]];
  magnetizedData = phaseData + 2 Pi magnet;
  {magnetizedData, fit,
 Show[{ListPlot[{phaseData, magnetizedData}, Joined -> True, 
  GridLines -> {Automatic, Automatic},
  Epilog -> {
    Line[{{lower, 0}, {lower, 300}}],
    Line[{{upper, 0}, {upper, 300}}]}],
 Plot[fit[x], {x, 0, Ceiling@(1.1*Length@phaseData)},
  PlotStyle -> Directive[Gray]]
 }]}
]

In summary,

  1. Find a linear fit line based on a "nice" part of the data
  2. For every point, find the integer multiple of $2\pi$ that brings the data closest to the fit line.
  3. Return {processed data, fit line object, plot}

It's a little rough but it gets the job done.

mag1 = magnetize[phaseUnwrap[data[[;; , 2]]], 50, 120];   
(*using data and function from the question post *) 
mag1[[3]]

image of plot output

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.