相位展开(phase unwrapping)

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

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

1. 令$x_{u}(n)$xu​(n) = $x_{w}(n)$xw​(n)
2. 从$x_{w}(n)$xw​(n)的第二个点开始，计算当前点与前一个点差
3. 若差大于$\frac{\pi}{2}$2π​，则$x_{w}(n)$xw​(n)当前点及后续所有点减去$\pi$π
4. 若差小于$-\frac{\pi}{2}$−2π​，则$x_{w}(n)$xw​(n)当前点及后续所有点加上$\pi$π
5. 重复上述步骤，直至遍历$x_{u}(n)$xu​(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$2π，所有的$\frac{\pi}{2}$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