## High accuracy FFT interpolation

The zero-padding FFT algorithm is a well-known method for applying a continuous delay to a signal. Yet it has the drawback that its accuracy is poor, since the interpolation error only decreases linearly with the number of samples.

This code snippet presents a simple modification of this algorithm that solves this accuracy problem. The proposed method has the following steps:

(1) Compute the FFT of the signal samples,

(2) repeat some of the spectral samples periodically,

(3) multiply these samples by some fixed weights,

(4) add zero padding, and

(5) compute the inverse FFT.

So, relative to the zero-padding FFT algorithm, this new method adds steps (2) and (3). The fact is that in this new method the interpolation error decreases EXPONENTIALLY. So, it is easy to obtain an accuracy below any numerical precision from a small number of samples. This new method is explained in detail in the paper

"Convolution-based trigonometric interpolation of band-limited signals", IEEE Transactions on Signal Processing, vol. 56, n. 11, pp. 5465-5477, Nov. 2008

The code that follows generates figure 2 in this paper. The proposed method is in function 'ProposedFFTUpsampling'. There you can clearly see the implementation of steps (1) to (5) above. It is worth comparing the zero-padding FFT algorithm in 'FraserFFTUpsampling' and the new method in 'ProposedFFTUpsampling'.

The code takes a few minutes to run due to the computation of the spectral coefficients for step (3). However, these coefficients must be computed only once ever for fixed T, B, P, and N.

```
function DemoProposed
% Author: J. Selva.
% First version: Sept. 2007.
%
% This m file contains the code for generating Fig. 2 in Sec. IV.A of paper
% "Convolution-based trigonometric interpolation of band-limited signals", IEEE
% Transactions on Signal Processing, vol. 56, n. 11, pp. 5465-5477, Nov. 2008. The
% equation numbers and references in this m file refer to this publication.
%
% Some comments:
%
% -The function 'KnabAPPeriodicApprox' computes the coefficients G(k/NT) in Sec. IV. This
% function takes several minutes, given that it involves a large number of numerical
% integrals which are performed in function KnabAPSpectrum. However, this must be done
% only once for fixed T, B, P and N, given that the G(k/NT) are independent of the sample
% values.
%
% 'RealTestSignal' implements a band-limited signal which consists of a sum of sinc pulses
% with arbitrary amplitudes and delays.
%
% 'FraserFFTUpsampling' implements the well-known zero-padding FFT
% algorithm. 'ProposedFFTUpsampling' implements the upsampling method in Sec. IV. It is
% worth comparing the code of both functions below.
warning off
T = 1; %Sampling period.
B = 0.5; %Two-sided bandwidth.
P = 18; %Pulse semi-length.
N = 1024; %Number of samples.
%Test signal composed of N+200 sinc pulses with random amplitudes and delays.
Sig = @(t)RealTestSignal(t,B,T,[-100,N+100],N+200);
%Compute coefficients Gn. This takes a while (5 minutes) because it is necessary to
%compute a large number of numerical integrals. However, since this step does not depend
%on the signal samples, it can be performed off-line and once ever.
disp(['Computing spectral coefficients. This takes a few minutes (7 or less), but ', ...
'must be done only',sprintf('\n'),'once ever for fixed T, B, P, and N.'])
% 'KnabAPPeriodicApprox' returns the frequency grid, specified by the first frequency fo
% and the frequency increment Df, the coefficients Gn, one for each grid frequency, and
% the truncation index kg.
[fo,Df,Gn,kg] = KnabAPPeriodicApprox(T,B,P,N);
sRef = Sig((0:N-1)*T); %Vector of signal samples with spacing T.
a = 10; %Oversampling factor.
sUpExact = Sig((0:a*N-1)*T/a); %Exact signal samples with oversampling a.
sUpFraser = FraserFFTUpsampling(sRef,a); %Interpolated samples using a zero-padded FFT.
sUpProposed = ProposedFFTUpsampling(sRef,T,B,12,a,Gn,kg); %Interpolated samples using
%the proposed method.
to = 0; %First instant in time grid.
Dt = T/a; %Grid time increment.
%Plot the interpolation error in dB for the zero-padding FFT method and for the proposed
%method.
plot(to+(0:a*N-1)*Dt,20*log10(abs([sUpExact - sUpFraser;sUpExact - sUpProposed])))
xlabel('Time t/T')
ylabel('Interpolation error (dB)')
annotation(gcf,'textarrow',[0.2911 0.2536],[0.8119 0.7429],...
'String',{'This is the zero-padding','FFT error'});
annotation(gcf,'arrow',[0.3321 0.06964],[0.4833 0.6333]);
annotation(gcf,'textarrow',[0.3446 0.2643],[0.419 0.3548],...
'String',{'This is the error of the','proposed FFT method. Note',...
'that the scale is in dB !!'});
warning on
%=====================================================================
%
% Conventional zero-padding FFT algorithm
function sUp = FraserFFTUpsampling(s,a)
N = length(s);
sF = a * fftshift(fft(s)); %Compute the FFT of the signal samples.
if mod(N,2) == 0 %Pad with zeros to obtain over-sampling factor a.
sF(end) = 0.5*sF(end);
sF = [zeros(1,N*(a-1)/2),sF,zeros(1,N*(a-1)/2)];
else
sF = [zeros(1,(N*(a-1)+1)/2),sF,zeros(1,(N*(a-1)-1)/2)];
end
sUp = ifft(ifftshift(sF)); %Compute the inverse FFT.
%=====================================================================
%
% Proposed FFT algorithm
function sUp = ProposedFFTUpsampling(s,T,B,P,a,Gn,kg)
N = length(s);
sF = fft(s); %Compute the FFT of the signal samples.
sF = sF(1+mod(-kg:kg,N)) .* Gn * (N*a); %Repeat some of the spectral samples periodically,
%and apply the window specified by the vector Gn.
sF = sF(:).';
if mod(N,2) == 0 %Pad with zeros to obtain over-sampling factor a.
sF = [zeros(1,N*a/2-kg),sF,zeros(1,N*a/2-kg-1)];
else
sF = [zeros(1,(N*a-2*kg-1)/2),sF,zeros(1,(N*a-2*kg-1)/2)];
end
sUp = ifft(ifftshift(sF)); %Compute the inverse FFT.
%=====================================================================
%
% Returns the spectral coefficients Gn and the frequency grid specification f0, Df, Gn.
function [f0,Df,Gn,n0] = KnabAPPeriodicApprox(T,B,P,N)
Df = 1/(N*T);
n0 = floor((1/T-B/2)/Df);
f0 = -Df*n0;
Bg = 2/T-B;
Gn = KnabAPSpectrum(f0+Df*(0:2*n0),B,Bg,(P+1/2)*T)/N;
%=====================================================================
%
% Auxiliary function that computes the spectral coefficients Gn.
function G = KnabAPSpectrum(f,B,Bg,TL)
% B: Two-sided signal bandwidth.
% Bg: Two-sided bandwidth of output pulse.
% Be: Two-sided bandwidth of fast-decay pulse.
% TL: Window's width between zeros in time domain.
%
% The two-sided bandwidth of the output pulse is Bg = B+2*Be.
Nf = numel(f);
G = zeros(size(f));
Be = (Bg-B)/2;
Ba = (Bg+B)/2;
Int = @(x)ConvWindow(x,Be,2*TL);
for k = 1:Nf
if f(k)-Ba/2 <= -Be/2 && Be/2 <= f(k)+Ba/2
G(k) = 1;
else
IInt = IntervalIntersection([-Be/2,Be/2],f(k)+[-Ba/2,Ba/2]);
if ~isempty(IInt)
G(k) = quad(Int,IInt(1),IInt(2),1e-16);
end
end
end
% Auxiliary function. It computes samples of the continuous convolution of a square pulse
% with the Kaiser window.
function W = ConvWindow(f,B,TL)
%Kaiser window.
%B: Two-sided window width.
%TL: Width between zeros in time domain.
W = zeros(size(f));
I = abs(f) < B/2;
W(I) = besseli(0,pi*B*(TL/2)*sqrt(1-(2*f(I)/B).^2))/(B*sinc(j*B*TL/2));
function [IInt,RelAInd,RelBInd] = IntervalIntersection(IA,IB)
if IA(2) < IB(1) || IB(2) < IA(1) %No overlap case.
IInt = [];
RelAInd = [];
RelBInd = [];
return
end
SwapIntervals = diff(IA) > diff(IB);
if SwapIntervals %Force IA to be the sorter interval.
Ipr = IA;
IA = IB;
IB = Ipr;
end
if IA(1) < IB(1)
IInt = [IB(1),IA(2)]; %There is intersection to the left.
RelAInd = IB(1)-IA(1);
RelBInd = 0;
if SwapIntervals
Relpr = RelAInd;
RelAInd = RelBInd;
RelBInd = Relpr;
end
return
end
if IA(2) <= IB(2)
IInt = IA; %IA lies inside IB.
RelAInd = 0;
RelBInd = IA(1)-IB(1);
if SwapIntervals
Relpr = RelAInd;
RelAInd = RelBInd;
RelBInd = Relpr;
end
return
end
IInt = [IA(1),IB(2)]; %There is intersection to the right.
RelAInd = 0;
RelBInd = IA(1)-IB(1);
if SwapIntervals
Relpr = RelAInd;
RelAInd = RelBInd;
RelBInd = Relpr;
end
%=====================================================================
%
% Computes samples of a signal formed by sincs with random amplitudes and delays.
function s = RealTestSignal(t,B,T,I,Np)
Sr = rand('seed');
rand('seed',0); %Always select the same random numbers in the sequel.
M = numel(t);
Ap = (rand(Np,1)-0.5)*2 ; %Random amplitudes.
tp = I(1)+rand(Np,1)*diff(I); %Random delays.
rand('seed',Sr)
st = size(t);
t = t(:);
s = zeros(numel(t),1);
for r = 1:Np
s = s + sinc(B*(t-tp(r)))*Ap(r); %Compute sum of sincs
end
s = reshape(s,st);
```