## How to speed up the sinc series

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)

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

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

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

May 18, 20111 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

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)

April 4, 20113 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

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

November 4, 20103 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);