unwrap_phase.m

来自「实现地震勘探中」· M 代码 · 共 41 行

function phase=unwrap_phase(phase)
% Unwrap phase as described in 
% Yanghua Wang, 2000, A robust seismic phase unwrapping method;
% Journal of Seismic Exploration, 93-104
% 
% Written by: E. R.: March 4, 2004
% Last updated:
%
%          phase=unwrap_phase(phase)   
% INPUT
% phase    wrapped phase (radians); a linear trend should have been removed.
% OUTPUT
% phase    unwrapped phase (radians)

[phase,ndims]=shiftdim(phase);  % Make sure that first dimension is not singleton

dphase=diff(phase);
idx=find(abs(dphase) > pi);     % Find phase changes > pi
if ~isempty(idx)
   dphase(idx)=mod(dphase(idx)+pi,2*pi)-pi;
end

bool=abs(dphase) > pi/2;

bool(end)=false;


if any(bool)
   bool(2:end-1)=bool(2:end-1) | (bool(1:end-2) & bool(3:end));
   bool=[0;bool];
   idx1=find(diff(bool) ~= 0);
   for ii=1:2:length(idx1)
      temp=mod(phase(idx1(ii+1)+1)-phase(idx1(ii)+1)+pi,2*pi)-pi;
%      dphase(idx1(ii))=mod(phase(idx1(ii+1)+1)-phase(idx1(ii)+1)+pi,2*pi)-pi;
      dphase(idx1(ii):idx1(ii+1))=temp/(idx1(ii+1)-idx1(ii));
   end
end
dphase=[phase(1,:);dphase];

phase=shiftdim(cumsum(dphase),-ndims);  % Undo dimension change, if there was one