## Interleave Zeros

July 10, 20111 comment Coded in Matlab
function out = interleave_zeros(a)

% Interleave vector with zeros.
% Only 1D matrices are supported.
%
% Usage: b = INTERLEAVE_ZEROS(a);
%
%        A is the vector that will have zeros inserted
%        B will contain the vector set A and the interleaved
%          zeros (twice the length)
%
% Author: sparafucile17 2/13/04

my_zeros = zeros(1, length(a));
out = matintrlv([a  my_zeros],2,length(a));

## Zero Crossing Counter

July 9, 20119 comments Coded in Matlab
function count = zero_crossings(x)
% Count the number of zero-crossings from the supplied time-domain
% input vector.  A simple method is applied here that can be easily
% ported to a real-time system that would minimize the number of
% if-else conditionals.
%
% Usage:     COUNT = zero_crossings(X);
%
%            X     is the input time-domain signal (one dimensional)
%            COUNT is the amount of zero-crossings in the input signal
%
% Author:    sparafucile17 06/27/04
%

% initial value
count = 0;

% error checks
if(length(x) == 1)
error('ERROR: input signal must have more than one element');
end

if((size(x, 2) ~= 1) && (size(x, 1) ~= 1))
error('ERROR: Input must be one-dimensional');
end

% force signal to be a vector oriented in the same direction
x = x(:);

num_samples = length(x);
for i=2:num_samples

% Any time you multiply to adjacent values that have a sign difference
% the result will always be negative.  When the signs are identical,
% the product will always be positive.
if((x(i) * x(i-1)) < 0)
count = count + 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

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

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

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

MRsy = pr(:,1);
MRsd = pr(:,2:end);

varargout = {st.EnergyRyy-RealTrace(Rsy',MRsy),-2*RealDiag(MRsy,Rdy'-MRsy'*Rsd),...

case 'value0'

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

vf = varargin{3};
Rsy = Correlation(st.icf,'value',vf);

end

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

function C = RealDiag(A,B)

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

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=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);
Rdd = RegGeoSum(M,vf,vf,1,1);

function st = MLEstimate(y,nVar,PFA)

p = Preprocessing(y,nVar,PFA);

st = MethodSerial1(p);

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

## Digital Comb Filter

June 21, 2011 Coded in Matlab
% Usage:     [B,A] = COMB(order, scalar);
%
%             ORDER is the number of samples delayed prior to add
%             SCALAR is the coefficient that will be applied to
%                the delayed signal path at the final summation block.
%
% Note, there are two types of comb filters.  A DC-blocker and a DC-passer.
% To get a DC-Blocker (tooth at DC), pass in a -1 for the scalar.
% To get a DC-Passer (+6dB at DC), pass in a +1 for the scalar.
%
% By default, if the scalar is not passed, a DC-Passer is assumed.
%
% Author: sparafucile17 03/16/04

% Validate that the proper argument count was supplied
error(nargchk(1, 2, nargin));

% Use scalar is passed as an argument, otherwise assume scalar=1;
if (nargin == 1)
scalar = 1;
else
scalar = varargin{1};
end

% Input has zeros to simulate a single delay of N samples
a = [ 1 ];
b = [ 1 zeros(1, N-1) scalar*1];

## Equal Loudness Curves (ISO226)

June 21, 20111 comment Coded in Matlab
function [spl, freq] = iso226(phon);
%
% Generates an Equal Loudness Contour as described in ISO 226
%
% Usage:  [SPL FREQ] = ISO226(PHON);
%
%         PHON is the phon value in dB SPL that you want the equal
%           loudness curve to represent. (1phon = 1dB @ 1kHz)
%         SPL is the Sound Pressure Level amplitude returned for
%           each of the 29 frequencies evaluated by ISO226.
%         FREQ is the returned vector of frequencies that ISO226
%           evaluates to generate the contour.
%
% Desc:   This function will return the equal loudness contour for
%         your desired phon level.  The frequencies evaulated in this
%         function only span from 20Hz - 12.5kHz, and only 29 selective
%         frequencies are covered.  This is the limitation of the ISO
%         standard.
%
%         In addition the valid phon range should be 0 - 90 dB SPL.
%         Values outside this range do not have experimental values
%         and their contours should be treated as inaccurate.
%
%         If more samples are required you should be able to easily
%         interpolate these values using spline().
%
% Author: sparafucile17 03/01/05

%                /---------------------------------------\
%%%%%%%%%%%%%%%%%          TABLES FROM ISO226             %%%%%%%%%%%%%%%%%
%                \---------------------------------------/
f = [20 25 31.5 40 50 63 80 100 125 160 200 250 315 400 500 630 800 ...
1000 1250 1600 2000 2500 3150 4000 5000 6300 8000 10000 12500];

af = [0.532 0.506 0.480 0.455 0.432 0.409 0.387 0.367 0.349 0.330 0.315 ...
0.301 0.288 0.276 0.267 0.259 0.253 0.250 0.246 0.244 0.243 0.243 ...
0.243 0.242 0.242 0.245 0.254 0.271 0.301];

Lu = [-31.6 -27.2 -23.0 -19.1 -15.9 -13.0 -10.3 -8.1 -6.2 -4.5 -3.1 ...
-2.0  -1.1  -0.4   0.0   0.3   0.5   0.0 -2.7 -4.1 -1.0  1.7 ...
2.5   1.2  -2.1  -7.1 -11.2 -10.7  -3.1];

Tf = [ 78.5  68.7  59.5  51.1  44.0  37.5  31.5  26.5  22.1  17.9  14.4 ...
11.4   8.6   6.2   4.4   3.0   2.2   2.4   3.5   1.7  -1.3  -4.2 ...
-6.0  -5.4  -1.5   6.0  12.6  13.9  12.3];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Error Trapping
if((phon < 0) | (phon > 90))
disp('Phon value out of bounds!')
spl = 0;
freq = 0;
else
%Setup user-defined values for equation
Ln = phon;

%Deriving sound pressure level from loudness level (iso226 sect 4.1)
Af=4.47E-3 * (10.^(0.025*Ln) - 1.15) + (0.4*10.^(((Tf+Lu)/10)-9 )).^af;
Lp=((10./af).*log10(Af)) - Lu + 94;

%Return user data
spl = Lp;
freq = f;
end

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

## Shelving Filter Design

May 29, 20111 comment Coded in Matlab
function [b, a]  = shelving(G, fc, fs, Q, type)

%
% Derive coefficients for a shelving filter with a given amplitude and
% cutoff frequency.  All coefficients are calculated as described in
% Zolzer's DAFX book (p. 50 -55).
%
% Usage:     [B,A] = shelving(G, Fc, Fs, Q, type);
%
%            G is the logrithmic gain (in dB)
%            FC is the center frequency
%            Fs is the sampling rate
%            Q adjusts the slope be replacing the sqrt(2) term
%            type is a character string defining filter type
%                 Choices are: 'Base_Shelf' or 'Treble_Shelf'
%
% Author:    sparafucile17 08/22/05
%

%Error Check
if((strcmp(type,'Base_Shelf') ~= 1) && (strcmp(type,'Treble_Shelf') ~= 1))
error(['Unsupported Filter Type: ' type]);
end

K = tan((pi * fc)/fs);
V0 = 10^(G/20);
root2 = 1/Q; %sqrt(2)

%Invert gain if a cut
if(V0 < 1)
V0 = 1/V0;
end

%%%%%%%%%%%%%%%%%%%%
%    BASE BOOST
%%%%%%%%%%%%%%%%%%%%
if(( G > 0 ) & (strcmp(type,'Base_Shelf')))

b0 = (1 + sqrt(V0)*root2*K + V0*K^2) / (1 + root2*K + K^2);
b1 =             (2 * (V0*K^2 - 1) ) / (1 + root2*K + K^2);
b2 = (1 - sqrt(V0)*root2*K + V0*K^2) / (1 + root2*K + K^2);
a1 =                (2 * (K^2 - 1) ) / (1 + root2*K + K^2);
a2 =             (1 - root2*K + K^2) / (1 + root2*K + K^2);

%%%%%%%%%%%%%%%%%%%%
%    BASE CUT
%%%%%%%%%%%%%%%%%%%%
elseif (( G < 0 ) & (strcmp(type,'Base_Shelf')))

b0 =             (1 + root2*K + K^2) / (1 + root2*sqrt(V0)*K + V0*K^2);
b1 =                (2 * (K^2 - 1) ) / (1 + root2*sqrt(V0)*K + V0*K^2);
b2 =             (1 - root2*K + K^2) / (1 + root2*sqrt(V0)*K + V0*K^2);
a1 =             (2 * (V0*K^2 - 1) ) / (1 + root2*sqrt(V0)*K + V0*K^2);
a2 = (1 - root2*sqrt(V0)*K + V0*K^2) / (1 + root2*sqrt(V0)*K + V0*K^2);

%%%%%%%%%%%%%%%%%%%%
%   TREBLE BOOST
%%%%%%%%%%%%%%%%%%%%
elseif (( G > 0 ) & (strcmp(type,'Treble_Shelf')))

b0 = (V0 + root2*sqrt(V0)*K + K^2) / (1 + root2*K + K^2);
b1 =             (2 * (K^2 - V0) ) / (1 + root2*K + K^2);
b2 = (V0 - root2*sqrt(V0)*K + K^2) / (1 + root2*K + K^2);
a1 =              (2 * (K^2 - 1) ) / (1 + root2*K + K^2);
a2 =           (1 - root2*K + K^2) / (1 + root2*K + K^2);

%%%%%%%%%%%%%%%%%%%%
%   TREBLE CUT
%%%%%%%%%%%%%%%%%%%%

elseif (( G < 0 ) & (strcmp(type,'Treble_Shelf')))

b0 =               (1 + root2*K + K^2) / (V0 + root2*sqrt(V0)*K + K^2);
b1 =                  (2 * (K^2 - 1) ) / (V0 + root2*sqrt(V0)*K + K^2);
b2 =               (1 - root2*K + K^2) / (V0 + root2*sqrt(V0)*K + K^2);
a1 =             (2 * ((K^2)/V0 - 1) ) / (1 + root2/sqrt(V0)*K + (K^2)/V0);
a2 = (1 - root2/sqrt(V0)*K + (K^2)/V0) / (1 + root2/sqrt(V0)*K + (K^2)/V0);

%%%%%%%%%%%%%%%%%%%%
%   All-Pass
%%%%%%%%%%%%%%%%%%%%
else
b0 = V0;
b1 = 0;
b2 = 0;
a1 = 0;
a2 = 0;
end

%return values
a = [  1, a1, a2];
b = [ b0, b1, b2];

## Peak/Notch Filter Design

May 21, 20113 comments Coded in Matlab
function [b, a]  = peaking(G, fc, Q, fs)

% Derive coefficients for a peaking filter with a given amplitude and
% bandwidth.  All coefficients are calculated as described in Zolzer's
% DAFX book (p. 50 - 55).  This algorithm assumes a constant Q-term
% is used through the equation.
%
% Usage:     [B,A] = peaking(G, Fc, Q, Fs);
%
%            G is the logrithmic gain (in dB)
%            FC is the center frequency
%            Q is Q-term equating to (Fb / Fc)
%            Fs is the sampling rate
%
% Author:    sparafucile17 08/22/05
%

K = tan((pi * fc)/fs);
V0 = 10^(G/20);

%Invert gain if a cut
if(V0 < 1)
V0 = 1/V0;
end

%%%%%%%%%%%%%%
%   BOOST
%%%%%%%%%%%%%%
if( G > 0 )

b0 = (1 + ((V0/Q)*K) + K^2) / (1 + ((1/Q)*K) + K^2);
b1 =        (2 * (K^2 - 1)) / (1 + ((1/Q)*K) + K^2);
b2 = (1 - ((V0/Q)*K) + K^2) / (1 + ((1/Q)*K) + K^2);
a1 = b1;
a2 =  (1 - ((1/Q)*K) + K^2) / (1 + ((1/Q)*K) + K^2);

%%%%%%%%%%%%%%
%    CUT
%%%%%%%%%%%%%%
else

b0 = (1 + ((1/Q)*K) + K^2) / (1 + ((V0/Q)*K) + K^2);
b1 =       (2 * (K^2 - 1)) / (1 + ((V0/Q)*K) + K^2);
b2 = (1 - ((1/Q)*K) + K^2) / (1 + ((V0/Q)*K) + K^2);
a1 = b1;
a2 = (1 - ((V0/Q)*K) + K^2) / (1 + ((V0/Q)*K) + K^2);

end

%return values
a = [  1, a1, a2];
b = [ b0, b1, b2];

## 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[.

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)],...
%

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