相位解卷绕-UNwrap

由于工作中涉及到相位求解及解卷绕，而且需要进行算法的上板，因此花了一点时间分析了MATLAB的UNwrap算法并基于此写了一个C语言版本，现做一个整理及分享。

UNwrap详解

一般在我们计算一个系统相频特性时，就要用到反正切函数提取相位，计算机中反正切函数规定，在一、二象限中的角度为0～pi，三四象限的角度为0～-pi。
反正切函数一共有两个：atan(double s)和atan2(double y, double x)，二者均返回弧度值，如果转换成角度需要额外计算。

atan(double s)：可以理解为二象限反正切，返回正切值double s = y/x所对应的弧度值，但是取值仅在$[-\frac{\pi }{2},\frac{\pi }{2}]$
atan2(double y, double x)：是四象限反正切，返回点(x,y)对应的弧度值，但是取值在$[-\pi ,\pi ]$之间

相位卷绕及解卷绕

在实际的应用处理过程中，反正切函数所提取的相位可能会发生跳变，跳变幅度为2π，这就是相位的卷绕。unwrap的作用就是解卷绕，使相位在π处不发生跳变，从而展示出真实的相位变化。
在实际应用过程中，UNwrap在检查到数据前后两点的差距超过π时就认为有跳变，即需要进行解卷绕。
举例模拟：
前后差值为：-0.5π，0，0.8π，跳变为0；
前后差值为：1.1π，1.9π，2.5π，跳变为（-2π），修正后为-0.9π，-0.1π，0.5π；
前后差值为：-3.1π，-3.8π，-4.6π，跳变为（2π*2），修正后为0.9π，0.2π，-0.6π；
…
过程拆分解析：

1. 前后两个数据的差值在[-π，π]之间，则跳变的大小为0；
2. 前后两个数据的差值大于π，则跳变大小为-2π*n，n为跳变的次数；
3. 前后两个数据的差值小于π，则跳变大小为2π*n，n为跳变的次数；
   基于上述分析我们可以发现跳变规律是按照前后差值进行计算的，且该变化规律是以0值为对称中心，只是跳变次数需要注意取负号。

规律总结：
4. (0，π]： 跳变为0 * 2π； (-π，0]： 跳变为0 * 2π；
5. (π，3π]： 跳变为-1 * 2π； (-3π，-π]： 跳变为1 * 2π；
6. (3π，5π]：跳变为-2 * 2π ； (-5π，-3π]：跳变为2 * 2π；
7. (5π，7π]：跳变为-3 * 2π； (-7π，-5π]：跳变为3 * 2π；
8. …
9.
算法设计思路：
可以将前后两个数的差值与2π相除得到的结果来确定跳变的次数，例如：
案例1. 差值为1.1π，其与2π相除的结果是0.55，整数部分是0，余数大于等于0.5则+1，最终跳变为1；
案例2. 差值为2.5π，其与2π相除的结果是1.25，整数部分是1，余数小于0.5则+0，最终跳变为1；
案例3. 差值为-3.1π，其与2π相除的结果是-1.55，整数部分是-1，余数大于0.5则-0，最终跳变为-2；

算法流程：
10. 计算实时相位差dp= phase_new - phase_prior；
11. 计算跳变跳变幅度dp_corr（即规律2，3，4，…）；
2.1. 计算dp和2π的比值dp_corr，该比值的整数部分即需要跳变的幅度值；
2.2. 计算dp_corr和1之间的余数，利用该余数的绝对值和0.5进行比较来修正跳变次数；
2.3. 大于0.5则四舍五入进一，小于0.5则截尾取整；
12. 利用实时相位差dp和阈值π比较（即规律0），再一次修正跳变幅度dp_corr；
13. 将实时相位减去2π*跳变幅度dp_corr进行相位修正（dp_corr需要进行累加，因为前一时刻如果需要进行跳变修改，那么后续都需要进行跳变修改）。

阈值为π的相位解卷绕C代码实现

全局变量设置，随着接收到新的数据不断更新

double                MssCusum = 0.0;           //unwrap中的累积和
double                MssPriorPhase = 0.0;      //UNwrap中的前一个相位数据，初始化为0

阈值为π的相位解卷绕

/* *
 *  @b Description
 *  @n
 *      Author：柠の
 *      Myself UNWrap
 *
 *  @Param[in]
 *      phase computed in real time
 *
 *  @retval
 *      Always return correct real-time result after Unwrap.
 */
double MmDemo_unwrap(double Phase)
{
    double       dp;          // phase variations
    double       dp_corr;
    double       dp_corr_rem; // remainder of dp_corr
    double       PhaseLatter; // Phase behind unwrap
    bool         roundDown;   //

    dp              =      Phase - MssPriorPhase;
    MssPriorPhase   =      Phase;

    /* Compute an integer describing how many times 2*pi we are off */
    // dp in [-pi, pi]: dp_corr = 0,
    // elseif dp in [-3*pi, 3*pi]: dp_corr = 1,
    // else if dp in [-5*pi, 5*pi]: dp_corr = 2, ...

    dp_corr         =      dp / 2 /M_PI;
    dp_corr_rem     =      fmod(dp_corr,1.0);
    if ( fabs(dp_corr_rem) <= 0.5 )
    {
        roundDown = 1;
    }
    else
    {
        roundDown = 0;
    }

    if (roundDown)                 //   如果小于等于0.5，那么将dp_corr截尾取整
    {
        if (dp_corr >= 0)
        {
            dp_corr = floor(dp_corr);
        }
        else
        {
            dp_corr = ceil(dp_corr);
        }
    }

    dp_corr = round(dp_corr);      // 四舍五入

    if (fabs(dp) < M_PI)           //unwrap的跳变阈值，如果需要自适应修改，将该参数转为函数的输入
    {
        dp_corr = 0;
    }

    MssCusum        =       MssCusum + dp_corr;
    PhaseLatter     =       Phase - (2*M_PI)*MssCusum;
    return PhaseLatter;
}

总结

上述就是我参考MATLAB的UNwrap算法进行修改后编写的C语言版的UNwrap算法，经过简单测试是没有问题的，希望各位批评指正。