相位展开(解卷绕，unwrap)的Python实现

究极小菜鸟

已于 2024-01-09 18:24:35 修改

阅读量1.6k

点赞数 11

文章标签： python 开发语言

于 2024-01-09 18:20:43 首次发布

本文链接：https://blog.csdn.net/King_Chanzy/article/details/135485767

版权

0 前言

由于求余弦函数的反函数得到的相位值总是在一个固定的区间内，也就是说反函数的值域在[0,π]之间，使得求出的相位值也是周期性变化的，所以需要进行相位展开或者是解卷绕的操作。

比如随便画一个余弦函数y=cos(x)的图像并计算相位值，得到下面的图：

可以发现此时的相位也是周期性变化的，但实际上的相位值应该是单调递增的，因此有些场景是需要知道真实的相位值的，所以需要对相位值进行展开或者是解卷绕(unwrap)的操作。

目前有些关于解卷绕的文章要钱，有些讲得不清楚，所以分享给大家~

1 理论部分

不难发现，要实现相位值进行展开或者是解卷绕，其实只需要将图1中相位值递减的部分翻折上去，递增的部分平移上去即可。具体怎么实现呢？

实现算法：将图1中的周期性相位函数（蓝色的线）称为：卷绕函数 $\psi(i)$ ；待展开的函数称为解卷绕函数 $\phi(i)$ （将蓝色线翻折得到的一条单调的线）

如果卷绕函数 $\psi(i)$ 和解卷绕函数 $\phi(i)$ 的梯度方向相同，则有： $\phi (i)=2\pi k+\psi (i)$ 。也就是前面说的，递增的部分，就平移上去，平移量一定是2π的整数倍。
如果卷绕函数 $\psi(i)$ 和解卷绕函数 $\phi(i)$ 的梯度方向不同，则有： $\phi (i)=2\pi k-\psi (i)$ 。也就是前面说的，递减的部分，就翻折上去。
其中，开始时k=1，且每经过一个波峰和波谷，k=k+1。
至于1和2的公式是怎么来的，可以自己画个图思考一下就很容易推出来啦

2 代码部分

卷绕函数图

首先第一张图的代码如下：x是自变量，y是函数值，p是相位值

font2 = {'size': 13}
x=np.arange(0,20,0.01)
y=np.cos(x)
p=np.arccos(y)
plt.subplot(2,1,1)
plt.plot(x,y)
plt.ylabel('y',font2)
plt.xlabel('x',font2)
plt.subplot(2,1,2)
plt.plot(x,p,'b')
plt.ylabel('phase',font2)
plt.xlabel('x',font2)
plt.show()

标记极值点

有了前面的理论，我们就可以写代码来实现了

首先是，我们需要知道k的值，因此需要把波峰和波谷都找出来（这里可以直接找极值点，如果是实验数据有波动的话，可能会有一些不是波峰波谷的极值点被标记出来，需要进行数据平滑和寻峰算法）。

def get_maxima(values: np.ndarray):
    """极大值点"""
    max_index = sg.argrelmax(values)[0]
    return max_index

def get_minima(values: np.ndarray):
    """极小值点"""
    min_index = sg.argrelmin(values)[0]  # 极小值的下标
    return min_index  # 返回极小值

max_index = get_maxima(p)
min_index = get_minima(p)

extremum=np.concatenate((max_index,min_index))
extremum=np.sort(extremum)
print(extremum)
plt.plot(x,p)
plt.scatter(x[extremum],p[extremum],s=100,c='red',marker='*')
plt.xlabel('time-index',font2)
plt.ylabel('phase',font2)
plt.title('normal phase value')
plt.grid()
plt.show()

结果如下所示：

相位展开

就是实现前面理论部分的算法，直接看代码吧~

dotnumber=len(extremum)
num=0
k=1
n=-1         #用于判断“+”和“-”
p1=np.array(p) #解卷绕的相位值
            
for i in range(0,len(p1)):
    if num<dotnumber-1:
        if extremum[num+1]>i>extremum[num]:
            p1[i]=2*k*np.pi+n*p1[i]
        elif extremum[num+1]==i:
            p1[i]=2*k*np.pi+n*p1[i]
            num=num+1
            if num%2==0:   #每经过一个波峰和波谷，k就加一
                k=k+1
            n=n*-1
    if i>extremum[dotnumber-1]:
        if num%2==0:
            p1[i]=2*k*np.pi-p1[i]
        else:
            p1[i]=2*k*np.pi+p1[i]

plt.plot(p,'b')
plt.plot(p1,'r')
plt.xlabel('time-index',font2)
plt.ylabel('phase',font2)
plt.title('normal phase value')
plt.legend([ 'wraping phase', 'unwraping phase'])
plt.grid()
plt.show()

得到结果如下：

可以发现，解卷绕后的相位就是一条单调的线了，成功~

需要指出的是：这里的代码仅仅是针对相位刚开始是递增的情况，如果刚开始是递减的情况，原理其实也是一样的，仅需要稍微修改一下解卷绕部分的代码就可以了。

完整代码如下：

import numpy as np
import matplotlib.pyplot as plt
import scipy.signal as sg

font2 = {'size': 13}
x=np.arange(0,20,0.01)
y=np.cos(x)
p=np.arccos(y)
plt.subplot(2,1,1)
plt.plot(x,y)
plt.ylabel('y',font2)
plt.xlabel('x',font2)
plt.subplot(2,1,2)
plt.plot(x,p,'b')
plt.ylabel('phase',font2)
plt.xlabel('x',font2)
plt.show()

def get_maxima(values: np.ndarray):
    """极大值点"""
    max_index = sg.argrelmax(values)[0]
    return max_index

def get_minima(values: np.ndarray):
    """极小值点"""
    min_index = sg.argrelmin(values)[0]  # 极小值的下标
    return min_index  # 返回极小值

max_index = get_maxima(p)
min_index = get_minima(p)

extremum=np.concatenate((max_index,min_index))
extremum=np.sort(extremum)
print(extremum)
plt.plot(x,p)
plt.scatter(x[extremum],p[extremum],s=100,c='red',marker='*')
plt.xlabel('time-index',font2)
plt.ylabel('phase',font2)
plt.title('normal phase value')
plt.grid()
plt.show()

dotnumber=len(extremum)
num=0
k=1
n=-1
p1=np.array(p)
            
for i in range(0,len(p1)):
    if num<dotnumber-1:
        if extremum[num+1]>i>extremum[num]:
            p1[i]=2*k*np.pi+n*p1[i]
        elif extremum[num+1]==i:
            p1[i]=2*k*np.pi+n*p1[i]
            num=num+1
            if num%2==0:
                k=k+1
            n=n*-1
    if i>extremum[dotnumber-1]:
        if num%2==0:
            p1[i]=2*k*np.pi-p1[i]
        else:
            p1[i]=2*k*np.pi+p1[i]

plt.plot(x,p,'b')
plt.plot(x,p1,'r')
plt.xlabel('time-index',font2)
plt.ylabel('phase',font2)
plt.title('normal phase value')
plt.legend([ 'wraping phase', 'unwraping phase'])
plt.grid()
plt.show()