Instantaneous phase calculation problems

Question

I am currently struggling with a problem when calculating the instantaneous phase of a wavefield using a 1D Hilbert transformation.

The scheme is as follows.

1D Hilbert transformation of wave field $U(x,z)$ $$ q_z(x,z) = H_z[U(x,z)] = \int_R U(x-\xi,z)\frac{d\xi}{\xi} $$

Calculation of the instantaneous phase $\phi(x,z)$

$$ \phi(x,z) = arctan\left(\frac{q_z(x,z)}{U(x,z}\right), $$ where $x$ and $z$ are the coordinates and $H_z$ is the 1D Hilbert transformation over the vertical direction.

My problem is, that $\phi(x,z)$ which is implemented using the intrinsic function $atan2()$ contains phase shifts [$\pi \rightarrow -\pi$].

Using the instantaneous phase for further calculations in which its gradient is derived hence leads to imaging artifacts as shown in the attached pictures:

TOP LEFT: input wavefield $U(x,z)$

TOP RIGHT: Hilbert transformed wavefield $q_z(x,z) = H_z[U(x,z)]$

BOTTOM LEFT: instantaneous phase $\phi(x,z)$

BOTTOM RIGHT: derived propagation angle with artifacts

.

My question: I am able to overcome this problem using Gauss filtering or Median filtering. Yet this does not wotrk properly for more complicated wavefields and higher frequencies. Is there another way to overcome those sharp transitions in the area of my wavefield? (not the noise around)

Thanks alot!

Answer 1

One way to overcome this problem is phase unwrapping, which is also implemented by the Matlab function unwrap. But I believe that the best option for obtaining the derivative of the phase is to avoid a direct computation, but instead start out from the complex function itself. To simplify notation, let's use $U$ (without arguments):

$$U=re^{j\phi}=a+ib$$ where $r$ is the magnitude and $\phi$ is the phase. The derivative of $U$ is given by

$$U'=r'e^{j\phi}+rj\phi'e^{j\phi}\\ \frac{U'}{U}=\frac{r'}{r}+j\phi'$$

Since $r'/r$ is real-valued, we get

$$\phi'=\text{Im}\left\{\frac{U'}{U}\right\}=\frac{ab'-a'b}{a^2+b^2}\tag{1}$$

where $a$ and $b$ are the real and imaginary parts of $U$, respectively. Using (1) instead of a direct derivative of the phase should solve the problem.