Signal fitting with subsample resolution

Update: The most recent version can be found here and a demonstration here.
The downloadable code is a "real-life" implementation that is a good choice for "real-life problems", but has grown to almost 700 lines. The code below implements the main functionality only, but is shorter.

Determines sub-sample delay and scaling factor between two signals and resamples / scales one signal to match the other.

For example in measurement automation, it is rather common to apply a test signal to a device-under-test and record the output. A delay is caused both from the latency of measurement instruments and the group delay of the device-under-test itself. In some applications, for example when analyzing the distortion properties of power amplifiers, it is necessary to align the signals within a small fraction of a sample length to avoid systematic errors.

The function fitSignalFFT() determines the delay with sub-sample resolution. It provides a delayed-and-scaled version of one signal that matches the other one.

This code snippet provides an extended version of the implementation in my "delay estimation by FFT" blog entry.

Resampling a signal with a fractional (subsample) delay assumes reconstruction and resampling of an equivalent continuous-time waveform. There are a few "gotchas" such as "ringing" related to the underlying assumptions (cyclic signal, ideal lowpass as reconstruction filter). Some related background can be found in this code snippet (resampling on a regular grid via FFT) and here (resampling at arbitrary time instants by direct evaluation of the Fourier series).

For non-periodic signals, use adequate zero-padding.

The picture shows an example run:
Two sine waves with different phase and amplitude are used as inputs. The function call returns an offset in samples and a scaling factor that maps the reference signal least-squares-optimally to the other signal.

Further, a delayed and scaled version of the reference signal is returned (green trace, overlapping the blue trace).

Below the example invocation.

ph = (0:15) * 2 * pi / 16;
ref = sin(ph);
sig = 1.23456 * sin(ph + 0.98765)
[coeff, shiftedRef, delta] = fitSignal_FFT(sig, ref);
figure(); hold on;
plot(ref, 'k+-');
h = plot(sig, 'b+-'); set(h, 'lineWidth', 5);
plot(shiftedRef, 'g+-')
legend('reference signal', 'signal', 'reference shifted / scaled to signal');
title('fitSignal\_FFT demo');

Copy the following code snippet into a file fitSignal_FFT.m.

% *******************************************************
% delay-matching between two signals (complex/real-valued)
% M. Nentwig
%
% * matches the continuous-time equivalent waveforms
%   of the signal vectors (reconstruction at Nyquist limit =>
%   ideal lowpass filter)
% * Signals are considered cyclic. Use arbitrary-length 
%   zero-padding to turn a one-shot signal into a cyclic one.
%
% * output:
%   => coeff: complex scaling factor that scales 'ref' into 'signal'
%   => delay 'deltaN' in units of samples (subsample resolution)
%      apply both to minimize the least-square residual      
%   => 'shiftedRef': a shifted and scaled version of 'ref' that 
%      matches 'signal' 
%   => (signal - shiftedRef) gives the residual (vector error)
%
% *******************************************************
function [coeff, shiftedRef, deltaN] = fitSignal_FFT(signal, ref)
    n=length(signal);
    % xyz_FD: Frequency Domain
    % xyz_TD: Time Domain
    % all references to 'time' and 'frequency' are for illustration only

    forceReal = isreal(signal) && isreal(ref);
    
    % *******************************************************
    % Calculate the frequency that corresponds to each FFT bin
    % [-0.5..0.5[
    % *******************************************************
    binFreq=(mod(((0:n-1)+floor(n/2)), n)-floor(n/2))/n;

    % *******************************************************
    % Delay calculation starts:
    % Convert to frequency domain...
    % *******************************************************
    sig_FD = fft(signal);
    ref_FD = fft(ref, n);

    % *******************************************************
    % ... calculate crosscorrelation between 
    % signal and reference...
    % *******************************************************
    u=sig_FD .* conj(ref_FD);
    if mod(n, 2) == 0
        % for an even sized FFT the center bin represents a signal
        % [-1 1 -1 1 ...] (subject to interpretation). It cannot be delayed. 
        % The frequency component is therefore excluded from the calculation.
        u(length(u)/2+1)=0;
    end
    Xcor=abs(ifft(u));

    %  figure(); plot(abs(Xcor));
    
    % *******************************************************
    % Each bin in Xcor corresponds to a given delay in samples.
    % The bin with the highest absolute value corresponds to
    % the delay where maximum correlation occurs.
    % *******************************************************
    integerDelay = find(Xcor==max(Xcor));
    
    % (1): in case there are several bitwise identical peaks, use the first one
    % Minus one: Delay 0 appears in bin 1
    integerDelay=integerDelay(1)-1;

    % Fourier transform of a pulse shifted by one sample
    rotN = exp(2i*pi*integerDelay .* binFreq);

    uDelayPhase = -2*pi*binFreq;
    
    % *******************************************************
    % Since the signal was multiplied with the conjugate of the
    % reference, the phase is rotated back to 0 degrees in case
    % of no delay. Delay appears as linear increase in phase, but
    % it has discontinuities.
    % Use the known phase (with +/- 1/2 sample accuracy) to 
    % rotate back the phase. This removes the discontinuities.
    % *******************************************************
    %  figure(); plot(angle(u)); title('phase before rotation');
    u=u .* rotN;
    
    % figure(); plot(angle(u)); title('phase after rotation');
    
    % *******************************************************
    % Obtain the delay using linear least mean squares fit
    % The phase is weighted according to the amplitude.
    % This suppresses the error caused by frequencies with
    % little power, that may have radically different phase.
    % *******************************************************
    weight = abs(u); 
    constRotPhase = 1 .* weight;
    uDelayPhase = uDelayPhase .* weight;
    ang = angle(u) .* weight;
    r = [constRotPhase; uDelayPhase] .' \ ang.'; %linear mean square
    
    %rotPhase=r(1); % constant phase rotation, not used.
    % the same will be obtained via the phase of 'coeff' further down
    fractionalDelay=r(2);
    
    % *******************************************************
    % Finally, the total delay is the sum of integer part and
    % fractional part.
    % *******************************************************
    deltaN = integerDelay + fractionalDelay;

    % *******************************************************
    % provide shifted and scaled 'ref' signal
    % *******************************************************
    % this is effectively time-convolution with a unit pulse shifted by deltaN
    rotN = exp(-2i*pi*deltaN .* binFreq);
    ref_FD = ref_FD .* rotN;
    shiftedRef = ifft(ref_FD);
    
    % *******************************************************
    % Again, crosscorrelation with the now time-aligned signal
    % *******************************************************
    coeff=sum(signal .* conj(shiftedRef)) / sum(shiftedRef .* conj(shiftedRef));
    shiftedRef=shiftedRef * coeff;

    if forceReal
        shiftedRef = real(shiftedRef);
    end
end