Interpolation from Nonuniform Samples

This snippet allows one to interpolate band-limited signals from a set of nonuniform samples with high accuracy. To see how it works see the paper

[1] J. Selva, "Functionally weighted Lagrange interpolation of bandlimited
signals from nonuniform samples," IEEE Transactions on Signal Processing,
vol. 57, no. 1, pp. 168-181, Jan 2009.

Also, some part of the code is explained in

[2] J. Selva, "Design of barycentric interpolator for uniform and nonuniform
 sampling grids", IEEE Trans. on Signal Processing, vol. 58, n. 3,
 pp. 1618-1627, March 2010.

To test the snippet, copy the code in a file named Demo.m and execute it. It is worth changing the truncation index P in line 20 to see how the accuracy and the interpolation range vary.

function Demo

%Author: J. Selva.
%Date: May 2011.

%The interpolation method in this snippet has been published in 
%
% [1] J. Selva, "Functionally weighted Lagrange interpolation of bandlimited
% signals from nonuniform samples," IEEE Transactions on Signal Processing, 
% vol. 57, no. 1, pp. 168-181, Jan 2009.
%
% Also, the barycentric implementation of the interpolator appears in 
%
% [2] J. Selva, "Design of barycentric interpolator for uniform and nonuniform 
% sampling grids", IEEE Trans. on Signal Processing, vol. 58, n. 3, 
% pp. 1618-1627, March 2010.

T = 1;   %Sampling period
B = 0.7; %Two-sided signal bandwidth. It must be B*T < 1.
P = 12;  %Truncation index. The number of samples taken is 2*P+1. The error 
         %decreases EXPONENTIALLY with P as exp(-pi*(1-B*T)*P). 
 	%(Try different values of P like 3, 7, 20, and 30.)

int = [-30,30]*T; %Plot x limits.

JitMax = 0.4; %Maximum jitter. This must be between 0 and 0.5*T. The samples
              %are taken at instants p*T+d[p], -P<=p<=P and 
	      %-JitMax <= d[p] <= JitMax.  

while 1
  
  st = TestSignal(int,ceil(60/B),B); %This defines a test signal formed 
              %by a random sum of sincs of bandwidth B. 
  
  d = 2*(rand(1,2*P+1)-0.5)*JitMax; %Vector of random jitter.

  disp('Sampling instants (t_k/T):')
  disp((-P:P)*T+d)
  
  sg = TestSignal(st,(-P:P)*T+d); %Sample signal at jittered instants

  t = int(1):0.1*T:int(2); %Time grid for figure.

  sRef = TestSignal(st,t); %Exact signal samples in grid t.
  sInt = NonuniformInterp(t,sg,d,B,T); %Signal samples interpolated from 
                           %d (jitter values) and sg (values at jittered 
			   %instants).

  subplot(2,1,1)
  plot(t,sRef,t,sInt) %Plot exact and interpolated signals

  xlabel('Normalized time (t/T)')
  title('Signal (blue) and interpolator (green)')
  grid on

  subplot(2,1,2)
  plot(t,20*log10(abs(sRef-sInt))) %Plot interpolation error in dB.
  xlabel('Normalized time (t/T)')
  ylabel('Interpolation error (dB)')
  grid on
  
  R = input('Press Enter for a new trial, q to quit: ','s');

  if any(strcmp(R,{'q','Q'}))
    break
  end

end

function v = NonuniformInterp(u,sg,d,B,T)

%Performs interpolation from nonuniform samples.
% T   is the period of the reference grid (instants (-P:P)*T).
% d   is the vector of jitter values relative to the reference grid. It must have an odd ...
%     number of elements. 
% B:  signal's two-sided bandwidth. It must be B*T<1. The error decreases exponentially as
%     exp(-pi*(1-B*T)P).
% u:  Vector of instants where the signal is to be interpolated.
% sg: Sample values at instants (-P:P)*T+d.

if mod(length(d),2) ~= 1
  error('The number of samples must be odd')
end

P = (length(d)-1)/2;

tg = (-P:P)*T+d(:).';

w = baryWeights(tg,T,B,P);

v =  DerBarycentricInterp3tVecV1(w,sg,tg,u,0);

function v = APPulse(t,B,TSL)

v = real(sinc(B*sqrt(t.^2-TSL^2)))/real(sinc(1i*pi*B*TSL));

function w = baryWeights(vt,T,B,P)

%Barycentric weights. See Eq. (23) and sequel in [2].

g = APPulse(vt,1/T-B,T*(P+1));
gam = gamma(vt/T+P+1) .* gamma(-vt/T+P+1) .* g;

N = length(vt);
LD = ones(1,N);

for k = 1:N
  LD(k) = prod(vt(k)-vt(1:k-1))* prod(vt(k)-vt(k+1:N));
end

w = gam ./ LD;
w = w / max(abs(w));

function out = DerBarycentricInterp3tVecV1(alpha,s,t,tau,nd)

%Vectorized implementation of barycentric interpolator in [2, Sec. IV].

s = s(:).';
t = t(:).';
sztau = size(tau);
tau = tau(:);

Nt = numel(t);
Ntau = numel(tau);

vD = zeros(Ntau,1);
LF = [ones(Ntau,1),zeros(Ntau,Nt-1)];
out = zeros(Ntau,nd+1);

for k = 1:Nt-1
  vD = vD .* (tau-t(k))+alpha(k)*LF(:,k);
  LF(:,k+1) = LF(:,k) .* (tau-t(k));
end

vD = vD .* (tau-t(Nt)) + alpha(Nt) * LF(:,Nt);

z = s(ones(Ntau,1),:);

for kd = 0:nd

  pr = z(:,end);
  z1 = z-pr(:,ones(1,Nt)); cN = zeros(Ntau,1);

  for k = 1:Nt-1
    cN = cN .* (tau-t(k))+z1(:,k)*alpha(k) .* LF(:,k);
  end
  cN = cN ./ vD;
  
  ztau = z(:,end)+(tau-t(end)) .* cN;
  out(:,kd+1) = ztau;
  
  if kd < nd
    z = [ (z(:,1:end-1)-ztau)./(t(ones(Nt,1),1:end-1)-tau(:,ones(1,Nt))) , ...
	  cN ] * (kd+1);
  end
  
end

out = reshape(out,sztau);

function varargout = TestSignal(varargin)

%Test signal formed by a random sum of sincs. 

if nargin == 3
  [int,Nsinc,B] = deal(varargin{:});
  st = [];
  st.B = B;
  st.Ampl = (rand(1,Nsinc)-0.5)*2;
  st.Del = int(1)+(int(2)-int(1))*rand(1,Nsinc);
  varargout = {st};
else
  [st,t] = deal(varargin{:});
  v = zeros(size(t));
  for k = 1:numel(t)
    v(k) = st.Ampl * sinc(st.B*(t(k)-st.Del)).';
  end
  varargout = {v};
end