## 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
```