DSPRelated.com
Code
   Code Snippets    Submitted by jselva

How to speed up the sinc series

Jesus Selva ● December 7, 2011 ● Coded in Matlab 
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;

FFT for arbitrary time and frequency grids (Chirp-z Transform)

Jesus Selva ● November 17, 2011 ● Coded in Matlab 
function Demo

while 1
  
  %Frequency grid with arbitrary start value, spacing, and number of elements. ...
  %Mathematically, the grid is
  %
  %   fo+Df*m, m=0, 1,...,M-1.
  %
  
  fo = 10*rand(1)-0.5; %Start value.
  Df = rand(1)+0.01;   %Frequency spacing. 
  M = round(rand(1)*1000); %Number of frequencies. 

  %Time grid with arbitrary start value, spacing, and number of elements. Mathematically, ...
  %the grid is
  %
  %   to+Dt*n, n=0, 1,...,N-1.
  %

  to = 10*rand(1)-0.5; %Start value; (instant of first sample).
  Dt = rand(1)+0.01;   %Time spacing. 
  N = round(rand(1)*1000); %Number of samples. 
  
  %Random vector of samples.
  x = randn(1,N)+1i*randn(1,N);
  
  %x(n) is the sample at instant to+Dt*(n-1), n=1,2,...,N. 

  %We proceed to compute the DFT of x at frequencies fo+m*Df, m=0, 1,..., M-1.

  %Compute DFT using the direct and chirp-z methods.
  
  tic, XD = DirectGridDFT(x,to,Dt,fo,Df,M); tmD = toc;
  tic, XF = chirpzDFT(x,to,Dt,fo,Df,M); tmF = toc;
  
  disp(['Timing direct method: ', num2str(tmD), ' sec.'])
  disp(['Timing chirp z: ',num2str(tmF),' sec.'])
  disp(['Relative difference: ', num2str(max(abs(XD-XF))/max(abs(XD)))])
  disp(' ')
  input('(Press a key for another trial or Ctrl-C) ')
  disp(' ')
end

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;

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

%Direct (and slow) version of the Fourier sum with arbitrary and regular time and
%frequency grids. Used for validation of chirpzDFT.

x = x(:).';

N = length(x);

X = zeros(1,M);

for k = 1:M
  for n = 1:N
    X(k) = X(k) + x(n)*exp(-1i*2*pi*(to+Dt*(n-1))*(fo+Df*(k-1)));
  end
end

ML Estimation of Several Frequencies Using the Nonuniform FFT

Jesus Selva ● July 8, 2011 ● Coded in Matlab 
function st = Demo

% Author: J. Selva.
% Date: July 2011.
% Version: 0.

st = [];

n = input('Options: 1 perform single estimation, 2 estimate RMS error: ');

switch n
  case 1
    st = TestSingleTrial;
  case 2
    TestRMSError;
end

function st = TestSingleTrial

%Scenario with 5 frequencies. 

M = 2048;
SNRInSpecdB = 27; %Signal-to-noise ratio after correlation. 

Nf = 5;
df = 1/M; %Minimum separation among frequencies. 

MaxAmp = 1; %Maximum amplitude.
PFA = 10.^-5;
nVar = M*MaxAmp^2/10^(SNRInSpecdB/10);
nThres = NoiseThres(M,nVar,PFA);
MinAmp = 2*nThres/M;

[vf,va] = RandomFreqAndAmpVectors(Nf,df,MinAmp,MaxAmp);

y = DataRealization(vf,va,M,nVar);

p = Preprocessing(y,nVar,PFA);

disp('True frequencies: ')
disp(vf)

disp('True amplitudes: ')
disp(va)

disp(['Post-correlation Signal-to-noise ratio: ',...
      num2str(SNRInSpecdB),' dB'])

s = input('To compute ML estimate, use (S)erial  or (M)anual method? ','s');

if strcmp(lower(s),'s')
  st = MethodSerial1(p);
else
  st = MethodManual1(p);
end

disp('ML frequency estimates: ')
disp(st.Freqs)

disp('ML amplitude estimates: ')
disp(st.Amps)

function TestRMSError

%Estimates the RMS error and computes the Cramer-Rao bounds.

format compact
format short g

M = 2048;
SNRInSpecdB = 27;

Nf = 5;
df = 1/M;

MaxAmp = 1;
PFA = 10.^-5;
nVar = M*MaxAmp^2/10^(SNRInSpecdB/10);
nThres = NoiseThres(M,nVar,PFA);
MinAmp = 2*nThres/M;

[vf,va] = RandomFreqAndAmpVectors(Nf,df,MinAmp,MaxAmp);

crb = CRBBound(vf,va,M,nVar);

Ntr = 0;
S2 = zeros(Nf,1);
for k = 1:1e5

  y = DataRealization(vf,va,M,nVar);

  p = Preprocessing(y,nVar,PFA);

  st = MethodSerial1(p);
  
  if length(st.Freqs) == 5
    Ntr = Ntr + 1;
    S2 = S2 + (vf-st.Freqs).^2;
  end
  
  if mod(k,100) == 0
    disp(['Number or trials: ', num2str(k)])
    disp(['             Frequencies:      ' num2str(vf.')])
    disp(['Deviation CRB bounds (dB):   ' num2str(20*log10(crb.'))])
    disp(['          RMS errors (dB):   ' num2str(10*log10(S2(:).'/Ntr))])
    disp(' ');
  end
end

function [vf,va] = RandomFreqAndAmpVectors(Nf,MinFreqSep,MinAmp,MaxAmp)

%Constructs random vectors of frequencies and amplitudes. The module of the largest
%amplitude is MaxAmp. There are no amplitudes with module below MinAmp, and the minimum
%separation among frequencies is MinFreqSep.

vf = sort(rand(Nf,1));
vf = vf + MinFreqSep*(cumsum(ones(Nf,1))-1);
vf = vf/(1+(Nf-1)*MinFreqSep) - 0.5;

pr = rand(Nf,1);
ind = floor(rand(1)/Nf)+1;
pr(ind) = 1;

va = (MinAmp+(MaxAmp-MinAmp)*pr) .* exp(1i*2*pi*rand(Nf,1));

function y = DataRealization(vf,a,M,nVar)

m1 = FirstIndex(M);

y = exp(1i*2*pi*(m1:m1+M-1).'*vf(:).')*a(:)+randn(M,2)*[1 ; 1i]*sqrt(nVar/2);

function p = Preprocessing(y,nVar,PFA)

%Performs the preprocessing. Fundamentally this consists in computing the FFT inside 
%"Correlation.m".
%The input parameters are
%
% y: Data vector.
% nVar: Noise variance estimate. 
% PFA: Probability of seeing a frequency when there isn't any. Usually a very low value
%      like 1e-5.

p = [];
p.M = length(y);
p.nVar = nVar;
p.PFA = PFA;
p.nThres = NoiseThres(p.M,p.nVar,p.PFA); %Noise threshold.

icf = Correlation(y,p.nThres); %Struct with data for interpolating the correlation.

p.cfd = Cost(real(y(:)'*y(:)),icf); %Struct with data for computing the ML cost function ...
                                    %and its derivatives. 

p.freqThres = 1/(10*p.M); %If two frequencies differ in less than ...
                          %freqThres  in  the  iterative  search, then   one  of  them  is
                          %discarded. This prevents ill-conditioning.

p.freqStopThres = 1/(p.M*1e4); %Threshold for stopping the ML cost function minimization.

function st = MethodSerial1(p)

%Adds frequencies and computes ML estimates sequentially.

vfEst = [];
ModifiedEstimation = logical(1);

while ModifiedEstimation
  
  [vf1,va1] = Correlation(p.cfd.icf,'gridMainLocalMaxima',vfEst);

  ModifiedEstimation = ~isempty(vf1);

  if ModifiedEstimation
    vfEst(end+1) = vf1(1);
    [listv,listLv,listCL,ErrRatio] = LookForMinimum(p.cfd,vfEst,...
	10.^-8,p.freqThres,p.freqStopThres);
    vfEst = sort(listv(:,end));
  end 
end

st = Report(p,vfEst);

function st = MethodManual1(p,varargin)

%Allows the user to set the initial frequency estimates on a plot.

if nargin == 1
  vfEst = [];
else
  vfEst = varargin{1};
end

vfEst = vfEst(:);

figure
subplot(2,1,1)
Correlation(p.cfd.icf,'PlotPeriodogramFit',vfEst);
yl = get(gca,'YLim');
subplot(2,1,2)
Correlation(p.cfd.icf,'PlotPeriodogramFitError',vfEst);

while 1
  pr1 = input('A/a -> Add freq., D/d -> Delete freq, X/x -> exit: ','s');
  
  switch lower(pr1)
    case 'a'
      [f1,pr] = ginput(1);
      vfEst = sort([vfEst; f1]);
      [listv,listLv,listCL,ErrRatio] = LookForMinimum(p.cfd,vfEst,10.^-8,... 
	  p.freqThres,p.freqStopThres);
      vfEst = sort(listv(:,end));

      subplot(2,1,1)
      Correlation(p.cfd.icf,'PlotPeriodogramFit',vfEst);

      subplot(2,1,2)
      Correlation(p.cfd.icf,'PlotPeriodogramFitError',vfEst);

    case 'd'
      [f1,pr] = ginput(1);
      [pr,ind] = min(abs(f1-vfEst));
      vfEst(ind) = [];

      subplot(2,1,1)
      Correlation(p.cfd.icf,'PlotPeriodogramFit',vfEst);

      subplot(2,1,2)
      Correlation(p.cfd.icf,'PlotPeriodogramFitError',vfEst);

    case 'x'
      
      st = Report(p,vfEst);
      return
  end
end

function thres = NoiseThres(M,nVar,PFA)

%Computes the noise threshold. The method will discard any correlation local maximum below
%this threshold.

thres = sqrt((M*nVar/2)*icdf('chi2',(1-PFA)^(1/M),2));

function m1 = FirstIndex(M)

m1 = -ceil(M/2);

function varargout = Correlation(varargin)

%Computes the correlation value and its derivatives by calling "Barycentric". It does
%other things like plotting. 

if ~isstruct(varargin{1})
  st = [];
  temp = {[],[],10.^-15,1.5};
  [temp{1:nargin}] = deal(varargin{:});
  [a,st.nThres,epsilon,OvFMin] = deal(temp{:});
  a = a(:).';
  
  st.M = length(a);
  st.K = 2^nextpow2(st.M*OvFMin);
  st.bary = Barycentric(0,SampleFourierCorrelation1DVer1(a,st.K),...
           -2*FirstIndex(st.M),1/st.K,whatP(1/st.K,-2*FirstIndex(st.M),epsilon));
       

  st.IndActive = find(abs(st.bary.a) > st.nThres);
  st.IndActive = unique([st.IndActive-1,st.IndActive,st.IndActive+1]);
  if st.IndActive(1) == 0
    st.IndActive(1) = [];
  end

  if st.IndActive(end) == st.K+1
    st.IndActive(end) = [];
  end
  
  st.aTActive = st.bary.a(st.IndActive);
  st.vfActive = centreFreq((st.IndActive-1)/st.K);
  varargout = {st};
  return
end 

[st,msg] = deal(varargin{1:2});
  
switch msg
    
  case 'value'
    varargout = {Barycentric(st.bary,varargin{3:end})};
    
  case 'gridValueInRange'
      
    [ind1,ind2] = deal(varargin{3:4});
    varargout = {st.bary.a(mod(ind1:ind2,st.K)+1)};
     
  case 'ampEstimates'
      
    vf = varargin{3};
    if length(vf) == 0
      varargout = {[]};
    end

    varargout = {RegGeoSum(st.M,vf,vf,0,0) \ Barycentric(st.bary,vf)};
      
  case 'gridValue'
    
    lInd = varargin{3};
    if nargin == 4
      vf = varargin{4};
    else
      vf = [];
    end
      
    if isempty(vf)
      varargout = {st.bary.a(mod(lInd,st.K)+1)};
    else
      vf1 = centreFreq(lInd/st.K);
      
      VRaa = RegGeoSum(st.M,vf,vf,0,0);
      VRa1a = RegGeoSum(st.M,vf1,vf,0,0);
      
      pr = VRa1a*(VRaa\Barycentric(st.bary,vf));
      pr = reshape(pr,size(lInd));

      varargout = {st.bary.a(mod(lInd,st.K)+1)-pr};
    end

  case 'gridFit'
    
    lInd = varargin{3};
    if nargin == 4
      vf = varargin{4};
    else
      vf = [];
    end
      
    if isempty(vf)
      varargout = {zeros(size(lInd))};
    else
      vf1 = centreFreq(lInd/st.K);
      
      VRaa = RegGeoSum(st.M,vf,vf,0,0);
      VRa1a = RegGeoSum(st.M,vf1,vf,0,0);
      
      pr = VRa1a*(VRaa\Barycentric(st.bary,vf));
      pr = reshape(pr,size(lInd));

      varargout = {pr};
    end
    
  case 'plotErrorIndB'

    if nargin == 2
      vf = [];
    else 
      vf = varargin{3};
    end
      
    s = 20*log10(abs(Correlation(st,'gridValue',0:st.K-1,vf)));
    gf = centreFreq((0:st.K-1)/st.K);
    [gf,ord] = sort(gf);
    s = s(ord);
      
    plot(gf,s,[-0.5,0.5],20*log10(abs(st.nThres)*[1,1]))
      
    Ms = max(s);
    set(gca,'YLim',[Ms-50,Ms+10])
    xlabel('Frequency (f)')
    ylabel('Error spectrum (dB)')
    grid on
    
    varargout = {};
      
  case 'PlotPeriodogramFit'
    
    vfEst = varargin{3};
    
    if isempty(vfEst)
      
       s0 = 20*log10(abs(Correlation(st,'gridValue',0:st.K-1)));
       gf = centreFreq((0:st.K-1)/st.K);
       [gf,ord] = sort(gf);
       s0 = s0(ord);
      
       plot(gf,s0)
       hold on
       nt = 20*log10(abs(st.nThres));
       plot([-0.5,0.5],nt([1,1]).',...
	   'LineWidth',2,'LineStyle','--','Color',[1 0 0])
       hold off
       text(-0.49,nt+2,'Noise threshold')
       
       xlabel('Frequency (f)')
       ylabel('Amplitude (dB)')
       title('Periodogram data, fit, and est. frequencies (stems)')
       return
    end
    
    vaEst = Correlation(st,'ampEstimates',vfEst);
    
    stem(vfEst,20*log10(st.M*abs(vaEst)),'g')
    
    hold on 
    
    s0 = 20*log10(abs(Correlation(st,'gridValue',0:st.K-1)));
    gf = centreFreq((0:st.K-1)/st.K);
    [gf,ord] = sort(gf);
    s0 = s0(ord);
      
    plot(gf,s0,'b')
    hold on
    nt = 20*log10(abs(st.nThres));
    plot([-0.5,0.5],nt([1,1]).',...
	'LineWidth',2,'LineStyle','--','Color',[1 0 0])
      
    s1 = 20*log10(abs(Correlation(st,'gridFit',0:st.K-1,vfEst)));
    s1 = s1(ord);
      
    plot(gf,s1,'r')

    Ms = max(s0);
    set(gca,'YLim',[Ms-50,Ms+10])
    
    xlabel('Frequency (f)')
    ylabel('Amplitude (dB)')
    title('Periodogram fit')
    
    hold off
    
    varargout = {};
    
  case 'PlotPeriodogramFitError'
    
    if nargin == 2
      vfEst = [];
    else
      vfEst = varargin{3};
    end
    
    vaEst = Correlation(st,'ampEstimates',vfEst);
    
    s0 = 20*log10(abs(Correlation(st,'gridValue',0:st.K-1,vfEst)));
    gf = centreFreq((0:st.K-1)/st.K);
    [gf,ord] = sort(gf);
    s0 = s0(ord);
     
    nt = 20*log10(abs(st.nThres));
    
    plot(gf,s0)
    hold on 
    plot([-0.5,0.5],nt([1,1]).',...
	'LineWidth',2,'LineStyle','--','Color',[1 0 0])
    hold off
    
    Ms = max(s0);
    set(gca,'YLim',[Ms-50,Ms+10])
    
    xlabel('Frequency (f)')
    ylabel('Amplitude (dB)')
    title('Periodogram error')
    
    varargout = {};
    
  case 'gridValueGivenAmps'
      
    [lInd,vf,aEst] = deal(varargin{3:5});
    vf1 = centreFreq(lInd/st.K);
    VRa1a = RegGeoSum(st.M,vf1,vf,0,0);
    
    pr1 = st.bary.a(mod(lInd,st.K)+1);
    pr1 = reshape(pr1,size(lInd));
    pr2 = VRa1a*aEst;
    pr2 = reshape(pr2,size(lInd));
      
    varargout = {pr1-pr2};
      
  case 'gridLocalMaxima'
      
    if nargin == 2

      pr = [-2,-1,0].';
      
      pr = abs(st.bary.a(mod(st.IndActive([1;1;1],:) + ...
             pr(:,ones(1,length(st.IndActive))),st.K)+1));

      pr = pr(1,:)<=pr(2,:) & pr(3,:)<=pr(2,:);
      
      vf = st.vfActive(pr);
      va = st.aTActive(pr);
   
      varargout = {vf,va};
      return
    end
      
    vf = varargin{3};
      
    if length(vf) == 0
      [varargout{1:2}] = Correlation(st,'gridLocalMaxima');
      return
    end
   
    amp = Correlation(st,'ampEstimates',vf);

    tmp = abs(Correlation(st,...
	    'gridValueGivenAmps',st.IndActive-1,vf,amp)) >= st.nThres;
     
    ind = st.IndActive(tmp);
    
    pr = [-2,-1,0].';
    
    ind1 = ind([1,1,1].',:) + pr(:,ones(1,length(ind)));
    
    pr = abs(Correlation(st,'gridValueGivenAmps',ind1,vf,amp));

    pr = pr(1,:)<=pr(2,:) & pr(3,:)<=pr(2,:);
    
    ind = ind(pr)-1;
   
    if isempty(ind)
      varargout = {[],[]};
    else
      varargout = {centreFreq(ind/st.K),Correlation(st,...
	  'gridValueGivenAmps',ind,vf,amp)};
    end

  case 'gridMainLocalMaxima'
      
    [vf,va] = Correlation(st,'gridLocalMaxima',varargin{3:end});
    if length(vf) == 0
      varargout = {[],[]};
      return
    end
	  
    [varargout{1:2}] = selectLocalMaximaVer3(vf,va,st.nThres,st.M);
	
      
  case 'gridMaximum'
      
    pr = Correlation(st,'gridValue',0:st.K-1,varargin{3:end});
    [ma,ind] = max(abs(pr));
      
    varargout = {(ind-1)/st.K,ma};
      
end

function tab = SampleFourierCorrelation1DVer1(a,K)  

M = length(a);
m1 = FirstIndex(M);

a = a(:).';
tab = fft([a(-m1+1:end),zeros(1,K-M),a(1:-m1)]);

function out = centreFreq(f)

out = f + ceil(-0.5-f);

function P = whatP(T,B,epsilon)

P = ceil(acsch(epsilon)/(pi*(1-B*T)));

function d = dist(x,y)

if x<y
  d = min([y-x,x+1-y]);
else
  d = min([x-y,y+1-x]);
end

function v = sincMask(x,y,M)

v = zeros(size(x));
for k = 1:length(x)
  x1 = max([dist(x(k),y)-1/(2*M),0]);
  v(k) = 1/max([M*sin(pi*x1),1]); 
end

function v = slBoundVer2(x,vam,vfe,nThres,M1)

vam = vam(:);
vfe = vfe(:);

sll = sincMask(vfe,x,M1).'*vam;

v = max([sll+nThres,1.1*sll]);

function varargout = selectLocalMaximaVer3(vf,va,nThres,M)

va1 = abs(va(:));
[pr,ord] = sort(va1);
ord = flipud(ord);

vf1 = vf(ord);
va1 = va1(ord);

n1 = length(va1);

if va1(1) < nThres 
  varargout = {[],[]};
  return
end

va2 = va1(1);
va2Abs = va1(1);
vf2 = vf1(1);

for n=2:n1
  if va1(n)<nThres
    continue
  end
  if va1(n) < slBoundVer2(vf1(n),va2Abs,vf2,nThres,M)
    continue
  end
  vf2(end+1) = vf1(n);
  va2(end+1) = va1(n);
  va2Abs(end+1) = abs(va(ord(n)));
end

varargout = {vf2,va2};

function out = Barycentric(varargin)

%Performs  barycentric interpolation using the method in [Bary].

if ~isstruct(varargin{1})
  st = [];
  [st.n0,st.a,st.B,st.T,st.P] = deal(varargin{:});
  
  st.N = length(st.a);
  st.alpha = baryWeights(st.T,st.B,st.P);
  st.alpha = st.alpha([1:st.P,st.P+2:end,st.P+1]);
  
  st.vt = (-st.P:st.P)*st.T; 
  st.vt = st.vt([1:st.P,st.P+2:end,st.P+1]);

  out = st;
  
else
  
  [st,v] = deal(varargin{1:2});
  
  if nargin == 3
    nd = varargin{3};
  else
    nd = 0;
  end
 
  v = v(:);
  Nv = numel(v);

  n = floor(v/st.T+0.5);
  u = v - n*st.T;
  n1 = n - st.n0;
    
  pr = [-st.P:-1,1:st.P,0];
  da = st.a(mod(n1(:,ones(1,2*st.P+1)) + pr(ones(Nv,1),:),st.N)+1);
  
  out = DerBarycentricInterp3Vec(st.alpha,da,st.vt,u,nd);
     
end

function w = baryWeights(T,B,P)

vt = (-P:P)*T;
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 v = APPulse(t,B,TSL)

v = real(sinc(B*sqrt(t.^2-TSL^2)))/real(sinc(1i*pi*B*TSL));

function out = DerBarycentricInterp3Vec(alpha,zL,t,tauL,nd)

NtauL = size(tauL,1);
LBtau = 200; %Block size for vectorized processing. Change this to adjust the memory
             %usage.

NBlock = 1+floor(NtauL/LBtau);

t = t(:).';
Nt = length(t);

out = zeros(NtauL,nd+1);

for r = 0:NBlock-1
  ind1 = 1+r*LBtau;
  ind2 = min([(r+1)*LBtau,NtauL]);
 
  Ntau = ind2-ind1+1;
  z = zL(ind1:ind2,:);
  tau = tauL(ind1:ind2);
 
  vD = zeros(Ntau,1);
 
  LF = [1,zeros(1,Nt-1)];
  LF = LF(ones(Ntau,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);

  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(ind1:ind2,kd+1) = ztau;
 
    if kd < nd
      pr1 = ztau;
      pr1 = z(:,1:end-1) - pr1(:,ones(1,Nt-1));
     
      pr2 = t(1:end-1);
      pr2 = pr2(ones(Ntau,1),:)-tau(:,ones(1,Nt-1));
     
      z = [pr1 ./ pr2, cN] * (kd+1);
     
    end
  end
 
end

function varargout = Cost(varargin)

%Computes the ML cost function differentials by calling "Correlation" and "RegGeoSum". 

if ~isstruct(varargin{1})
  st = [];
  [st.EnergyRyy,st.icf] = deal(varargin{:});
  varargout = {st};
  return
end

st = varargin{1};

if nargin > 2
  vf = varargin{3};
else
  vf = [];
end

switch varargin{2}
  
  case 'value'
    
    vf = varargin{3};
    Rss = RegGeoSum(st.icf.M,vf,vf,0,0);
    Rsd = RegGeoSum(st.icf.M,vf,vf,0,1);
    Rdd = RegGeoSum(st.icf.M,vf,vf,1,1);
    
    pr = Correlation(st.icf,'value',vf,1);
    Rsy = pr(:,1);
    Rdy = pr(:,2);
    
    pr = Rss\[Rsy,Rsd];
    MRsy = pr(:,1);
    MRsd = pr(:,2:end);
    
    varargout = {st.EnergyRyy-RealTrace(Rsy',MRsy),-2*RealDiag(MRsy,Rdy'-MRsy'*Rsd),...
	2*RealHadamard((Rdd'-Rsd'*MRsd).',MRsy*MRsy')};

  case 'value0'
  
    if length(vf) == 0
      varargout = {st.EnergyRyy};
      return
    end

    vf = varargin{3};
    Rss = RegGeoSum(st.icf.M,vf,vf,0,0);
    Rsy = Correlation(st.icf,'value',vf);
    
    [varargout{1}] = st.EnergyRyy-RealTrace(Rsy',Rss\Rsy);
    
end

function v = RealHadamard(A,B)

v = real(A).*real(B)-imag(A).*imag(B);

function C = RealDiag(A,B)

% Same as real(diag(A*B)) but using less flops.

C = RealHadamard(A.',B);

if size(C,1) == 1
  C = C.';
else
  C = sum(C).';
end

function C = RealTrace(A,B)

% Same as real(trace(A*B)) but using less flops.

C=RealHadamard(A.',B);
C=sum(C(:));

function [listv,listLv,listCL,ErrRatio] = ...
    LookForMinimum(L,v0,nThres,freqThres,freqStopThres,varargin)

%Executes a descent algorithm using a method similar to that in [ML].

if ~isempty(varargin) 
  nItMax = varargin{1};
else
  nItMax = 20;
end

v0 = v0(:);
Nf = numel(v0);

listv = [v0,zeros(Nf,nItMax)];
listLv = [Cost(L,'value0',v0),zeros(1,nItMax)];
listCL = [NaN,zeros(1,nItMax)];

succ = 1;
nIt = 1;

while succ && nIt <= nItMax
  [succ,v,Lv,CL] = NewtonIt(L,listv(:,nIt),freqThres);
  
  if succ 
    nIt = nIt + 1;
    listv(:,nIt) = v;
    listLv(nIt) = Lv;
    listCL(nIt) = CL;
  else
    break
  end
    
  if max(abs(listv(:,nIt)-listv(:,nIt-1))) < freqStopThres
    break 
  end
end

ErrRatio = (Cost(L,'value0')-listLv(end))/listLv(end);
listv(:,nIt+1:end) = [];
listLv(nIt+1:end) = [];
listCL(nIt+1:end) = [];

function [v,vf] = validVector(vf,thres)

vf = sort(vf(:));
if length(vf)==1 
  v = 1;
  return
end

m = vf(2:end)-vf(1:end-1) < thres;

if any(m)
  v = 0;
  vf([logical(0); m]) = [];
else
  v = 1;
end

function [succ,v1,Lv1,nCutsLeft] = NewtonIt(L,v0,thres)
  
Lv1 = [];
nCutsLeft = [];

[vv,v1] = validVector(v0,thres);

if ~vv
  succ = 0;
  return
end

[Lv0,g0,h0] = Cost(L,'value',v0);

d = -h0\g0;

nCutsLeft = 20;

while nCutsLeft > 0
  [vv,v1] = validVector(v0+d,thres);
  if ~vv
    succ = 0;
    v1 = v0;
  end
  
  Lv1 = Cost(L,'value0',v0+d);
  
  if Lv1 < Lv0
    succ = 1;
    v1 = v0 + d;
    break;
  end
  
  nCutsLeft = nCutsLeft - 1;
  d = d/2;
end

if nCutsLeft == 0
  succ = 0;
end

function out = RegGeoSum(varargin)

%Computes the value of a geometric sum or its differentials using the exact formula. It
%takes care of the preventable singularity at zero.

switch nargin
  case 2
    [M,vf] = deal(varargin{:});
    nd = 0;
    msg = 'value';
    
  case 3
    [M,vf,nd] = deal(varargin{:});
    msg = 'value';
  case 5
    [M,vf1,vf2,nd1,nd2] = deal(varargin{:});
    msg = 'corr';
end

switch msg
  
  case 'value'
    out = RegGeoSum1(M,vf,nd);
    
  case 'corr'
    vf1 = vf1(:);
    Nf1 = length(vf1);
    vf2 = vf2(:);
    Nf2 = length(vf2);
    
    out = RegGeoSum1(M,-vf1(:,ones(1,Nf2))+vf2(:,ones(1,Nf1)).',nd1+nd2)*(-1)^nd1;
end

function v = RegGeoSum1(M,vf,nd)

m1 = -ceil(M/2);

ind = abs(vf)>eps^(1/3);

vf1 = vf(ind);

v = zeros(size(vf));

switch nd
  
  case 0

    v(ind) = -(exp(1i*2*pi*m1*vf1)-exp(1i*2*pi*(m1+M)*vf1))./(-1+exp(1i*2*pi*vf1));
    
    v(~ind) = M;
  
  case 1
    
    v(ind) = (-2*1i*pi*(  exp(1i*2*pi*vf1*(1 + m1))*(-1 + m1) - exp(1i*2*pi*vf1*m1)*m1  ...
	-exp(1i*pi*2*vf1*(1 + M + m1))*(-1 + M + m1) + ...
       exp(1i*2*pi*vf1*(M + m1))*(M + m1)))./(-1 + exp(1i*2*pi*vf1)).^2;
   
    v(~ind) = pi*1i*M*(-1+M+2*m1);
    
  case 2
    
    v(ind) = (4*(exp(1i*2*pi*vf1*(2 + m1))*(-1 + m1)^2 + exp(1i*2*pi*vf1*m1)*m1^2 - ...
       exp(1i*2*pi*vf1*(2 + M + m1))*(-1 + M + m1)^2 - ...
       exp(1i*2*pi*vf1*(M + m1))*(M + m1)^2 + ...
       exp(1i*2*pi*vf1*(1 + m1))*(1 - 2*(-1 + m1)*m1) + ... 
       exp(1i*2*pi*vf1*(1 + M + m1))* ...
        (-1 + 2*M^2 + 2*(-1 + m1)*m1 + M*(-2 + 4*m1)))*pi^2)./ ...
       (-1 + exp(1i*2*pi*vf1)).^3;
   
    v(~ind) = -(2*pi^2/3)*M*(1+6*m1^2+6*m1*(-1+M)-3*M+2*M^2);
    
end

function st = Report(p,vfEst)

st = [];
st.Freqs = vfEst;
st.Amps = Correlation(p.cfd.icf,'ampEstimates',vfEst);

function FreqsDevBounds = CRBBound(vf,va,M,nVar)

%Computes the Cramer-Rao bound of the frequencies.

va = va(:);
vf = vf(:);

Raa = RegGeoSum(M,vf,vf,0,0);
Rad = RegGeoSum(M,vf,vf,0,1);
Rdd = RegGeoSum(M,vf,vf,1,1);

FreqsDevBounds = sqrt((nVar/2)*diag(inv(real((conj(va)*va.') .* (Rdd-Rad'*(Raa\Rad))))));

function st = MLEstimate(y,nVar,PFA)

p = Preprocessing(y,nVar,PFA);

st = MethodSerial1(p);

st.CRBBoundsOnFreqDevs = CRBBound(st.Freqs,st.Amps,length(y),nVar);

Farrow 2-D Image Interpolation

Jesus Selva ● June 15, 2011 ● Coded in Matlab 
function Demo

%Author: J. Selva.
%Date: June 2011.

close all
format short g
format compact

T = 1; %Sampling period in both dimensions.
B = 1/(1.223*T); %This is the two-sided bandwidth in both dimensions. It is a typical ...
                 %value in SAR image co-registration.
P = 18; %Kernel semi-length.
Q = 10; %Polynomial order in Farrow structure.

%The next call to "ChebApproxPolys" computes all Farrow coefficients.
t0 = -P*T-T/2; %The parameters t0, Nt, 2*P+1 specify the 2*P+1 intervals in which Knab's ...
               %pulse is to be approximated by polynomials. 
Nt = 2*P+1;
Dt = T;
IOut = [-T/2,T/2];
  
%Perform Chebyshev interpolation. The parameters T, B and  P are passed to KnabAPPulse.
K = ChebApproxPolys('KnabAPPulse',t0,Nt,Dt,IOut,40,Q,T,B,P).';

%Evaluate the test image in a regular grid with spacing 1.
Lx = 100;
Ly = 100;

x0 = (0:Lx-1)';
x0 = x0(:,ones(1,Lx));

y0 = x0.';
x0 = x0(:);
y0 = y0(:);

disp('Computing image in a regular grid...')

tic;
Image = reshape(TestImage5(x0,y0,B),[Lx,Ly]);
disp(['(Timing: ',num2str(toc), ' sec.)'])
disp(' ')

clear x0 y0

%Specify NP points on the image at random.
NP = 10000;

x = P*T+rand(1,NP)*(Lx-2*P)*T;
y = P*T+rand(1,NP)*(Ly-2*P)*T;

%Evaluate image directly.
disp('Computing image at selected positions exactly...')
tic;
vRef = TestImage5(x,y,B);
disp(['(Timing: ',num2str(toc), ' sec.)'])
disp(' ')

%Interpolate the image using the 2D-Farrow structure. 
disp('Computing image at selected positions using Farrow 2-D...')
tic;
v = Farrow2DInterp(Image,K,x,y);
disp(['(Timing: ',num2str(toc), ' sec.)'])
disp(' ')

%Show the error. 
disp(['Maximum interpolation error: ', ...
      num2str(20*log10(max(abs(vRef-v))/max(abs(vRef)))),' dB'])

disp(' ')
disp('Display image and interpolation error:')
disp(' ')

%Evaluate the test image in a regular grid with spacing 1.
Lx = 100;
Ly = Lx;

x1 = (P*T:0.5*T:(Lx-P)*T)';
x1 = x1(:,ones(1,length(x1)));

y1 = x1.';
x1 = x1(:);
y1 = y1(:);

tic;
disp('Computing image in over-sampled grid ...')
vRef =  TestImage5(x1,y1,B);
disp(['(Timing: ',num2str(toc), ' sec.)'])
disp(' ')

tic;
disp('Computing signal in over-sampled grid using 2-D Farrow ...')
v = Farrow2DInterp(Image,K,x1,y1);
disp(['(Timing: ',num2str(toc), ' sec.)'])
disp(' ')

x1 = (P*T:0.5*T:(Lx-P)*T)';
y1 = x1.';

disp(['Maximum interpolation error: ', ...
      num2str(20*log10(max(abs(vRef-v))/max(abs(vRef)))),' dB'])

vRef = reshape(vRef,[length(x1),length(y1)]);
v = reshape(v,[length(x1),length(y1)]);

mesh(x1,y1,abs(vRef))

xlabel('x coordinate')
ylabel('y coordinate')
zlabel('Image''s abs. value')
title('Test Image (abs. value)')

figure

mesh(x1,y1,20*log10(abs(v-vRef)/max(abs(vRef(:)))))
xlabel('x coordinate')
ylabel('y coordinate')
zlabel('Interp. error (dB)')
title('Interpolation error (dB)')

function P = ChebApproxPolys(Func,t0,Nt,Dt,IOut,NPoints,Q,varargin)

%This function constructs Chebyshev interpolation polynomials in a set of consecutive ...
%intervals.
%
% Func is the function to be interpolated.
% t0 is the first instant of the first interval.
% Nt is the number of intervals.
% IOut is the interval in which the output polynomials will be evaluated.
% Dt is the interval width. The intervals are [t0,t0+Dt], [t0+Dt,t0+2*Dt],....
%    [t0+(Nt-1)*Dt,t0+Nt*Dt].
% NPoints is the number of points in which Func is evaluated, for each interval.
% Q is the output polynomial order.
% varargin contains extra parameters which are passed to varargin.

%Construct matrix that maps function values to polynomial coefficients. 
%The algorithm is that in section 9.8 of
% 
% W. H. Press et al., Numerical Recipes in C. Cambridge University
% Press, 1997.

M = MapPolyMatrix([-1,1],IOut,Q) * ChebCoefMatrix(Q) * ...
    BlockChebInvMatrix(NPoints,Q);

% Get Chebyshev roots.
r = ChebRoots(NPoints);
r = r(:);

P = zeros(Q,Nt);

%Compute polynomial for each interval.
for k = 0:Nt-1,
  v1 = feval(Func,MapIaToIb(r,[-1,1],t0+[k,k+1]*Dt),varargin{:});
  P(:,k+1) = double(M * v1(:));
end

function M = MapPolyMatrix(IIn,IOut,Q)

%For a polynomial with Q coefficients which is evaluated in interval IIn, this function ...
%gives the matrix that transforms its coefficients, so that it can be evaluated in the ...
%associated points in interval IOut. 

ayx = (IIn(1)-IIn(2))/(IOut(1)-IOut(2));
byx = IIn(1)-ayx*IOut(1);

if byx == 0
  M = diag(ayx.^(0:Q-1));
else
  M = zeros(Q);
  for q = 0:Q-1
    for r = 0:q
      M(r+1,q+1) = nchoosek(q,r) * ayx^r * byx^(q-r);
    end
  end
end

function M = ChebCoefMatrix(N)

M = zeros(N);

M(1) = 1;

if N > 1
  M(2,2) = 1;
end

for n = 3:N,
  M(2:n,n) = 2*M(1:n-1,n-1);
  M(1:n-2,n) = M(1:n-2,n) - M(1:n-2,n-2);
end

function M = BlockChebInvMatrix(N,Q)

r = ChebRoots(N);
r= r(:);

M = zeros(N,Q);

if Q >=1
  M(:,1) = 1;
end

if Q >= 2
  M(:,2) = r;
end

for q = 3:Q,
  M(:,q) = 2 * r .* M(:,q-1) - M(:,q-2);
end

M = M.'/N;
M(2:end,:) = M(2:end,:) * 2;

function R = ChebRoots(N)

R = cos(pi*((1:N)-1/2)/N);

function y = MapIaToIb(x,Ia,Ib)

%Linear transformation that maps interval Ia onto interval Ib.

y = (Ib(1)+Ib(2))/2 + (x-(Ia(1)+Ia(2))/2) * ((Ib(2)-Ib(1))/(Ia(2)-Ia(1))); 

function [g,BoundAP] = KnabAPPulse(t,T,B,P)

%Knab interpolation pulse in 
%
%  J. J. Knab, 'Interpolation of band-limited functions using the Approximate
%  Prolate series', IEEE Transactions on Information Theory, vol. IT-25, 
%  no. 6, pp. 717-720, Nov 1979. 
%
%Grid points -P*T:T:P*T.
%Interpolation interval: [-T/2,T/2].

g = sinc(t/T) .* sinc((1/T-B)*sqrt(t.^2-(P*T)^2)) / sinc(j*(1-B*T)*P);

if nargout == 2
  BoundAP = 1/sinh(P*pi*(1-B*T));
end

function v = TestImage5(x,y,B)

%Test image formed by a sum of random sincs.

Seed = rand('state');

rand('state',362436069);

L = 100; %Image size in both dimensions.

Np = 5000;
p = rand(Np,2)*L;
A = exp(j*2*pi*rand(Np,1));

rand('state',Seed);

v = zeros(size(x));

for k = 1:length(x(:))
  v(k) = A(:).' * (sinc(B*(x(k)-p(:,1))) .* sinc(B*(y(k)-p(:,2))));
end

function vx = Farrow2DInterp(Image,Kernel,x,y)

%Farrow interpolation in two dimensions using interlaced FFTs.

vxSize = size(x);

[LK,NK] = size(Kernel);
[Lx,Ly] = size(Image);

nx = round(x);
ny = round(y);
ux = x-nx;
uy = y-ny;

np = 1+nx+ny*Lx;
np = np(:);

clear x y nx ny

P = (LK-1)/2;

NFFTx = 2^ceil(log2(Lx+LK-1));
NFFTy = 2^ceil(log2(Ly+LK-1));

KernelFx = fft(Kernel,NFFTx);
KernelFy = fft(Kernel,NFFTy);

ImageFx = fft(Image,NFFTx);

vx = zeros(size(np))+j*zeros(size(np));
for kx = NK:-1:1,
 
  ImageCx = ifft(ImageFx .* KernelFx(:,kx(ones(1,Ly))),NFFTx); 
  ImageCx = ImageCx(1+P:Lx+LK-1-P,:);
  
  ImageCxFy = fft(ImageCx,NFFTy,2);
  
  vy = zeros(size(np))+j*zeros(size(np));
  for ky = NK:-1:1,
    ImageCxCy = ifft(ImageCxFy .* KernelFy(:,ky(ones(1,Lx))).',NFFTy,2);
    ImageCxCy = ImageCxCy(:,1+P:Ly+LK-1-P);
  
    vy(:) = vy(:) .* uy(:) + ImageCxCy(np);
  end
  
  vx(:) = vx(:) .* ux(:) + vy(:);
end

vx = reshape(vx,vxSize);

Interpolation from Nonuniform Samples

Jesus Selva ● May 18, 2011●1 comment ● Coded in Matlab 
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

ML frequency estimation using the FFT

Jesus Selva ● April 18, 2011 ● Coded in Matlab 
function Demo(L)

%Author: J. Selva. 
%Date: April 2011.

%This demo constructs trials of an undamped exponential plus noise, and then computes the
%maximum likelihood estimate of its frequency, amplitude, and phase. Only one FFT of size
%2^ceil(log2(st.M)+1) is computed, plus a few operations. The complexity is that given by
%the FFT. The algoritm is explained in 
%
%         http://arxiv.org/abs/1104.3069
%

format long g
format compact

M = 1024; %Number of samples.
gdB = 20; %SNR at matched filter output.
ARef = 1.3; %Amplitude.
fRef = 0.372; %Frequency in [-0.5,0.5[.
PhiRef = 0.43; %Phase (rad).

if nargin > 0
  for k = 1:2:length(L)
    eval([L{k} '=L{k+1};']);
  end
end

%Define a random vector with 1024 elements. 

SNRdB =  gdB - 10*log10(M); %Signal-to-noise ratio (dB).

%Complex sinusoidal amplitude, frequency, and phase.

NVar = ARef^2*10.^(-SNRdB/10); %Noise variance. 

m1 = -ceil(M/2);

while 1
  
  disp(' ')
  disp(['M = ', num2str(M), ', gdB = ', num2str(gdB), ' dB, ARef = ', num2str(ARef), ...
	', fRef = ', num2str(fRef), ', PhiRef = ', num2str(PhiRef), ' rad'])
  disp(' ')
  
  %Set up a trial.
  a = ARef*exp(1i*2*pi*fRef*(m1:m1+M-1).'+1i*PhiRef)+(randn(M,1)+1i*randn(M,1))*sqrt(NVar/2);

  %Call the nufft function. The output is a struct with some variables for later use.

  st = nufft(a);

  %Compute the DFT peak. 
  [f1,A1,Phi1,L1,Deg,nIt] = MLEstimate(st);

  %The output is the following, 
  %
  %   f1, A1, and Phi1 are the ML estimates of frequency, amplitude and phase, respectively. 
  %
  %   L1 is the ML cost function value at f1 and its first two  derivatives. L1(1) should be
  %     large   since it is   the square value of   the peak amplitude.    L1(2) is the cost
  %     function   derivative,  and  it  should  be   small  since  we  are  at   a critical
  %     point. Finally  L1(2) should be a  negative number with large  amplitude since f1 is
  %     the abcissa of a maximum (second derivative).
  %
  %   Deg: If Deg==1 then the Newton iteration was cut too many times. 
  %   nIt: Number of Newton iterations. 

  disp('Exact and estimated frequency, amplitude, and phase:')
  disp([fRef,ARef,PhiRef;f1,A1,Phi1])
  disp(['Number of iterations: ' num2str(nIt)])
  %
  %Plot the spectrum close to the ML estimate
  %
  f = fRef+(-10:.05:10)/st.K;
  sp = nufft(st,f);
  plot(f,20*log10(abs(sp)/M))

  hold on
  stem(f1,20*log10(A1),'r')
  hold off

  grid on
  set(gca,'YLim',[-40,10])
  xlabel('f')
  ylabel('Amplitude (dB)')
  text(fRef+2/st.K,5,{['Freq. ' num2str(f1)],['Ampl. ' num2str(A1)],...
      ['Phase. (rad) ' num2str(Phi1)]})
  %

  R = input('Press Enter for a new trial, q to quit: ','s');

  if any(strcmp(R,{'q','Q'}))
    break
  end

end

%------------------------------------------------------
%------------------------------------------------------

function [f1,A1,Phi1,L1,Deg,nIt] = MLEstimate(st)

%This function finds out the peak of a DFT. Its input is the struct produced by the nufft function. 

%Find the spectral peak in the frequency grid.
[pr11,ind] = max(abs(st.aF).^2);

%Compute its frequency.
f0 = (ind-1)/st.K;
f0 = f0+ceil(-0.5-f0);

%L0 contains the differentials of the ML cost function of orders 0,1, 2 ...
%at f0. 
L0 = LDiff(nufft(st,f0,2));

nItMax = 20;
err = 10.^-12;
nIt = 0;

%Call Iter (Newton iteration) until there is no change or the solution is degenerate.
while 1
  [Deg,pr11,f1,L1] = Iter(st,f0,L0);
  
  if Deg || abs(f1-f0) <= err || nIt >= nItMax
    break
  end
  
  f0 = f1;
  L0 = L1;
  nIt = nIt+1;
end

pr = nufft(st,f1);
A1 = abs(pr)/st.M;
Phi1 = angle(pr); 

function L = LDiff(c)

%Computes the differentials of the ML cost function from the nufft differentials.

L = [real(c(1)*conj(c(1))); 2*real(conj(c(1))*c(2)); ...
  2*real(conj(c(2))*c(2)+conj(c(1))*c(3))];

function [Deg,nCut,f1,L1] = Iter(st,f0,L0)

%Performs the Newton iteration, cutting the Newton direction if necessary.

f1 = f0-L0(2)/L0(3); %Newton iteration.
f1 = f1+ceil(-0.5-f1); %Make it fall in [-0.5,0.5].
L1 = LDiff(nufft(st,f1,2)); %Compute the differentials at f1.

nCutMax = 20; %The difference f1-f0 will be reduced in the sequel if there is no ...
              %improvement in the cost function
nCut = 0;
d = 1;

%In this loop the difference f1-f0 is halved a maximum of nCutMax times until there is ...
%an increase in the cost function

while L1(1) < L0(1) && nCut < nCutMax
  d = d/2;
  f1 = f0-d*L0(2)/L0(3);
  f1 = f1+ceil(-0.5-f1);
  L1 = LDiff(nufft(st,f1,2)); 
  nCut = nCut+1;
end

Deg = nCut >= nCutMax;

%----------------------------------------------------------------
%----------------------------------------------------------------

%Non-uniform FFT function. This code is discussed in 
%   http://www.dsprelated.com/showcode/159.php

function varargout = nufft(varargin)

%Author: J. Selva
%Date: April 2011.
%
%nufft computes the FFT at arbitrary frequencies using barycentric interpolation.
%
%For the interpolation method, 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.
%
%Usage:
%
% -First call nufft with the vector of samples as argument,
%
%       st = nufft(a); %a is the vector of samples.
%
%    The output is an struct. The field st.aF is the DFT of vector a, sampled with spacing
%    1/st.K. 
%
%    The DFT is defined using
%
%        A_k = \sum_{n=m_1}^{m_1+M-1} a_n e^{-j 2 \pi n f}
%
%    where m_1=-ceil(M/2) and M the number of samples. 
%
% -Then call nufft with sintax 
%
%        out = nufft(st,f,nd,method)
%
%   where
%
%     f: Vector of frequencies where the DFT is evaluated. Its elements must follow
%        abs(f(k))<=0.5
%
%     nd: Derivative order. nufft computes derivatives up to order nd. At the output 
%         out(p,q) is the (q-1)-order derivative of the DFT at frequency f(p). The first
%         column of out contains the zero-order derivatives, i.e, the values of the DFT at
%         frequencies in vector f. 
%     method: If omitted, it is 'baryVec'. Four methods have been implemented:
%
%         +'baryLoop': The output is interpolated using the barycentric method and a loop
%           implementation.
%
%         +'baryVec': The output is interpolated using the barycentric method and a
%          vectorized implementation.
%
%         +'directLoop': The output is computed using the DFT sum directly and a loop
%           implementation.
%  
%         +'directVec': The output is computed using the DFT sum directly and a vectorized
%           implementation.

 
if ~isstruct(varargin{1})
  st = [];
  st.a = varargin{1};
  st.a = st.a(:);
  
  st.M = length(st.a);
  st.m1 = -ceil(st.M/2);
  st.K =  2^ceil(log2(st.M)+1);
  st.aF = fft([st.a(-st.m1+1:st.M); zeros(st.K-st.M,1); st.a(1:-st.m1)]);
  
  errb = 10.^-13; %This is the interpolation accuracy. You can set any larger value, and ...
      %this would reduce st.P. The number of interpolation samples is 2*st.P+1.
  
  st.T = 1/st.K;
  st.B = -2*st.m1;
  st.P = ceil(acsch(errb)/(pi*(1-st.B*st.T)));
  st.vt = MidElementToEnd((-st.P:st.P)*st.T);
  
  st.alpha =  MidElementToEnd(baryWeights(st.T,st.B,st.P));
  
  varargout = {st};
  
  return
end

[st,f] = deal(varargin{1:2});
  
nd = 0;
   
if nargin > 2
  nd = varargin{3};
end
  
method = 'baryVec';
  
if nargin > 3
  method = varargin{4};
end
  
Nf = length(f);
out = zeros(Nf,nd+1);

switch method
    
  case 'baryLoop' %Interpolated method using loops

    for k = 1:length(f)

      x = f(k);
      
      n = floor(x/st.T+0.5);
      u = x - n * st.T;

      da = MidElementToEnd(st.aF(1+mod(n-st.P:n+st.P,st.K)).');

      out(k,:) = DerBarycentricInterp3(st.alpha,da,st.vt,u,nd);

    end

  case 'baryVec' %Vectorized interpolated method
    
    f = f(:);
    Nf = length(f);
    
    n = floor(f/st.T+0.5);
    u = f - n * st.T;
    
    pr = [-st.P:-1 , 1:st.P , 0];
    
    ms = st.aF(1+mod(n(:,ones(1,2*st.P+1)) + pr(ones(Nf,1),:),st.K));
    
    if length(f) == 1
      ms = ms.';
    end
    
    out = DerBarycentricInterp3Vec(st.alpha,ms,st.vt,u,nd);
    
  case 'directLoop' % Direct method using loops
      
    for k = 1:length(f)
      out(k,1) = 0;
	
      for r = st.m1:st.m1+st.M-1
	out(k,1) = out(k,1)+exp(-1i*2*pi*r*f(k))*st.a(r-st.m1+1);
      end
	
      for kd = 1:nd
	out(k,kd+1) = 0;
	
	for r = st.m1:st.m1+st.M-1
	  out(k,kd+1) = out(k,kd+1) + ...
	      ((-1i*2*pi*r).^kd .* exp(-1i*2*pi*r*f(k)))*st.a(r-st.m1+1);
	end
      
      end
    end  
    
  case 'directVec' %Vectorized direct method
	
    for k = 1:length(f)
      out(k,1) = exp(-1i*2*pi*(st.m1:st.m1+st.M-1)*f(k))*st.a;
      
      for kd = 1:nd
	out(k,kd+1) = ...
	    ((-1i*2*pi*(st.m1:st.m1+st.M-1)).^kd .* ...
	    exp(-1i*2*pi*(st.m1:st.m1+st.M-1)*f(k)))*st.a;
      end
      
    end
   
end
 
varargout = {out};

function v = MidElementToEnd(v)

ind = ceil(length(v)/2);
v = [v(1:ind-1),v(ind+1:end),v(ind)];

function v = APPulse(t,B,TSL)

v = real(sinc(B*sqrt(t.^2-TSL^2)))/real(sinc(1i*pi*B*TSL));

function w = baryWeights(T,B,P)

vt = (-P:P)*T;
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 = DerBarycentricInterp3(alpha,s,t,tau,nd)

vD = 0;
Nt = length(t);
LF = [1,zeros(1,Nt-1)];
out = zeros(1,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;

for kd = 0:nd

  z1 = z-z(end); cN = 0;

  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(1:end-1)-tau) , cN ] * (kd+1);
  end
  
end

function out = DerBarycentricInterp3Vec(alpha,zL,t,tauL,nd)

NtauL = size(tauL,1);
LBtau = 200; %Block size for vectorized processing. Change this to adjust the memory 
             %usage. 

NBlock = 1+floor(NtauL/LBtau);

Nt = length(t);

out = zeros(NtauL,nd+1);

for r = 0:NBlock-1
  ind1 = 1+r*LBtau;
  ind2 = min([(r+1)*LBtau,NtauL]);
  
  Ntau = ind2-ind1+1;
  z = zL(ind1:ind2,:);
  tau = tauL(ind1:ind2);
  
  vD = zeros(Ntau,1);
  
  LF = [1,zeros(1,Nt-1)];
  LF = LF(ones(Ntau,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);

  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(ind1:ind2,kd+1) = ztau;
  
    if kd < nd
      pr1 = ztau;
      pr1 = z(:,1:end-1) - pr1(:,ones(1,Nt-1));
      
      pr2 = t(1:end-1);
      pr2 = pr2(ones(Ntau,1),:)-tau(:,ones(1,Nt-1));
      
      z = [pr1 ./ pr2, cN] * (kd+1);
      
    end
  end
  
end

FFT at arbitrary frequencies (non-uniform FFT)

Jesus Selva ● April 4, 2011●3 comments ● Coded in Matlab 
Following you find the code of a demo script (Demo.m) and that of a function (nufft.m). Just place them in separate m files (Demo.m and nufft.m) and execute the first. 

------------------------------------------------------
Code of file Demo.m
------------------------------------------------------

echo on 

%Author: J. Selva. 
%Date: April 2011.

%Define a random vector with 1024 elements. 

M = 1024;
a = rand(M,1)+1i*rand(M,1);

%Call the nufft function. The output is a struct with some variable for later use.

st = nufft(a)

%The field st.aF is the fft of size st.K.

%Define a vector of random frequencies.

f = rand([2*M,1])-0.5;

%Compute the DFT at frequencies in f and the derivatives of first and second order using
%loops and direct evaluation. This takes a while
%

tic; outDL = nufft(st,f,2,'directLoop'); tDL = toc;

%The function returns the DFT values in the first column, the first derivatives in the
%second column, and so on.

outDL(1:10,:)

%Do the same computation using interpolation. This is fast.
%
tic; outIL = nufft(st,f,2); tIL = toc;

%Compare the execution times
%

[tDL,tIL]

%The executions are faster if the code is vectorized. This is an example
%
tic; outDV = nufft(st,f,2,'directVec'); tDV = toc;
tic; outIV = nufft(st,f,2,'baryVec'); tIV = toc;

[tDV, tIV]

%For larger data sizes the difference is more significant. To  check this, let us repeat
%the same experiment with M=1024*5.
%

M = 1024*5;
a = rand(M,1)+1i*rand(M,1);
st = nufft(a);

f = rand([2*M,1])-0.5;

tic; outDV = nufft(st,f,2,'directVec'); tDV = toc; %This line takes a while.
tic; outIV = nufft(st,f,2); tIV = toc;

%These are the timings:

[tDV,tIV]

%Finally, let  us check the accuracy of the interpolated values

max(abs(outDV-outIV))./max(abs(outDV))

echo off

----------------------------------------------------------
----------------------------------------------------------

Code of file nufft.m. 

----------------------------------------------------------
----------------------------------------------------------

function varargout = nufft(varargin)

%Author: J. Selva
%Date: April 2011.
%
%nufft computes the FFT at arbitrary frequencies using barycentric interpolation.
%
%For the interpolation method, 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.
%
%Usage:
%
% -First call nufft with the vector of samples as argument,
%
%       st = nufft(a); %a is the vector of samples.
%
%    The output is an struct. The field st.aF is the DFT of vector a, sampled with spacing
%    1/st.K. 
%
%    The DFT is defined using
%
%        A_k = \sum_{n=m_1}^{m_1+M-1} a_n e^{-j 2 \pi n f}
%
%    where m_1=-ceil(M/2) and M the number of samples. 
%
% -Then call nufft with sintax 
%
%        out = nufft(st,f,nd,method)
%
%   where
%
%     f: Vector of frequencies where the DFT is evaluated. Its elements must follow
%        abs(f(k))<=0.5
%
%     nd: Derivative order. nufft computes derivatives up to order nd. At the output 
%         out(p,q) is the (q-1)-order derivative of the DFT at frequency f(p). The first
%         column of out contains the zero-order derivatives, i.e, the values of the DFT at
%         frequencies in vector f. 
%     method: If omitted, it is 'baryVec'. Four methods have been implemented:
%
%         +'baryLoop': The output is interpolated using the barycentric method and a loop
%           implementation.
%
%         +'baryVec': The output is interpolated using the barycentric method and a
%          vectorized implementation.
%
%         +'directLoop': The output is computed using the DFT sum directly and a loop
%           implementation.
%  
%         +'directVec': The output is computed using the DFT sum directly and a vectorized
%           implementation.

 
if ~isstruct(varargin{1})
  st = [];
  st.a = varargin{1};
  st.a = st.a(:);
  
  st.M = length(st.a);
  st.m1 = -ceil(st.M/2);
  st.K =  2^ceil(log2(st.M)+1);
  st.aF = fft([st.a(-st.m1+1:st.M); zeros(st.K-st.M,1); st.a(1:-st.m1)]);
  
  errb = 10.^-13; %This is the interpolation accuracy. You can set any larger value, and ...
      %this would reduce st.P. The number of interpolation samples is 2*st.P+1.
  
  st.T = 1/st.K;
  st.B = -2*st.m1;
  st.P = ceil(acsch(errb)/(pi*(1-st.B*st.T)));
  st.vt = MidElementToEnd((-st.P:st.P)*st.T);
  
  st.alpha =  MidElementToEnd(baryWeights(st.T,st.B,st.P));
  
  varargout = {st};
  
  return
end

[st,f] = deal(varargin{1:2});
  
nd = 0;
   
if nargin > 2
  nd = varargin{3};
end
  
method = 'baryVec';
  
if nargin > 3
  method = varargin{4};
end
  
Nf = length(f);
out = zeros(Nf,nd+1);

switch method
    
  case 'baryLoop' %Interpolated method using loops

    for k = 1:length(f)

      x = f(k);
      
      n = floor(x/st.T+0.5);
      u = x - n * st.T;

      da = MidElementToEnd(st.aF(1+mod(n-st.P:n+st.P,st.K)).');

      out(k,:) = DerBarycentricInterp3(st.alpha,da,st.vt,u,nd);

    end

  case 'baryVec' %Vectorized interpolated method
    
    f = f(:);
    Nf = length(f);
    
    n = floor(f/st.T+0.5);
    u = f - n * st.T;
    
    pr = [-st.P:-1 , 1:st.P , 0];
    
    ms = st.aF(1+mod(n(:,ones(1,2*st.P+1)) + pr(ones(Nf,1),:),st.K));
    
    if length(f) == 1
      ms = ms.';
    end
    
    out = DerBarycentricInterp3Vec(st.alpha,ms,st.vt,u,nd);
    
  case 'directLoop' % Direct method using loops
      
    for k = 1:length(f)
      out(k,1) = 0;
	
      for r = st.m1:st.m1+st.M-1
	out(k,1) = out(k,1)+exp(-1i*2*pi*r*f(k))*st.a(r-st.m1+1);
      end
	
      for kd = 1:nd
	out(k,kd+1) = 0;
	
	for r = st.m1:st.m1+st.M-1
	  out(k,kd+1) = out(k,kd+1) + ...
	      ((-1i*2*pi*r).^kd .* exp(-1i*2*pi*r*f(k)))*st.a(r-st.m1+1);
	end
      
      end
    end  
    
  case 'directVec' %Vectorized direct method
	
    for k = 1:length(f)
      out(k,1) = exp(-1i*2*pi*(st.m1:st.m1+st.M-1)*f(k))*st.a;
      
      for kd = 1:nd
	out(k,kd+1) = ...
	    ((-1i*2*pi*(st.m1:st.m1+st.M-1)).^kd .* ...
	    exp(-1i*2*pi*(st.m1:st.m1+st.M-1)*f(k)))*st.a;
      end
      
    end
      
end
 
varargout = {out};

function v = MidElementToEnd(v)

ind = ceil(length(v)/2);
v = [v(1:ind-1),v(ind+1:end),v(ind)];

function v = APPulse(t,B,TSL)

v = real(sinc(B*sqrt(t.^2-TSL^2)))/real(sinc(1i*pi*B*TSL));

function w = baryWeights(T,B,P)

vt = (-P:P)*T;
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 = DerBarycentricInterp3(alpha,s,t,tau,nd)

vD = 0;
Nt = length(t);
LF = [1,zeros(1,Nt-1)];
out = zeros(1,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;

for kd = 0:nd

  z1 = z-z(end); cN = 0;

  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(1:end-1)-tau) , cN ] * (kd+1);
  end
  
end

function out = DerBarycentricInterp3Vec(alpha,zL,t,tauL,nd)

NtauL = size(tauL,1);
LBtau = 200; %Block size for vectorized processing. Change this to adjust the memory 
             %usage. 

NBlock = 1+floor(NtauL/LBtau);

Nt = length(t);

out = zeros(NtauL,nd+1);

for r = 0:NBlock-1
  ind1 = 1+r*LBtau;
  ind2 = min([(r+1)*LBtau,NtauL]);
  
  Ntau = ind2-ind1+1;
  z = zL(ind1:ind2,:);
  tau = tauL(ind1:ind2);
  
  vD = zeros(Ntau,1);
  
  LF = [1,zeros(1,Nt-1)];
  LF = LF(ones(Ntau,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);

  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(ind1:ind2,kd+1) = ztau;
  
    if kd < nd
      pr1 = ztau;
      pr1 = z(:,1:end-1) - pr1(:,ones(1,Nt-1));
      
      pr2 = t(1:end-1);
      pr2 = pr2(ones(Ntau,1),:)-tau(:,ones(1,Nt-1));
      
      z = [pr1 ./ pr2, cN] * (kd+1);
      
    end
  end
  
end

Efficient update of interpolation FIR coefficients

Jesus Selva ● November 10, 2010 ● Coded in Matlab 
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));

High accuracy FFT interpolation

Jesus Selva ● November 4, 2010●3 comments ● Coded in Matlab 
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);