【DSP】理解并用MATLAB手动实现unwrap()

unwrap函数详解

• 一般在计算一个系统相频特性时，就要用到反正切函数提取相位，计算机中反正切函数规定：在一、二象限中的角度为0 ~ π，三、四象限的角度为0 ~ -π
• 但实际得到的结果会发生相位跳变，跳变幅度为2π，这就叫做相位的卷绕
• unwrap的作用就是解卷绕，使相位在π处不发生跳变，从而反映出真实的相位变化

unwarp标准系统函数

• unwrap(pha, tol, dim)
  • pha为数列或矩阵
  • tol为判断的阈值
  • dim可以设置为操作对象
    • 若对每一列进行操作，那参数可填1或者省略
    • 若对每一行进行操作，那参数需填2

手动实现unwrap

拟一组数据

• 假设此时我们拥有一组数据：Data = [-59, -92, -126, -158, 166, 133, 96, 61, 26, -7, -18, -59, -92, -126, -158, 166, 133, 96, 61, 26, -7, -18];

• 

• 如上图所示，这些数据呈现一种类似周期性的循环，但这并不是真的周期函数，只是因为当他的数值跨越 ±π 时，就会根据其变化形式，条编导边界重新计算

• 这种状态我们一般称之为相位的卷绕或相位的折叠，用wrap表示，则相位饿解卷绕被称为unwrap

分析拟定数据

• 数据本身呈现递减趋势，从-59°开始不断递减
• 数据中在 -158° -> 166° 出现第一次跳变，跳变后继续下降，直到下一次跳变的出现

修订拟定数据 - 解卷绕

1. 恢复原本的度数 -158° -> 166°
   • -158°到边界点的距离 = |-158° - (-180°)| = 22°
   • 边界点到166°的距离 = |180° - 166°| = 14°
   • 故整体变化度数 = 22° + 14° = 36°
   • 原本度数 = -158° - 36° = -194°
2. 观察跳变数据
   • 166 - (-194°) = 360°，恰好为360°，从-180°跳到了180°，将函数值增加了360°
3. 得出本质
   • 根据函数值的变化，wrap的本质就是：当相位超出边界时，补偿360°并继续计算，而unwrap就是为了消除这种补偿
4. 后续数据
   • 后续数据也应该随着跳变点处恢复原本数据而进行改变，即在遇到下一个跳变点之前，所有数据均应该减去360°得到原始数值
   • 在遇到第二个跳变点时，由于进一步的跳变令补偿再次增加了360°，故第二次跳变之后的值到第三次跳变到来之前，所有数据应该在原本基础上修正360°*2 = 720°

整理公式

• 根据上述分析unwrap与wrap的关系可以表达为
  $$unwrap(x)=wrap(x)\pm 360\ast wrapNumber(x)$$
• 公式中的±分别对应数据递增或递减的情况，将其使用wrapCycle表示，公式可变为
  $$unwrap(x)=wrap(x)-wrapCycle\ast wrapNumber(x)$$
• 公式中每个对应的数据点均为其对应一个wrapCycle表示其跳变，公式可变为
  $$unwrap(x)=wrap(x)-wrapCycle(x)\ast wrapNumber(x)$$

MATLAB手动实现

%% unwrap
% unwrap(x) = wrap(x) - wrapCycle*wrapNum(x)

% data0
Data = [-59, -92, -126, -158, 166, 133, 96, 61, 26, -7, -18, -59, -92, -126, -158, 166, 133, 96, 61, 26, -7, -18];
subplot(3, 1, 1);
plot(1:length(Data), Data);
title('Ram simulated data')
grid on;

% use system function unwrap
unwrap_Data = unwrap(deg2rad(Data), [], 2);
unwrap_finallyData = rad2deg(unwrap_Data);
subplot(3, 1, 2);
plot(1:length(Data), unwrap_finallyData);
title('Use system function unwrap')
grid on;

% try to simulate unwrap
length_Data = length(Data);
finallyData = [];
CycleNumber = [];
wrapCycle = [];
% calculata CycleNum
for i = 1 : length_Data
    if i > 1
        if (Data(i) - Data(i-1)) > 180 
            CycleNumber(i)= CycleNumber(i-1) + 1;
        elseif (Data(i) - Data(i-1)) < -180
            CycleNumber(i) = CycleNumber(i-1) + 1;
        else
            CycleNumber(i) = CycleNumber(i-1);
        end
    else
        CycleNumber(i) = 0;
    end
end
% calculata wrapCycle
for i = 1 : length_Data
    if i > 1
        if (Data(i) - Data(i-1)) > 180 
            wrapCycle(i)= 360;
        elseif (Data(i) - Data(i-1)) < -180
            wrapCycle(i) = -360;
        else
            wrapCycle(i) = wrapCycle(i-1);
        end
    else
        wrapCycle(i) = 0;
    end        
end
% updata Data
for i = 1 : length_Data
    finallyData(i) = Data(i) - wrapCycle(i)*CycleNumber(i);
end
subplot(3, 1, 3);
plot(1:length(Data), finallyData);
title('try to simulate unwrap')
grid on;

MATLAB代码解释

1. % data0
   • 拟定一组含有跳变的数据值，并将其进行绘制，作为对比
2. % use system function unwrap
   • 使用系统标准unwrap函数进行解卷绕
   • 由于拟定数据表示角度，故在调用unwrap()函数需要对数据进行deg2rad角度到弧度的转换
   • 使用unwrap()解卷绕之后，为方便对比，对数据进行rad2deg弧度到角度的转换
   • 进行绘制作为对比
3. % try to simulate unwrap
   • 手动实现unwrap
4. % calculata CycleNum
   • 第一个数据跳变次数为0
   • 之后所有的数据是否跳变取决于与前一个数据的差值是否超过180°
   • 当跳变出现，说明针对前序数据发生一次跳变，即新跳变数据的跳变数比前一次跳变的跳变数多一
5. % calculata wrapCycle
   • 第一个数据的跳变大小为0
   • 之后所有的数据跳变大小取决于与前一个数据的差值是否超过180°
     • 若差值在[-180, 180]，跳变大小与前序数据跳变大小相等
     • 若差值 > 180， 跳变大小为360°进行解卷绕
     • 若差值 < -180， 跳变大小为-360°进行解卷绕
6. % updata Data
   • 更新手动实现的unwrap数据结果
   • 进行绘制作为对比

MATLAB结果对比

• 数据对比

• 图像对比