Aligning two signals with subsample-accuracy

Summary

A Matlab / Octave function to align two signals with subsample accuracy

Introduction

Delay of a signal can be implemented as phase shift in the frequency domain. The mechanism is useful both to estimate the delay between a signal and a reference, and to time-shift a signal itself.
The details of the implementation are described in my blog article. This page contains a somewhat extended version of the function, that also works with complex signals.
The function assumes cyclic signals. This is usually not a limitation, because signals can be zero-padded as needed.
The function takes a signal and a reference.
It estimates and compensates:

The delay (full and fractional samples
An arbitrary scaling factor (real-valued or complex)

Usage

shifted_and_scaled_reference=mn_fitSignal_FFT(signal, reference);
errorsignal=signal-shifted_and_scaled_reference;

Download

mn_fitSignal_FFT.m

Complete example

Download
Two shifted pulses (“signal”, “reference”) are generated from the same function.
The fit_signal_FFT function shifts and scales the reference (blue), so that it matches the signal (green) optimally. The remaining difference (red) is very small.

Example output

close all; clear all;
addpath('../matlab/mn_blocks');
st=0.1;

pbase1=-10:st:10;
p1=exp(1.1+1.2i)*exp(-pbase1 .^2 +1i*pbase1);

pbase2=pbase1+13.4567*st;
p2=exp(1.3+1.4i)*exp(-pbase2 .^2 +1i*pbase2);

[coeff, p3]=mn_fitSignal_FFT(p1, p2);
figure(); hold on;
plot(abs(p1), 'g');
plot(abs(p2), 'b');
plot(abs(p3), 'k+'); 
plot(abs(p1-p3), 'r+');

legend('signal', 'reference', 'shifted and scaled reference', 'remaining error');
sum(abs(p1-p3))

Source code

% *******************************************************
% © Markus Nentwig, 2007
% This file is provided WITHOUT ANY WARRANTY
% *******************************************************

% returns complex coeff that minimizes the remaining energy
% in signal-coeff*ref
% returns also shifted and scaled reference, so that
% signal-shiftedRef is minimal
function [coeff, shiftedRef, deltaN]=mn_fitSignal_FFT(signal, ref)
  n=length(signal);
  
  % *******************************************************
  % 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 ...]. It cannot be delayed. The frequency component
    % is 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));
  % 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.
  % Also allow for a constant term to account for a constant
  % phase rotation in the time domain.
  % *******************************************************
  weight=abs(u); 
  constRotPhase=1 .* weight;
  uDelayPhase=uDelayPhase.* weight;
  ang=angle(u).* weight;
  r=transpose([constRotPhase; uDelayPhase])\transpose(ang); %linear mean square
							    
  %rotPhase=r(1); % constant phase rotation, not used here, obtained through coeff below; 
  fractionalDelay=r(2);
  
  % *******************************************************
  % Finally, the total delay is the sum of integer part and
  % fractional part.
  % *******************************************************
  deltaN=integerDelay+fractionalDelay;

  % *******************************************************
  % time-align the reference (phase shift in frequency domain)
  % *******************************************************
  % Fourier transform of a 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
  % gives the scaling factor
  % *******************************************************
  coeff=sum(signal .* conj(shiftedRef)) / sum(shiftedRef .* conj(shiftedRef));
  shiftedRef=shiftedRef * coeff;

© Markus Nentwig 2007-2008
The content of this page is provided without any warranty and may not be reproduced without permission.

Comments? Questions?

Please send me a mail! mnentwig@elisanet.fi