DSPRelated.com
Code

Delay estimation revisited

Markus Nentwig June 9, 20125 comments Coded in Matlab

Estimates delay and scaling factor between two signals by iteration ( Code snippet: )

Introduction

Quite many DSP-related tasks involve an unknown time delay between signals that needs to be estimated, and maybe corrected. The code snippet presents a "quick-and-dirty" solution - it's not intended for implementation on a DSP, but for testbenches and the like, where robustness and simplicity are more important than efficiency.
The topic is well-researched, for example the references in [1] should be a good starting point.
A Matlab implementation that is stated to be maximum-likelyhood (which would be the "gold standard" for delay estimation) can be found in [2].

Figure 1 shows an example signal that has been delayed and scaled. The code snippet estimates delay and scaling factor, "un-delays" and scales back so that the output ideally matches the input.
The algorithm is simple but robust. It relies on FFT-based crosscorrelation, for example as illustrated in [3].
More sophisticated frequency domain methods such as [4] are more accurate and efficient, but the difficulty of resolving the +/- pi phase ambiguity may cause them to fail for some input signals, depending on power spectrum and group delay profile. The presented code snippet finds the time delay that maximizes crosscorrelation, regardless of whether "the" delay between the signals over frequency is well-defined or not.

Signals do not need to be coarse-aligned or even arrive in the right order (see "cyclic signals" section below).


Figure 1: delayed / scaled signals

Algorithm

The method isn't particularly imaginative - and that's an advertised feature, as little can go wrong:
Crosscorrelate, find the highest peak, make an educated guess where the maximum is hiding between the bins, time-shift and repeat.
The flowchart gives a general overview.


Figure 2: Algorithm

The actual iteration loop is somewhat more complex, as finite numerical accuracy gives bitwise identical values of the (shallow) correlation peak over some delay range. If this is observed for two points in the search window, the algorithm will continue to improve the third point, or exit if all three are identical. For a required accuracy up to 10-5 Tsample, this is usually not an issue.
The "educated guess" fits a square polynomial through the three points of the search window and finds the maximum.
After some iterations, numerical accuracy will cause nonsense guesses. If the predicted location of the maximum is outside the search window, it is discarded and linear interpolation is used instead. Otherwise, the algorithm will converge, if slowly.

The polynomial equation is derived here (input file for Maxima computer algebra system (open-source) ).
Methods are known for the interpolation of peaks in FFT data [5], [6] (note, the problem appears similar but is probably not exactly the same). Those could be used to speed up convergence, but I haven't investigated them, as the numerical accuracy of the input data appears to be the main limitation.

Cyclic signals

The signals are treated as cyclic. The first and last sample point have no special meaning, i.e. circshifting both input signals by the same amount will give the same result.
Delaying a signal with sharp transients may reveal ringing that is usually invisible, as long as one stays on the sampling grid:
Nyquist's vestigial sideband theorem [7] promises "no intersample interference", but it keeps quiet about the times between samples...

The delay is returned in the interval [-n/2..n/2], with a signal length of n.

Code snippet

The code snippet includes a demo function that sets up test signals with a known delay. It calls iterDelayEst(...), applies the estimated delay / scaling factor and compares. The result is shown in Figure 3.


Figure 3: Output

Conclusion

The function should return delay and scaling factor to maximize the crosscorrelation, and not be fussy about it.
It's a cleaned-up rewrite, please let me know (comment!) in case of problems.

References

[1] Y. Zhang; W. H. Abdulla: A Comparative Study of Time-Delay Estimation Techniques Using Microphone Arrays

[2] Moddemeijer, R., An information theoretical delay estimator

[3] J. del Peral-Rosado et al: Preliminary Analysis of the Positioning Capabilities of the Positioning Reference Signal of 3GPP LTE

[4] M. Nentwig: Delay estimation by FFT

[5] R. Lyons: Accurate Measurement of a Sinusoid's Peak Amplitude Based on FFT Data

[6] Tutorial: Interpolating the peak location of a FFT

[7] Lecture notes: Pulse shaping and equalization

% ****************************************************************
% find sub-sample delay and scaling factor between two cyclic signals 
% to maximize crosscorrelation
% Markus Nentwig, 120609_v1.1
% ****************************************************************
function iterDelayEstDemo();
    close all;
    n = 1024;    
    % n = 1024 * 256; disp('*** test: long signal enabled ***');

    % ****************************************************************
    % random signal
    % ****************************************************************
    fd = randn(1, n) + 1i * randn(1, n);
    
    % ****************************************************************
    % lowpass filter
    % ****************************************************************
    f = binFreq(n);
    fd(abs(f) > 0.045) = 0;
    s1 = real(ifft(fd)) * sqrt(n);
    
    % ****************************************************************
    % create delayed 2nd signal
    % ****************************************************************
    dTest_samples = 12.3456;
    cTest = 1.23456;
    % cTest = cTest + 1i; disp('*** test: complex coeff enabled ***'); 
    % cTest = -cTest; disp('*** test: negative coeff enabled ***'); 

    s2 = cTest * cs_delay(s1, 1, dTest_samples);
    % s2 = s2 + 0.5*randn(size(s2)); disp('*** test: noise enabled ***');
    
    % ****************************************************************
    % estimate delay
    % ****************************************************************
    [delay_samples, coeff] = iterDelayEst(s1, s2);

    % ****************************************************************
    % correct it
    % ****************************************************************
    s2a = cs_delay(s1, 1, delay_samples);
    s2b = s2a * coeff;

    figure(); hold on;
    h = plot(real(s1), 'k'); set(h, 'lineWidth', 3);
    h = plot(real(s2), 'b'); set(h, 'lineWidth', 3);
    h = plot(real(s2a), 'r'); set(h, 'lineWidth', 1);
    h = plot(real(s2b), 'm'); set(h, 'lineWidth', 1);
    xlim([1, numel(s1)]);
    xlabel('samples');
    legend('s1', 's2', 's2 un-delayed', 's2 un-delayed and scaled');
    title('test signals');
    
    format long;
    disp('nominal delay of s2 relative to s1')';
    dTest_samples
    disp('iterDelayEst() returned:');
    delay_samples
    disp('original scaling factor:');
    cTest
    disp('estimated scaling factor:');
    coeff
end

% ****************************************************************
% estimates delay and scaling factor 
% ****************************************************************
function [delay_samples, coeff] = iterDelayEst(s1, s2)
    
    s1 = s1(:) .'; % force row vectors
    s2 = s2(:) .';
    rflag = isreal(s1) && isreal(s2);
    n = numel(s1);
    halfN = floor(n/2);
    assert(numel(s2) == n, 'signals must have same length');

    % ****************************************************************
    % constants
    % ****************************************************************    
    % exit if uncertainty below threshold
    thr_samples = 1e-7;

    % exit after fixed number of iterations
    nIter = 25;

    % frequency domain representation of signals
    fd1 = fft(s1);
    fd2 = fft(s2);    

    % first round: No delay was applied
    tau = [];
    fd2Tau = fd2; % delayed s2 in freq. domain
    
    % frequency corresponding to each FFT bin -0.5..0.5
    f = binFreq(n);

    % uncertainty plot data
    e = [];

    % normalization factor
    nf = sqrt((fd1 * fd1') * (fd2 * fd2')) / n; % normalizes to 1
    
    % search window: 
    % known maximum and two surrounding points
    x1 = -1;
    x2 = -1;
    x3 = -1;
    y1 = -1;
    y2 = -1;
    y3 = -1;
    
    % ****************************************************************
    % iteration loop
    % ****************************************************************
    for count = 1:nIter
    
        % ****************************************************************
        % crosscorrelation with time-shifted signal
        % ****************************************************************
        xcorr = abs(ifft(fd2Tau .* conj(fd1)))/ nf;

        % ****************************************************************
        % detect peak
        % ****************************************************************
        if isempty(tau)

            % ****************************************************************
            % startup
            % initialize with three adjacent bins around peak
            % ****************************************************************
            ix = find(xcorr == max(xcorr));
            ix = ix(1); % use any, if multiple bitwise identical peaks
            
            % indices of three bins around peak
            ixLow = mod(ix-1-1, n) + 1; % one below
            ixMid = ix;
            ixHigh = mod(ix-1+1, n) + 1; % one above

            % delay corresponding to the three bins
            tauLow = mod(ixLow -1 + halfN, n) - halfN;
            tauMid = mod(ixMid -1 + halfN, n) - halfN;         
            tauHigh = mod(ixHigh -1 + halfN, n) - halfN; 

            % crosscorrelation value for the three bins
            xcLow = xcorr(ixLow);
            xcMid = xcorr(ixMid);
            xcHigh = xcorr(ixHigh);
            
            x1 = tauLow;
            x2 = tauMid;
            x3 = tauHigh;
            y1 = xcLow;
            y2 = xcMid; 
            y3 = xcHigh;
        else
            % ****************************************************************
            % only main peak at first bin is of interest
            % ****************************************************************
            tauMid = tau;
            xcMid = xcorr(1);

            if xcMid > y2
                % ****************************************************************
                % improve midpoint
                % ****************************************************************
                if tauMid > x2
                    % midpoint becomes lower point
                    x1 = x2;
                    y1 = y2;
                else
                    % midpoint becomes upper point
                    x3 = x2;
                    y3 = y2;
                end
                x2 = tauMid;
                y2 = xcMid;
            
            elseif tauMid < x2
                % ****************************************************************
                % improve low point
                % ****************************************************************
                assert(tauMid >= x1); % bitwise identical is OK
                assert(tauMid > x1 || xcMid > y1); % expect improvement
                x1 = tauMid;
                y1 = xcMid;
            elseif tauMid > x2 
                % ****************************************************************
                % improve high point
                % ****************************************************************
                assert(tauMid <= x3); % bitwise identical is OK                
                assert((tauMid < x3) || (xcMid > y3)); % expect improvement
                x3 = tauMid;
                y3 = xcMid;
            else
                assert(false, '?? evaluated for existing tau ??');
            end
        end

        % ****************************************************************
        % calculate uncertainty (window width)
        % ****************************************************************
        eIter = abs(x3 - x1);
        e = [e, eIter];
        if eIter < thr_samples
            % disp('threshold reached, exiting');
            break;
        end

        if y1 == y2 || y2 == y3
            % reached limit of numerical accuracy on one side
            usePoly = 0;
        else
            % ****************************************************************
            % fit 2nd order polynomial and find maximum
            % ****************************************************************
            num = (x2^2-x1^2)*y3+(x1^2-x3^2)*y2+(x3^2-x2^2)*y1;
            denom = (2*x2-2*x1)*y3+(2*x1-2*x3)*y2+(2*x3-2*x2)*y1;
            if denom ~= 0
                tau = num / denom;
                % is the point within [x1, x3]? 
                usePoly = ((tau > x1) && (tau < x3));
            else
                usePoly = 0;
            end            
        end
        if ~usePoly
            % revert to linear interpolation on the side with the
            % less-accurate outer sample 
            % Note: There is no guarantee that the side with the more accurate
            % outer sample is the right one, as the samples aren't 
            % placed on a regular grid!
            % Therefore, iterate to improve the "worse" side, which will
            % eventually become the "better side", and iteration converges.
            
            tauLow = (x1 + x2) / 2;
            tauHigh = (x2 + x3) / 2;
            if y1 < y3
                o = [tauLow, tauHigh];
            else
                o = [tauHigh, tauLow];                
            end
            % don't choose point that is identical to one that is already known
            tau = o(1);
            if tau == x1 || tau == x2 || tau == x3
                tau = o(2);
                if tau == x1 || tau == x2 || tau == x3
                    break;
                end
            end
        end

        % ****************************************************************
        % advance 2nd signal according to location of maximum
        % phase shift in frequency domain - delay in time domain
        % ****************************************************************
        fd2Tau = fd2 .* exp(2i * pi * f * tau);
    end % for

    % ****************************************************************
    % plot the uncertainty (window size) over the number of iterations
    % ****************************************************************
    if true
        figure(); semilogy(e, '+-'); grid on;
        xlabel('iteration');
        title('uncertainty in delay');
    end

    % ****************************************************************
    % the delay estimate is the final location of the delay that 
    % maximized crosscorrelation (center of window).
    % ****************************************************************
    delay_samples = x2;

    % ****************************************************************
    % Coefficient: Turn signal 1 into signal 2
    % ****************************************************************
    coeff = fd2Tau * fd1' ./ (fd1 * fd1');

    % ****************************************************************
    % chop roundoff error, if input signals are known to be 
    % real-valued.
    % ****************************************************************
    if rflag
        coeff = real(coeff);
    end
end

% ****************************************************************
% frequency corresponding to FFT bin
% ****************************************************************
function f = binFreq(n)
    f = (mod(((0:n-1)+floor(n/2)), n)-floor(n/2))/n;
end

% ****************************************************************
% delay by phase shift
% needed only for demo code
% ****************************************************************
function waveform = cs_delay(waveform, rate_Hz, delay_s)
    rflag = isreal(waveform);
    
    n = numel(waveform);
    cycLen_s = n / rate_Hz;
    nCyc = delay_s / cycLen_s();
    
    f = 0:(n - 1);
    f = f + floor(n / 2);
    f = mod(f, n);
    f = f - floor(n / 2);
    
    phase = -2 * pi * f * nCyc;
    rot = exp(1i*phase);
    
    waveform = ifft(fft(waveform) .* rot);
    if rflag
        waveform = real(waveform);
    end
end