相位展开(phase unwrapping)

我在做科研项目时，遇到过使用相位信息，但是相位被折叠(wrapped)的情况，相位计算公式如下：
$$\Theta =arctan(\frac{b}{a})$$
其中b、a分别代表虚部和实部，由于反正切函数的区间是$(-\frac{\pi }{2},\frac{\pi }{2})$，故实际计算相位时会将结果折叠在这个区间内。如下图，这张图使用的是matlab的atan函数计算并绘制的，可见，相位范围被折叠在$(-\frac{\pi }{2},\frac{\pi }{2})$区间内，且在某些点出现了跳跃。

仔细观察跳跃点，不难发现，其跳跃要么是$+\pi$，要么是$-\pi$，通过前后点的差值可以判断是否发生跳跃。现在为了说明怎么相位展开，先作以下符号说明：$x_w(n)$表示被折叠的相位信号，$x_u(n)$表示展开的相位信号，则相位展开步骤如下：

1. 令$x_u(n)$ = $x_w(n)$
2. 从$x_w(n)$的第二个点开始，计算当前点与前一个点差
3. 若差大于$\frac{\pi }{2}$，则$x_u(n)$当前点及后续所有点减去$\pi$
4. 若差小于$-\frac{\pi }{2}$，则$x_u(n)$当前点及后续所有点加上$\pi$
5. 重复上述步骤，直至遍历$x_u(n)$
   对应的matlab代码如下：

xu = xw;
for i = 2:length(xw)
        diff = xw(i) - xw(i-1);
        if diff > pi/2
            xu(i:end) = xu(i:end) - pi;
        elseif diff < -pi/2
            xu(i:end) = xu(i:end) + pi;
        end
end

展开后相位如下图所示，可见展开后的相位是连续的。

另：若使用atan2计算相位，上述代码需做一些改变:所有的$\pi$–>$2\pi$，所有的$\frac{\pi }{2}$–>$\pi$，具体如下代码所示。

xu = xw;
for i = 2:length(xw)
        diff = xw(i) - xw(i-1);
        if diff > pi
            xu(i:end) = xu(i:end) - 2*pi;
        elseif diff < -pi
            xu(i:end) = xu(i:end) + 2*pi;
        end
end