微分-解卷绕-积分算法提升matlab unwrap函数解卷绕范围

相位调制的信号，采用IQ数字解调，或者硬件上正交解调的相位，经反正切计算后，phase相位值落在-π到π范围内而不连续。因此，matlab的unwrap函数修正向量phase中的弧度相位角，当phase的连续元素之间的绝对跳跃大于或等于π弧度的默认跳跃公差时，将±2π的倍数相加。

这种算法，当真实的绝对跳跃在π范围内是准确的，大于π时将会出现误差。最近的一篇论文分析并解决了这个问题。论文：

Large Dynamic Range Optical Fiber Distributed Acoustic Sensing (DAS) With Differential-Unwrapping-Integral Algorithm | IEEE Journals & Magazine | IEEE Xplore

接下来我用仿真对论文中提出的算法进行了分析，matlab代码和figure结合分析：

一.图1.1表示仿真中假设需要解调的相位phase,是一个幅值A=10，频率为1KHz的正弦信号，采样率为50KHz，总的采样时间为1s，初始幅值，一阶微分为0。通过phase_wrap = wrapToPi(phase)函数，图1.2表示函数将phase卷绕到-π-π范围内。再通过phase_wrap_unwrap = unwrap(phase_wrap)函数，解卷绕，图1.3表示解卷绕后的phase_wrap_unwrap 和真实相位phase吻合，即解卷绕后的结果是正确的。

clear,
clc,
f_vibration = 1e3;
f_sampling = 50e3;
t_total = 1;
t_axis = 0:1/f_sampling:(t_total-1/f_sampling);
A = 10;
phase = A*cos(2*pi*f_vibration*t_axis);
phase = phase - phase(1);
figure(11),
set(figure(11),'unit','centimeters','position', [5 3 20 15]);
subplot(411),plot(t_axis,phase),xlim([0,0.01]),
xlabel('Time(s)'),ylabel('Amp (rad)'),title('phase A=10'),

phase_wrap = wrapToPi(phase);
subplot(412),plot(t_axis,phase_wrap),xlim([0,0.01]),
xlabel('Time(s)'),ylabel('Amp (rad)'),title('phase wrap'),

phase_wrap_unwrap = unwrap(phase_wrap);
subplot(413),plot(t_axis,phase,'--','LineWidth',1.2),xlim([0,0.01]),hold on
plot(t_axis,phase_wrap_unwrap ),xlim([0,0.01]),hold off
xlabel('Time(s)'),ylabel('Amp (rad)'),title('1、phase 2、phase wrap unwrap '),


phase_diff = [0,diff(phase)];
subplot(414),plot(t_axis,phase_diff),xlim([0,0.01]),hold on,
plot(t_axis,pi*ones(length(t_axis),1),'--'),xlim([0,0.01]),
plot(t_axis,-pi*ones(length(t_axis),1),'--'),xlim([0,0.01]),hold off,
xlabel('Time(s)'),ylabel('delta(Amp) (rad)'),title('phase diff'),

 二.然而该正弦信号幅值A=50时，解卷绕的相位信号phase_wrap_unwrap 和phase有差别，如图2.3所示。原因是有些时刻的相位变化大于π，相位变化可用一阶微分表示，phase_diff = [0,diff(phase)]，如图2.4，某些位置一阶微分绝对值大于π，phase_wrap_unwrap将会在这些地方积累整数倍2π的误差。如图1.4，一阶微分绝对值小于π，则不会出现误差。

f_vibration = 1e3;
f_sampling = 50e3;
t_total = 1;
t_axis = 0:1/f_sampling:(t_total-1/f_sampling);
A = 50;
phase = A*cos(2*pi*f_vibration*t_axis);
phase = phase - phase(1);
figure(12),
set(figure(12),'unit','centimeters','position', [5 3 20 15]);
subplot(411),plot(t_axis,phase),xlim([0,0.01]),
xlabel('Time(s)'),ylabel('Amp (rad)'),title('phase A=50'),

phase_wrap = wrapToPi(phase);
subplot(412),plot(t_axis,phase_wrap),xlim([0,0.01]),
xlabel('Time(s)'),ylabel('Amp (rad)'),title('phase wrap'),

phase_wrap_unwrap = unwrap(phase_wrap);
subplot(413),plot(t_axis,phase,'--','LineWidth',1.2),xlim([0,0.01]),hold on
plot(t_axis,phase_wrap_unwrap ),xlim([0,0.01]),hold off
xlabel('Time(s)'),ylabel('Amp (rad)'),title('1、phase 2、phase wrap unwrap '),


phase_diff = [0,diff(phase)];
subplot(414),plot(t_axis,phase_diff),xlim([0,0.01]),hold on,
plot(t_axis,pi*ones(length(t_axis),1),'--'),xlim([0,0.01]),
plot(t_axis,-pi*ones(length(t_axis),1),'--'),xlim([0,0.01]),hold off,
xlabel('Time(s)'),ylabel('delta(Amp) (rad)'),title('phase diff'),

 三、针对相位跳跃大于π无法解卷绕的问题，论文提出将一阶微分phase_diff作为新的原相位函数，此时的相位跳跃是二阶微分函数phase_diff2 = [0,diff(phase_diff)]，如图3.1，phase_diff2绝对值不大于π。满足此条件的情况下，卷绕函数phase_wrap的微分phase_wrap_diff 等于微分的卷绕phase_diff，如图3.2。则phase_wrap_diff的解卷绕phase_wrap_diff_unwrap，等于原函数的微分phase_diff，如图3.3。将函数phase_wrap_diff_unwrap积分后得到的结果phase_wrap_diff_unwrap_int 等于原函数phase，如图3.4。

上述这段话有些绕口，总结起来就是一句绕口令：

如果相位函数的二阶微分绝对值不大于π，则该函数的卷绕的一阶微分，解卷绕后积分得到的结果，等于该函数。

f_vibration = 1e3;
f_sampling = 50e3;
t_total = 1;
t_axis = 0:1/f_sampling:(t_total-1/f_sampling);
A = 50;
phase = A*cos(2*pi*f_vibration*t_axis);
phase = phase - phase(1);
phase_diff = [0,diff(phase)];
phase_diff2 = [0,diff(phase_diff)];
figure(13),
set(figure(13),'unit','centimeters','position', [5 3 20 15]);
subplot(411),plot(t_axis,phase_diff2),xlim([0,0.01]),hold on,
plot(t_axis,pi*ones(length(t_axis),1),'--'),xlim([0,0.01]),
plot(t_axis,-pi*ones(length(t_axis),1),'--'),xlim([0,0.01]),hold off,
xlabel('Time(s)'),ylabel('delta(Amp)2 (rad)'),title('phase diff2 A=50'),

phase_diff_wrap = wrapToPi(phase_diff);
subplot(412),plot(t_axis,phase_diff_wrap,'--','LineWidth',1.2),xlim([0,0.01]),hold on,
phase_wrap = wrapToPi(phase);
phase_wrap_diff = wrapToPi([0,diff(phase_wrap)]);
plot(t_axis,phase_wrap_diff),xlim([0,0.01]),hold off,
xlabel('Time(s)'),ylabel('delta(Amp) (rad)'),title('1、phase diff wrap 2、phase wrap diff '),

phase_wrap_diff_unwrap =unwrap(phase_wrap_diff);
subplot(413),plot(t_axis,phase_diff,'--','LineWidth',1.2),xlim([0,0.01]),hold on,
plot(t_axis,phase_wrap_diff_unwrap),xlim([0,0.01]),hold off,
xlabel('Time(s)'),ylabel('delta(Amp) (rad)'),title('1、phase diff 2、phase wrap diff unwrap'),

phase_wrap_diff_unwrap_int = zeros(length(t_axis),1);
for i = 1:(length(t_axis)-1)
    phase_wrap_diff_unwrap_int(i+1) = phase_wrap_diff_unwrap_int(i) + phase_wrap_diff_unwrap(i+1);
end
subplot(414),plot(t_axis,phase,'--','LineWidth',1.2),xlim([0,0.01]),hold on,
plot(t_axis,phase_wrap_diff_unwrap_int),xlim([0,0.01]),hold off,title('1、phase 2、phase wrap diff unwrap int'),
xlabel('Time(s)'),ylabel('Amp (rad)')

 四、将正弦信号幅值A进一步增加到250，此时相位函数的二阶微分绝对值大于π，一阶微分也不能恢复原函数，需要求助于二阶微分。

f_vibration = 1e3;
f_sampling = 50e3;
t_total = 1;
t_axis = 0:1/f_sampling:(t_total-1/f_sampling);
A = 250;
phase = A*cos(2*pi*f_vibration*t_axis);
phase = phase - phase(1);
phase_diff = [0,diff(phase)];
phase_diff2 = [0,diff(phase_diff)];
figure(14),
set(figure(14),'unit','centimeters','position', [5 3 20 15]);
subplot(411),plot(t_axis,phase_diff2),xlim([0,0.01]),hold on,
plot(t_axis,pi*ones(length(t_axis),1),'--'),xlim([0,0.01]),
plot(t_axis,-pi*ones(length(t_axis),1),'--'),xlim([0,0.01]),hold off,
xlabel('Time(s)'),ylabel('delta(Amp)2 (rad)'),title('phase diff2 A=250'),

phase_diff_wrap = wrapToPi(phase_diff);
subplot(412),plot(t_axis,phase_diff_wrap,'--','LineWidth',1.2),xlim([0,0.01]),hold on,
phase_wrap = wrapToPi(phase);
phase_wrap_diff = wrapToPi([0,diff(phase_wrap)]);
plot(t_axis,phase_wrap_diff),xlim([0,0.01]),hold off,
xlabel('Time(s)'),ylabel('delta(Amp) (rad)'),title('1、phase diff wrap 2、phase wrap diff '),

phase_wrap_diff_unwrap =unwrap(phase_wrap_diff);
subplot(413),plot(t_axis,phase_diff,'--','LineWidth',1.2),xlim([0,0.01]),hold on,
plot(t_axis,phase_wrap_diff_unwrap),xlim([0,0.01]),hold off,
xlabel('Time(s)'),ylabel('delta(Amp) (rad)'),title('1、phase diff 2、phase wrap diff unwrap'),

phase_wrap_diff_unwrap_int = zeros(length(t_axis),1);
for i = 1:(length(t_axis)-1)
    phase_wrap_diff_unwrap_int(i+1) = phase_wrap_diff_unwrap_int(i) + phase_wrap_diff_unwrap(i+1);
end
subplot(414),plot(t_axis,phase,'--','LineWidth',1.2),xlim([0,0.01]),hold on,
plot(t_axis,phase_wrap_diff_unwrap_int),xlim([0,0.01]),hold off,title('1、phase 2、phase wrap diff unwrap int'),
xlabel('Time(s)'),ylabel('Amp (rad)'),

五、如果相位函数的三阶微分绝对值不大于π，则该函数的卷绕的二阶微分，解卷绕后二阶积分得到的结果，等于该函数。三阶，四阶以及N阶以此类推。

f_vibration = 1e3;
f_sampling = 50e3;
t_total = 1;
t_axis = 0:1/f_sampling:(t_total-1/f_sampling);
A = 250;
phase = A*cos(2*pi*f_vibration*t_axis);
phase = phase - phase(1);
phase_diff = [0,diff(phase)];
phase_diff2 = [0,diff(phase_diff)];
phase_diff3 = [0,diff(phase_diff2)];
figure(15),
set(figure(15),'unit','centimeters','position', [5 3 20 18]);
subplot(511),plot(t_axis,phase_diff3),xlim([0,0.01]),hold on,
plot(t_axis,pi*ones(length(t_axis),1),'--'),xlim([0,0.01]),
plot(t_axis,-pi*ones(length(t_axis),1),'--'),xlim([0,0.01]),hold off,
xlabel('Time(s)'),ylabel('delta(Amp)3 (rad)'),title('phase diff3 A=250'),

phase_diff2_wrap = wrapToPi(phase_diff2);
subplot(512),plot(t_axis,phase_diff2_wrap,'--','LineWidth',1.2),xlim([0,0.01]),hold on,
phase_wrap = wrapToPi(phase);
phase_wrap_diff = wrapToPi([0,diff(phase_wrap)]);
phase_wrap_diff2 = wrapToPi([0,diff(phase_wrap_diff)]);
plot(t_axis,phase_wrap_diff2),xlim([0,0.01]),hold off,
xlabel('Time(s)'),ylabel('delta(Amp)2 (rad)'),title('1、phase diff2 wrap 2、phase wrap diff2 '),

phase_wrap_diff2_unwrap =unwrap(phase_wrap_diff2);
subplot(513),plot(t_axis,phase_diff2,'--','LineWidth',1.2),xlim([0,0.01]),hold on,
plot(t_axis,phase_wrap_diff2_unwrap),xlim([0,0.01]),hold off,
xlabel('Time(s)'),ylabel('delta(Amp)2 (rad)'),title('1、phase diff2 2、phase wrap diff2 unwrap'),

phase_wrap_diff2_unwrap_int = zeros(length(t_axis),1);
for i = 1:(length(t_axis)-1)
    phase_wrap_diff2_unwrap_int(i+1) = phase_wrap_diff2_unwrap_int(i) + phase_wrap_diff2_unwrap(i+1);
end
subplot(514),plot(t_axis,phase_diff,'--','LineWidth',1.2),xlim([0,0.01]),hold on,
plot(t_axis,phase_wrap_diff2_unwrap_int),xlim([0,0.01]),hold off,title('1、phase diff  2、phase wrap diff2 unwrap int'),
xlabel('Time(s)'),ylabel('delta(Amp) (rad)'),

phase_wrap_diff2_unwrap_int2 = zeros(length(t_axis),1);
for i = 1:(length(t_axis)-1)
    phase_wrap_diff2_unwrap_int2(i+1) = phase_wrap_diff2_unwrap_int2(i) + phase_wrap_diff2_unwrap_int(i+1);
end
subplot(515),plot(t_axis,phase,'--','LineWidth',1.2),xlim([0,0.01]),hold on,
plot(t_axis,phase_wrap_diff2_unwrap_int2),xlim([0,0.01]),hold off,title('1、phase 2、phase wrap diff2 unwrap int2'),
xlabel('Time(s)'),ylabel('Amp (rad)'),

  六、将正弦信号幅值A进一步增加到500，如图6.1，此时相位函数的三阶微分，初略看绝对值小于π，但是起始位置绝对值大于π，最后结果就不能还原。

论文提出的算法涉及微分，积分操作，这里的不连续积分需要得到函数的初值，而我们事先不知道初值，用0代替的t=0之前的相位值。因此如图4.1所示，在t=0时刻，原相位函数在2阶微分上是不连续的，因此3阶微分起始位置绝对值大于π。随时间增加，误差逐渐累计，这个误差有可能可以用多项式拟合的方法消除。

f_vibration = 1e3;
f_sampling = 50e3;
t_total = 1;
t_axis = 0:1/f_sampling:(t_total-1/f_sampling);
A = 500;
phase = A*cos(2*pi*f_vibration*t_axis);
phase = phase - phase(1);
phase_diff = [0,diff(phase)];
phase_diff2 = [0,diff(phase_diff)];
phase_diff3 = [0,diff(phase_diff2)];
figure(16),
set(figure(16),'unit','centimeters','position', [5 3 20 18]);
subplot(511),plot(t_axis,phase_diff3),xlim([0,0.01]),hold on,
plot(t_axis,pi*ones(length(t_axis),1),'--'),xlim([0,0.01]),
plot(t_axis,-pi*ones(length(t_axis),1),'--'),xlim([0,0.01]),hold off,
xlabel('Time(s)'),ylabel('delta(Amp)2 (rad)'),title('phase diff3 A=500'),

phase_diff2_wrap = wrapToPi(phase_diff2);
subplot(512),plot(t_axis,phase_diff2_wrap,'--','LineWidth',1.2),xlim([0,0.01]),hold on,
phase_wrap = wrapToPi(phase);
phase_wrap_diff = wrapToPi([0,diff(phase_wrap)]);
phase_wrap_diff2 = wrapToPi([0,diff(phase_wrap_diff)]);
plot(t_axis,phase_wrap_diff2),xlim([0,0.01]),hold off,
xlabel('Time(s)'),ylabel('delta(Amp)2 (rad)'),title('1、phase diff2 wrap 2、phase wrap diff2 '),

phase_wrap_diff2_unwrap =unwrap(phase_wrap_diff2);
subplot(513),plot(t_axis,phase_diff2,'--','LineWidth',1.2),xlim([0,0.01]),hold on,
plot(t_axis,phase_wrap_diff2_unwrap),xlim([0,0.01]),hold off,
xlabel('Time(s)'),ylabel('delta(Amp)2 (rad)'),title('1、phase diff2 2、phase wrap diff2 unwrap'),

phase_wrap_diff2_unwrap_int = zeros(length(t_axis),1);
for i = 1:(length(t_axis)-1)
    phase_wrap_diff2_unwrap_int(i+1) = phase_wrap_diff2_unwrap_int(i) + phase_wrap_diff2_unwrap(i+1);
end
subplot(514),plot(t_axis,phase_diff,'--','LineWidth',1.2),xlim([0,0.01]),hold on,
plot(t_axis,phase_wrap_diff2_unwrap_int),xlim([0,0.01]),hold off,title('1、phase diff  2、phase wrap diff2 unwrap int'),
xlabel('Time(s)'),ylabel('delta(Amp) (rad)'),

phase_wrap_diff2_unwrap_int2 = zeros(length(t_axis),1);
for i = 1:(length(t_axis)-1)
    phase_wrap_diff2_unwrap_int2(i+1) = phase_wrap_diff2_unwrap_int2(i) + phase_wrap_diff2_unwrap_int(i+1);
end
subplot(515),plot(t_axis,phase,'--','LineWidth',1.2),xlim([0,0.01]),hold on,
plot(t_axis,phase_wrap_diff2_unwrap_int2),xlim([0,0.01]),hold off,title('1、phase 2、phase wrap diff2 unwrap int2'),
xlabel('Time(s)'),ylabel('Amp (rad)'),

 七、如果起始的几个点，相位的N阶微分均小于π，或者已知起始几个点的真实相位值，可以解决微分积分引入的初始相位误差。

f_vibration = 1e3;
f_sampling = 50e3;
t_total = 1;
t_axis = 0:1/f_sampling:(t_total-1/f_sampling);
A = 500;
phase = A*cos(2*pi*f_vibration*t_axis);
phase = phase - phase(1);
phase_diff = diff(phase);
phase_diff2 = diff(phase_diff);
phase_diff3 = diff(phase_diff2);
figure(17),
set(figure(17),'unit','centimeters','position', [5 3 20 18]);
subplot(511),plot(t_axis(4:end),phase_diff3),xlim([0,0.01]),hold on,
plot(t_axis,pi*ones(length(t_axis),1),'--'),xlim([0,0.01]),
plot(t_axis,-pi*ones(length(t_axis),1),'--'),xlim([0,0.01]),hold off,
xlabel('Time(s)'),ylabel('delta(Amp)3 (rad)'),title('phase diff3 A=500'),

phase_diff2_wrap = wrapToPi(phase_diff2);
subplot(512),plot(t_axis(3:end),phase_diff2_wrap,'--','LineWidth',1.2),xlim([0,0.01]),hold on,
phase_wrap = wrapToPi(phase);
phase_wrap_diff = wrapToPi(diff(phase_wrap));
phase_wrap_diff2 = wrapToPi(diff(phase_wrap_diff));
plot(t_axis(3:end),phase_wrap_diff2),xlim([0,0.01]),hold off,
xlabel('Time(s)'),ylabel('delta(Amp)2 (rad)'),title('1、phase diff2 wrap 2、phase wrap diff2 '),

phase_wrap_diff2_unwrap =unwrap(phase_wrap_diff2);
phase_wrap_diff2_unwrap = phase_wrap_diff2_unwrap - phase_wrap_diff2_unwrap(1) +  phase_diff2(1);

subplot(513),plot(t_axis(3:end),phase_diff2,'--','LineWidth',1.2),xlim([0,0.01]),hold on,
plot(t_axis(3:end),phase_wrap_diff2_unwrap),xlim([0,0.01]),hold off,
xlabel('Time(s)'),ylabel('delta(Amp)2 (rad)'),title('1、phase diff2 2、phase wrap diff2 unwrap'),

phase_wrap_diff2_unwrap_int = zeros(length(t_axis)-1,1);
phase_wrap_diff2_unwrap_int(1) = phase_diff(1);
for i = 1:(length(t_axis)-2)
    phase_wrap_diff2_unwrap_int(i+1) = phase_wrap_diff2_unwrap_int(i) + phase_wrap_diff2_unwrap(i);
end

subplot(514),plot(t_axis(2:end),phase_diff,'--','LineWidth',1.2),xlim([0,0.01]),hold on,
plot(t_axis(2:end),phase_wrap_diff2_unwrap_int),xlim([0,0.01]),hold off,title('1、phase diff  2、phase wrap diff2 unwrap int'),
xlabel('Time(s)'),ylabel('delta(Amp) (rad)'),
phase_wrap_diff2_unwrap_int2 = zeros(length(t_axis),1);
phase_wrap_diff2_unwrap_int2(1) = phase(1);
for i = 1:(length(t_axis)-1)
    phase_wrap_diff2_unwrap_int2(i+1) = phase_wrap_diff2_unwrap_int2(i) + phase_wrap_diff2_unwrap_int(i);
end
subplot(515),plot(t_axis,phase,'--','LineWidth',1.2),xlim([0,0.01]),hold on,
plot(t_axis,phase_wrap_diff2_unwrap_int2),xlim([0,0.01]),hold off,title('1、phase 2、phase wrap diff2 unwrap int2'),
xlabel('Time(s)'),ylabel('Amp (rad)'),

八、以此类推，下面这段matlab程序实现了三阶微分积分的反卷积，正弦相位幅值A=3000。

f_vibration = 1e3;
f_sampling = 50e3;
t_total = 1;
t_axis = 0:1/f_sampling:(t_total-1/f_sampling);
A = 3000;
phase = A*cos(2*pi*f_vibration*t_axis);
phase = phase - phase(1);
phase_diff = diff(phase);
phase_diff2 = diff(phase_diff);
phase_diff3 = diff(phase_diff2);
phase_diff4 = diff(phase_diff3);


phase_wrap = wrapToPi(phase);
phase_wrap_diff = wrapToPi(diff(phase_wrap));
phase_wrap_diff2 = wrapToPi(diff(phase_wrap_diff));
phase_wrap_diff3 = wrapToPi(diff(phase_wrap_diff2));

phase_wrap_diff3_unwrap =unwrap(phase_wrap_diff3);
phase_wrap_diff3_unwrap = phase_wrap_diff3_unwrap - phase_wrap_diff3_unwrap(1) +  phase_diff3(1);

phase_wrap_diff3_unwrap_int = zeros(length(t_axis)-2,1);
phase_wrap_diff3_unwrap_int(1) = phase_diff2(1);
for i = 1:(length(t_axis)-3)
    phase_wrap_diff3_unwrap_int(i+1) = phase_wrap_diff3_unwrap_int(i) + phase_wrap_diff3_unwrap(i);
end

phase_wrap_diff3_unwrap_int2 = zeros(length(t_axis)-1,1);
phase_wrap_diff3_unwrap_int2(1) = phase_diff(1);
for i = 1:(length(t_axis)-2)
    phase_wrap_diff3_unwrap_int2(i+1) = phase_wrap_diff3_unwrap_int2(i) + phase_wrap_diff3_unwrap_int(i);
end

phase_wrap_diff3_unwrap_int3 = zeros(length(t_axis),1);
phase_wrap_diff3_unwrap_int3(1) = phase(1);
for i = 1:(length(t_axis)-1)
    phase_wrap_diff3_unwrap_int3(i+1) = phase_wrap_diff3_unwrap_int3(i) + phase_wrap_diff3_unwrap_int2(i);
end

figure(18),
set(figure(18),'unit','centimeters','position', [5 3 20 6]);
plot(t_axis,phase,'--','LineWidth',1.2),xlim([0,0.01]),hold on,
plot(t_axis,phase_wrap_diff3_unwrap_int3),xlim([0,0.01]),hold off,title('1、phase 2、phase wrap diff3 unwrap int3 A=3000'),
xlabel('Time(s)'),ylabel('Amp (rad)'),

 九、总结：微分-解卷绕-积分算法提升了matlab unwrap函数解卷绕范围。但这种算法对初始相位要求比较苛刻。论文中似乎并没有提到这一问题。

论文中还定量讨论了可恢复相位幅值与阶次之间的关系，噪声与阶次的关系，阶次的选择问题，多频相位信号的恢复要求，等等。