DSPRelated.com
Code
   Code Snippets    How to speed up the sinc series

How to speed up the sinc series

Jesus Selva ● December 7, 2011 ● Coded in Matlab

Though the sinc series is the key for processing signals in the discrete domain, its interpolation performance is poor due to the slow decrease rate of the sinc tails. The solution to this problem has been available for decades but is fairly unknown in the signal processing field.

(This description is also available in pdf format)

The solution in the sequel is an implementation of the method in

 J. J. Knab, ``Interpolation of band-limited functions using the approximate prolate
 series,'' IEEE Trans. Inf. Theory, vol. IT-25, pp. 717-720, Nov. 1979.

Given a signal s(t) with two-sided bandwidth B, the Sampling Theorem provides the reconstruction formula

\tiny \dpi{200} s(nT+u)=\sum_{p=-\infty}^\infty s((n-p)T)\mathrm{sinc}(u/T+p),      (1)

where T is the sampling period, n is an integer, and u is the shift relative to the sampling grid of t, i.e, t=nT+u with -T/2 ≤ u <T/2. If (1) is truncated at indices ±P,

\tiny \dpi{200} s(nT+u)\approx\sum_{p=-P}^P s((n-p)T)\mathrm{sinc}(u/T+p),     (2)

the accuracy is poor due to the tails of the neglected sinc pulses. To improve accuracy, we can make good use of the sampling inefficiency as follows.

Consider a pulse w(t) whose main time lobe in inside the range ]-PT-T/2,PT+T/2[, i.e, w(t) is approximately zero if |t|≥PT+T/2, where P is the index at which we plan to truncate the series. For a small delay d, |d|≤T/2, this  condition implies that w(t-d) is approximately zero whenever |t|≥ (P+1)T. Additionally, if w(t) has two-sided bandwidth Bw and (B+Bw)T<1, then (1) can be applied to s(nT+u)w(u-d), i.e,


\tiny \dpi{200} s(nT+u)w(u-d)=\sum_{p=-\infty}^\infty s((n-p)T)w(-pT-d)\mathrm{sinc}(u/T+p).      (3)

and the truncated formula is

\tiny \dpi{200} s(nT+u)w(u-d)\approx\sum_{p=-P}^P s((n-p)T)w(-pT-d)\mathrm{sinc}(u/T+p).       (4)

In contrast with (2), this formula is accurate since the tails of the neglected terms have been damped by the samples w(-pT-d), |p|>P. Finally, if w(0)=1 and we set d=u, we get

\tiny \dpi{200} s(nT+u)\approx\sum_{p=-P}^P s((n-p)T)w(-pT-u)\mathrm{sinc}(u/T+p).     (5)

This is the accurate interpolator that we were looking for. To simplify notation, define the pulse

\tiny \dpi{200} g(t)=w(-t)\mathrm{sinc}(t/T).           (6)

Then (5) reads

\tiny \dpi{200} s(nT+u)\approx\sum_{p=-P}^P s((n-p)T)g(pT+u).       (7)


This formula says that we must replace the sinc pulse by g(t) to obtain accuracy. As to the pulse w(t), the one used in Knab's paper is the inverse Fourier transform of the Kaiser-Bessel window,

\tiny \dpi{200} w(t)=\frac{\mathrm{sinc}\Big((1-BT) \sqrt{(t/T)^2-P^2}\Big)} { \mathrm{sinc}(j(1-BT) P\big)}.

This window has an excellent performance. Actually, while the error for the sinc pulse decreases as 1/P, that of g(t) decreases as exp(-π P(1-BT)) !!. So, for the latter, small P values give very high accuracy. This behavior is due to the fact that w(t) is an approximation of the first prolate wave function, which is the band-limited signal with the highest energy concentration in a finite interval.
 
The following code tests the accuracy of g(t) versus that of sinc(t/T) for several P.

function Demo

close all

T = 1; %Sampling period.
B = 0.7; %Two-sided bandwidth.

%Plot error spectrum for P = 10, 20, and 30.

PlotErrorSpectrum(T,B,10);

figure

PlotErrorSpectrum(T,B,20);

figure 

PlotErrorSpectrum(T,B,30);

function PlotErrorSpectrum(T,B,P)

% Frequency grid specification in figure.
f0 = 0;
M = (2*P+1)*10;
Df = B/(2*M);

f = f0+Df*(0:M-1);

ESinc = zeros(size(f));
Eg = zeros(size(f));

%Look for the maximum error among different shifts u for each frequency in f. 
for u = -0.5:1/((2*P+1)*5):0.5

  HRef = exp(1i*2*pi*f*u); %Ideal spectrum.
  HSinc = chirpzDFT(sinc((-P:P)+u/T),-P*T,T,f0,Df,M); %Sinc spectrum.
  Hg = chirpzDFT(gPulseCoefs(u,T,B,P),-P*T,T,f0,Df,M); %g-pulse spectrum.

  ESinc = max([ESinc;abs(HRef-HSinc)]); %Compute current max. error for sinc.
  Eg = max([Eg;abs(HRef-Hg)]); %Compute current max. error for g. 

end

plot(f,20*log10([ESinc;Eg].'))

xlabel('Freq. (f*T)')
ylabel('Maximum error (dB)')

legend('sinc pulse','g pulse','Location','NorthOutside')

pr = get(gca,'YTick');

text(0,2*pr(end)-pr(end-1),['B*T = ',num2str(B*T),', P = ',num2str(P)])

title('Error spectrum (maximum for all possible shifts u)')

grid on

function c = gPulseCoefs(u,T,B,P)

t = (-P:P)*T+u;
c = sinc(t/T).*real(sinc((1-B*T)*sqrt((t/T).^2-P^2)))/real(sinc(1i*(1-B*T)*P));

function X = chirpzDFT(x,to,Dt,fo,Df,M)

%Author: J. Selva.
%Date: November 2011.
%
%This function computes the DFT for two regular time and frequency grids with arbitrary
%starting points and spacings, using the Chirp-z Transform. Specifically, it computes
%
%               N
%       X(k) = sum  x(n)*exp(-j*2*pi*(fo+Df*(k-1))*(to+Dt*(n-1))), 1 <= k <= M.
%              n=1
%
%Input parameters:
%
%  x: Signal samples. 
% to: Instant of first sample in x.
% Dt: Time spacing between consecutive samples of x.
% fo: First frequency in which the Fourier sum is computed.
% Df: Spacing of the desired frequency grid.
%  M: Desired number of output samples. 
%
%The algorithm is explained in Sec. 9.6.2, page 656 of
%
% A.V. Oppenheim, R. W. Schafer, J. R. Buck, "Discrete-time signal processing," second
% edition, 1998.
%

x = x(:).';

N = numel(x);
P = 2^ceil(log2(M+N-1));

%The next three lines do not depend on the vector x, and so they can be pre-computed if 
%the time and frequency grids do not change, (i.e, for computing the transform of 
%additional vectors). Thus, this algorithm just involves two FFTs for fixed grids and 
%three if the grids change. 

a = exp(-1i*2*pi*((1:N)*Dt*(fo-Df)+Dt*Df*(1:N).^2/2));
b = exp(-1i*2*pi*((to-Dt)*(fo+(0:M-1)*Df)+Dt*Df*(1:M).^2/2));
phi = fft(exp(1i*2*pi*Dt*Df*(1-N:M-1).^2/2),P); %FFT of chirp pulse. 

%Weigh x using a and perform FFT convolution with phi. 
X = ifft( fft(x.*a,P) .* phi );

%Truncate the convolution tails and weigh using b. 
X = X(N:M+N-1) .* b;