Efficient update of interpolation FIR coefficients

A common task  in many applications  is to apply  an arbitrary continuous  time shift to a
signal. For a fixed shift this can be easily done in Matlab using a filter design function
like  cremez. However,  if  the shift changes  from  sample to  sample,  then  it is  very
inefficient to re-compute the coefficients. The code that follows  presents a simple trick
to  update the  coefficients in  a   few operations  based   on the so-called  barycentric
interpolation.

In short,  suppose  we have computed  a  vector of  FIR  coefficients "cRef"  that is  able to
introduce a time shift "tRef", and assume that we would like to re-use it for a new time shift
"t". Then, we can do the following,

(1) Define a vector "tg"  containing the instants corresponding to the coefficients in
"cRef", i.e, the sampling instants. This could be a vector like [-P,-P+1,...,0,1,2,...P].

(2) Force the sum of the coefficients cRef to be equal to one, i.e, in Matlab do

  A = sum(cRef);
  c = cRef / A; %c will contain the new coefficients.

(3) This is the trick: adjust each coefficient using a barycentric shift. In Matlab,

  c = c .* (tRef-tg) ./ (t-tg);

(3) Make c have total sum equal to A

  c = A * c / sum(c);

That's all!

The explanation as to why it works so well can be found in

  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.

The code that follows evaluates the performance of this procedure by measuring the maximum spectral ripple. It takes a few minutes to run because it computes the optimal Remez coefficients for a large number of time shifts u. The output figure shows that the error of the barycentric coefficients is the same as that of the equi-ripple coeficients, in a range of length equal to two sampling periods.

This kind of interpolation has additional advantages,

- It is accurate even if the sampling instants are not uniformly spaced.

- It allows one to obtain also the derivatives of any order with low computational cost.

function BaryExample

%Author: J. Selva. 2010.
%
%This function compares the performance of the bary-shifted coefficients with that of the ...
%optimal equirriple coefficients. For a complete discussion see

%  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.

%This code was used to generate the example in Sec. IV.A of this paper.

T = 1; %Sampling period.
B = 0.7; %Two-sided bandwidth.
P = 5; %Filter semi-length. The total length is 2*P+1.

Lf = 0:B/400:B/2; %Frequency grid for evaluating the maximum ripple.

%Reference time shift, normalized to T.
unRef = 0.5;

%Reference filter coefficients. They introduce a fractional shift 0.5*T.

cRef = cremez(2*P,2*[0,B/2],{@fresp,T,P,T*unRef});

%Instants corresponding to each sample, relative to the central sample.
tg = fliplr((-P:P)*T);

%Grid of normalized (to T) time shifts in which the maximum ripple will be evaluated.
Lun = -P-2:.01:P+2;

%These variables will contain the maximum ripple for each time shift.
LBary = zeros(length(Lun),1);
LRemez = zeros(length(Lun),1);

for k= 1:length(Lun)
  LBary(k) = ErrorInFreq(BarycentricShift(cRef,tg,unRef,Lun(k),T),T,P,Lf,Lun(k));
  LRemez(k) = ErrorInFreq(cremez(2*P,2*[0,B/2],{@fresp,T,P,T*Lun(k)}),T,P,Lf,Lun(k));
end

%Plot the results. 

plot(Lun,mag2db(LRemez),'--',Lun,mag2dB(LBary),'-')

set(gca,'YLim',[-100,0])
xlabel('Normalized time shift (t/T)')
ylabel('Maximum spectral ripple (dB)')
grid on

annotation(gcf,'textarrow',[0.339285714285714 0.251785714285714],...
    [0.861904761904766 0.883333333333333],'TextEdgeColor','none',...
    'TextBackgroundColor',[1 1 1],...
    'String',{'Performance of bary-shifted','coefficients'});

annotation(gcf,'textarrow',[0.521428571428571 0.521428571428571],...
    [0.565666666666667 0.495238095238095],'TextEdgeColor','none',...
    'TextBackgroundColor',[1 1 1],...
    'String',{'For time shifts in [-T,T]','the performance is roughly','the same'});

annotation(gcf,'textarrow',[0.223214285714283 0.214285714285714],...
    [0.577571428571438 0.864285714285716],'TextEdgeColor','none',...
    'TextBackgroundColor',[1 1 1],...
    'String',{'Performance of optimal','equi-ripple coefficients'});

%=========================================================================================

%Applies the barycentric shift.

%cRef is the vector of FIR coefficients that interpolate at instant unRef*T from samples ...
%at instants in tg. 
%
%un*T is the desired interpolation instant. 
% 
%c is the new set of FIR coefficients.

function c = BarycentricShift(cRef,tg,unRef,un,T)

tRef = unRef * T;
t = un * T;

A = sum(cRef);
cRef = cRef / A;

c = cRef .* (tRef-tg) ./ (t-tg);

c = A * c / sum(c);

%=========================================================================================

%Evaluates the maximum spectral ripple of the FIR interpolator. For large P, it is more 
%efficient to employ the FFT.

function E = ErrorInFreq(c,T,P,Lf,un)

E = 0;
for k = 1:length(Lf)
  
  s0 = exp(1i*2*pi*Lf(k)*T*(-P:P));
  
  vRef = exp(1i*2*pi*Lf(k)*T*un);
  v = c*fliplr(s0).';
  E = max([E, abs(v-vRef)]);
end

%Auxiliary function for cremez.

function [DH,DW] = fresp(N,F,GW,W,T,P,u)

if ischar(N)
  DH = 'real';
  DW = [];
  return
end
  
DH = exp(1i*2*pi*(u-P*T)*GW/2);

DW = ones(size(DH));