利用两路信号的相位差对时域延迟进行校正

clc;clear all;close all;

%--------------------------------------------------------------------------
%设置参数
FS      = 9000;
TS      = 1/FS;
N       = 16384;
t       = (0:N-1)*TS;

CHAN_NUM      = 64;             %数据的路数
UNWRAP_MATLAB = 1;              %是否使用matlab内置函数unwrap
PI_N          = pi/16;
CLK_MHZ       = FS/CHAN_NUM;
delta_f12     = PI_N*FS/(2*pi)/CHAN_NUM;
delta_f23     = CLK_MHZ;

for mm=1:CHAN_NUM*4
    DELAY = mm;
    f1    = 500;
    f2    = f1 + delta_f12;
    f3    = f2 + delta_f23;
    
    sig1  = cos(2*pi*f1*t);
    sig2  = cos(2*pi*f2*t);
    sig3  = cos(2*pi*f3*t);
    
    %----------------------------------------------------------------------
    % 生成AD1和AD2的采样信号
    sig_AD1 = [sig1;sig2;sig3];
    sig_AD2 = [zeros(1,DELAY) sig1(1,1:end-DELAY);zeros(1,DELAY) sig2(1,1:end-DELAY);zeros(1,DELAY) sig3(1,1:end-DELAY)];
    
    %----------------------------------------------------------------------
    % 计算相位
    sig_AD1_phase    = zeros(3,1);
    sig_AD2_phase    = zeros(3,1);
    delta_AD12_phase = zeros(3,1);
    for m=1:3
        y_fft                 = fft(sig_AD1(m,:));
        [max_value,max_index] = max(abs(y_fft(1,1:N/2)));
        sig_AD1_phase(m,1)    = angle(y_fft(1,max_index));
        
        y_fft                 = fft(sig_AD2(m,:));
        [max_value,max_index] = max(abs(y_fft(1,1:N/2)));
        sig_AD2_phase(m,1)    = angle(y_fft(1,max_index));
    end
    
    %----------------------------------------------------------------------
    % 手动 对相位差unwrap
    if UNWRAP_MATLAB==0
        for m=1:3
           if(sig_AD2_phase(m,1)-sig_AD1_phase(m,1) >=pi)
               delta_AD12_phase(m,1) = sig_AD2_phase(m,1)-sig_AD1_phase(m,1)-2*pi;
           elseif(sig_AD2_phase(m,1)-sig_AD1_phase(m,1) <=-pi)
               delta_AD12_phase(m,1) = sig_AD2_phase(m,1)-sig_AD1_phase(m,1)+2*pi;
           else
               delta_AD12_phase(m,1) = sig_AD2_phase(m,1)-sig_AD1_phase(m,1);
           end
        end
        
        delta = zeros(2,1);
        for m=1:2
           if(delta_AD12_phase(m+1,1)-delta_AD12_phase(m,1) >=pi)
               delta_phi(m,1) = delta_AD12_phase(m+1,1)-delta_AD12_phase(m,1)-2*pi;
           elseif(delta_AD12_phase(m+1,1)-delta_AD12_phase(m,1) <=-pi)
               delta_phi(m,1) = delta_AD12_phase(m+1,1)-delta_AD12_phase(m,1)+2*pi;
           else
               delta_phi(m,1) = delta_AD12_phase(m+1,1)-delta_AD12_phase(m,1);
           end
        end
    %使用内置函数
    else
        delta_AD12_phase = unwrap(sig_AD2_phase-sig_AD1_phase);
        for m=1:2
            delta_phi(m,1) = delta_AD12_phase(m+1,1)-delta_AD12_phase(m,1);
        end
    end
    
    %----------------------------------------------------------------------
    delay_get = round(delta_phi(1,1)/PI_N*CHAN_NUM);
    delay_log(mm,1) = delta_phi(1,1)/PI_N*CHAN_NUM;
end