Code

High accuracy FFT interpolation

Jesus Selva November 4, 20103 comments Coded in Matlab

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);