Jesus Selva ● July 8, 2011
● Coded in
Matlab function st = Demo
% Author: J. Selva.
% Date: July 2011.
% Version: 0.
st = [];
n = input('Options: 1 perform single estimation, 2 estimate RMS error: ');
switch n
case 1
st = TestSingleTrial;
case 2
TestRMSError;
end
function st = TestSingleTrial
%Scenario with 5 frequencies.
M = 2048;
SNRInSpecdB = 27; %Signal-to-noise ratio after correlation.
Nf = 5;
df = 1/M; %Minimum separation among frequencies.
MaxAmp = 1; %Maximum amplitude.
PFA = 10.^-5;
nVar = M*MaxAmp^2/10^(SNRInSpecdB/10);
nThres = NoiseThres(M,nVar,PFA);
MinAmp = 2*nThres/M;
[vf,va] = RandomFreqAndAmpVectors(Nf,df,MinAmp,MaxAmp);
y = DataRealization(vf,va,M,nVar);
p = Preprocessing(y,nVar,PFA);
disp('True frequencies: ')
disp(vf)
disp('True amplitudes: ')
disp(va)
disp(['Post-correlation Signal-to-noise ratio: ',...
num2str(SNRInSpecdB),' dB'])
s = input('To compute ML estimate, use (S)erial or (M)anual method? ','s');
if strcmp(lower(s),'s')
st = MethodSerial1(p);
else
st = MethodManual1(p);
end
disp('ML frequency estimates: ')
disp(st.Freqs)
disp('ML amplitude estimates: ')
disp(st.Amps)
function TestRMSError
%Estimates the RMS error and computes the Cramer-Rao bounds.
format compact
format short g
M = 2048;
SNRInSpecdB = 27;
Nf = 5;
df = 1/M;
MaxAmp = 1;
PFA = 10.^-5;
nVar = M*MaxAmp^2/10^(SNRInSpecdB/10);
nThres = NoiseThres(M,nVar,PFA);
MinAmp = 2*nThres/M;
[vf,va] = RandomFreqAndAmpVectors(Nf,df,MinAmp,MaxAmp);
crb = CRBBound(vf,va,M,nVar);
Ntr = 0;
S2 = zeros(Nf,1);
for k = 1:1e5
y = DataRealization(vf,va,M,nVar);
p = Preprocessing(y,nVar,PFA);
st = MethodSerial1(p);
if length(st.Freqs) == 5
Ntr = Ntr + 1;
S2 = S2 + (vf-st.Freqs).^2;
end
if mod(k,100) == 0
disp(['Number or trials: ', num2str(k)])
disp([' Frequencies: ' num2str(vf.')])
disp(['Deviation CRB bounds (dB): ' num2str(20*log10(crb.'))])
disp([' RMS errors (dB): ' num2str(10*log10(S2(:).'/Ntr))])
disp(' ');
end
end
function [vf,va] = RandomFreqAndAmpVectors(Nf,MinFreqSep,MinAmp,MaxAmp)
%Constructs random vectors of frequencies and amplitudes. The module of the largest
%amplitude is MaxAmp. There are no amplitudes with module below MinAmp, and the minimum
%separation among frequencies is MinFreqSep.
vf = sort(rand(Nf,1));
vf = vf + MinFreqSep*(cumsum(ones(Nf,1))-1);
vf = vf/(1+(Nf-1)*MinFreqSep) - 0.5;
pr = rand(Nf,1);
ind = floor(rand(1)/Nf)+1;
pr(ind) = 1;
va = (MinAmp+(MaxAmp-MinAmp)*pr) .* exp(1i*2*pi*rand(Nf,1));
function y = DataRealization(vf,a,M,nVar)
m1 = FirstIndex(M);
y = exp(1i*2*pi*(m1:m1+M-1).'*vf(:).')*a(:)+randn(M,2)*[1 ; 1i]*sqrt(nVar/2);
function p = Preprocessing(y,nVar,PFA)
%Performs the preprocessing. Fundamentally this consists in computing the FFT inside
%"Correlation.m".
%The input parameters are
%
% y: Data vector.
% nVar: Noise variance estimate.
% PFA: Probability of seeing a frequency when there isn't any. Usually a very low value
% like 1e-5.
p = [];
p.M = length(y);
p.nVar = nVar;
p.PFA = PFA;
p.nThres = NoiseThres(p.M,p.nVar,p.PFA); %Noise threshold.
icf = Correlation(y,p.nThres); %Struct with data for interpolating the correlation.
p.cfd = Cost(real(y(:)'*y(:)),icf); %Struct with data for computing the ML cost function ...
%and its derivatives.
p.freqThres = 1/(10*p.M); %If two frequencies differ in less than ...
%freqThres in the iterative search, then one of them is
%discarded. This prevents ill-conditioning.
p.freqStopThres = 1/(p.M*1e4); %Threshold for stopping the ML cost function minimization.
function st = MethodSerial1(p)
%Adds frequencies and computes ML estimates sequentially.
vfEst = [];
ModifiedEstimation = logical(1);
while ModifiedEstimation
[vf1,va1] = Correlation(p.cfd.icf,'gridMainLocalMaxima',vfEst);
ModifiedEstimation = ~isempty(vf1);
if ModifiedEstimation
vfEst(end+1) = vf1(1);
[listv,listLv,listCL,ErrRatio] = LookForMinimum(p.cfd,vfEst,...
10.^-8,p.freqThres,p.freqStopThres);
vfEst = sort(listv(:,end));
end
end
st = Report(p,vfEst);
function st = MethodManual1(p,varargin)
%Allows the user to set the initial frequency estimates on a plot.
if nargin == 1
vfEst = [];
else
vfEst = varargin{1};
end
vfEst = vfEst(:);
figure
subplot(2,1,1)
Correlation(p.cfd.icf,'PlotPeriodogramFit',vfEst);
yl = get(gca,'YLim');
subplot(2,1,2)
Correlation(p.cfd.icf,'PlotPeriodogramFitError',vfEst);
while 1
pr1 = input('A/a -> Add freq., D/d -> Delete freq, X/x -> exit: ','s');
switch lower(pr1)
case 'a'
[f1,pr] = ginput(1);
vfEst = sort([vfEst; f1]);
[listv,listLv,listCL,ErrRatio] = LookForMinimum(p.cfd,vfEst,10.^-8,...
p.freqThres,p.freqStopThres);
vfEst = sort(listv(:,end));
subplot(2,1,1)
Correlation(p.cfd.icf,'PlotPeriodogramFit',vfEst);
subplot(2,1,2)
Correlation(p.cfd.icf,'PlotPeriodogramFitError',vfEst);
case 'd'
[f1,pr] = ginput(1);
[pr,ind] = min(abs(f1-vfEst));
vfEst(ind) = [];
subplot(2,1,1)
Correlation(p.cfd.icf,'PlotPeriodogramFit',vfEst);
subplot(2,1,2)
Correlation(p.cfd.icf,'PlotPeriodogramFitError',vfEst);
case 'x'
st = Report(p,vfEst);
return
end
end
function thres = NoiseThres(M,nVar,PFA)
%Computes the noise threshold. The method will discard any correlation local maximum below
%this threshold.
thres = sqrt((M*nVar/2)*icdf('chi2',(1-PFA)^(1/M),2));
function m1 = FirstIndex(M)
m1 = -ceil(M/2);
function varargout = Correlation(varargin)
%Computes the correlation value and its derivatives by calling "Barycentric". It does
%other things like plotting.
if ~isstruct(varargin{1})
st = [];
temp = {[],[],10.^-15,1.5};
[temp{1:nargin}] = deal(varargin{:});
[a,st.nThres,epsilon,OvFMin] = deal(temp{:});
a = a(:).';
st.M = length(a);
st.K = 2^nextpow2(st.M*OvFMin);
st.bary = Barycentric(0,SampleFourierCorrelation1DVer1(a,st.K),...
-2*FirstIndex(st.M),1/st.K,whatP(1/st.K,-2*FirstIndex(st.M),epsilon));
st.IndActive = find(abs(st.bary.a) > st.nThres);
st.IndActive = unique([st.IndActive-1,st.IndActive,st.IndActive+1]);
if st.IndActive(1) == 0
st.IndActive(1) = [];
end
if st.IndActive(end) == st.K+1
st.IndActive(end) = [];
end
st.aTActive = st.bary.a(st.IndActive);
st.vfActive = centreFreq((st.IndActive-1)/st.K);
varargout = {st};
return
end
[st,msg] = deal(varargin{1:2});
switch msg
case 'value'
varargout = {Barycentric(st.bary,varargin{3:end})};
case 'gridValueInRange'
[ind1,ind2] = deal(varargin{3:4});
varargout = {st.bary.a(mod(ind1:ind2,st.K)+1)};
case 'ampEstimates'
vf = varargin{3};
if length(vf) == 0
varargout = {[]};
end
varargout = {RegGeoSum(st.M,vf,vf,0,0) \ Barycentric(st.bary,vf)};
case 'gridValue'
lInd = varargin{3};
if nargin == 4
vf = varargin{4};
else
vf = [];
end
if isempty(vf)
varargout = {st.bary.a(mod(lInd,st.K)+1)};
else
vf1 = centreFreq(lInd/st.K);
VRaa = RegGeoSum(st.M,vf,vf,0,0);
VRa1a = RegGeoSum(st.M,vf1,vf,0,0);
pr = VRa1a*(VRaa\Barycentric(st.bary,vf));
pr = reshape(pr,size(lInd));
varargout = {st.bary.a(mod(lInd,st.K)+1)-pr};
end
case 'gridFit'
lInd = varargin{3};
if nargin == 4
vf = varargin{4};
else
vf = [];
end
if isempty(vf)
varargout = {zeros(size(lInd))};
else
vf1 = centreFreq(lInd/st.K);
VRaa = RegGeoSum(st.M,vf,vf,0,0);
VRa1a = RegGeoSum(st.M,vf1,vf,0,0);
pr = VRa1a*(VRaa\Barycentric(st.bary,vf));
pr = reshape(pr,size(lInd));
varargout = {pr};
end
case 'plotErrorIndB'
if nargin == 2
vf = [];
else
vf = varargin{3};
end
s = 20*log10(abs(Correlation(st,'gridValue',0:st.K-1,vf)));
gf = centreFreq((0:st.K-1)/st.K);
[gf,ord] = sort(gf);
s = s(ord);
plot(gf,s,[-0.5,0.5],20*log10(abs(st.nThres)*[1,1]))
Ms = max(s);
set(gca,'YLim',[Ms-50,Ms+10])
xlabel('Frequency (f)')
ylabel('Error spectrum (dB)')
grid on
varargout = {};
case 'PlotPeriodogramFit'
vfEst = varargin{3};
if isempty(vfEst)
s0 = 20*log10(abs(Correlation(st,'gridValue',0:st.K-1)));
gf = centreFreq((0:st.K-1)/st.K);
[gf,ord] = sort(gf);
s0 = s0(ord);
plot(gf,s0)
hold on
nt = 20*log10(abs(st.nThres));
plot([-0.5,0.5],nt([1,1]).',...
'LineWidth',2,'LineStyle','--','Color',[1 0 0])
hold off
text(-0.49,nt+2,'Noise threshold')
xlabel('Frequency (f)')
ylabel('Amplitude (dB)')
title('Periodogram data, fit, and est. frequencies (stems)')
return
end
vaEst = Correlation(st,'ampEstimates',vfEst);
stem(vfEst,20*log10(st.M*abs(vaEst)),'g')
hold on
s0 = 20*log10(abs(Correlation(st,'gridValue',0:st.K-1)));
gf = centreFreq((0:st.K-1)/st.K);
[gf,ord] = sort(gf);
s0 = s0(ord);
plot(gf,s0,'b')
hold on
nt = 20*log10(abs(st.nThres));
plot([-0.5,0.5],nt([1,1]).',...
'LineWidth',2,'LineStyle','--','Color',[1 0 0])
s1 = 20*log10(abs(Correlation(st,'gridFit',0:st.K-1,vfEst)));
s1 = s1(ord);
plot(gf,s1,'r')
Ms = max(s0);
set(gca,'YLim',[Ms-50,Ms+10])
xlabel('Frequency (f)')
ylabel('Amplitude (dB)')
title('Periodogram fit')
hold off
varargout = {};
case 'PlotPeriodogramFitError'
if nargin == 2
vfEst = [];
else
vfEst = varargin{3};
end
vaEst = Correlation(st,'ampEstimates',vfEst);
s0 = 20*log10(abs(Correlation(st,'gridValue',0:st.K-1,vfEst)));
gf = centreFreq((0:st.K-1)/st.K);
[gf,ord] = sort(gf);
s0 = s0(ord);
nt = 20*log10(abs(st.nThres));
plot(gf,s0)
hold on
plot([-0.5,0.5],nt([1,1]).',...
'LineWidth',2,'LineStyle','--','Color',[1 0 0])
hold off
Ms = max(s0);
set(gca,'YLim',[Ms-50,Ms+10])
xlabel('Frequency (f)')
ylabel('Amplitude (dB)')
title('Periodogram error')
varargout = {};
case 'gridValueGivenAmps'
[lInd,vf,aEst] = deal(varargin{3:5});
vf1 = centreFreq(lInd/st.K);
VRa1a = RegGeoSum(st.M,vf1,vf,0,0);
pr1 = st.bary.a(mod(lInd,st.K)+1);
pr1 = reshape(pr1,size(lInd));
pr2 = VRa1a*aEst;
pr2 = reshape(pr2,size(lInd));
varargout = {pr1-pr2};
case 'gridLocalMaxima'
if nargin == 2
pr = [-2,-1,0].';
pr = abs(st.bary.a(mod(st.IndActive([1;1;1],:) + ...
pr(:,ones(1,length(st.IndActive))),st.K)+1));
pr = pr(1,:)<=pr(2,:) & pr(3,:)<=pr(2,:);
vf = st.vfActive(pr);
va = st.aTActive(pr);
varargout = {vf,va};
return
end
vf = varargin{3};
if length(vf) == 0
[varargout{1:2}] = Correlation(st,'gridLocalMaxima');
return
end
amp = Correlation(st,'ampEstimates',vf);
tmp = abs(Correlation(st,...
'gridValueGivenAmps',st.IndActive-1,vf,amp)) >= st.nThres;
ind = st.IndActive(tmp);
pr = [-2,-1,0].';
ind1 = ind([1,1,1].',:) + pr(:,ones(1,length(ind)));
pr = abs(Correlation(st,'gridValueGivenAmps',ind1,vf,amp));
pr = pr(1,:)<=pr(2,:) & pr(3,:)<=pr(2,:);
ind = ind(pr)-1;
if isempty(ind)
varargout = {[],[]};
else
varargout = {centreFreq(ind/st.K),Correlation(st,...
'gridValueGivenAmps',ind,vf,amp)};
end
case 'gridMainLocalMaxima'
[vf,va] = Correlation(st,'gridLocalMaxima',varargin{3:end});
if length(vf) == 0
varargout = {[],[]};
return
end
[varargout{1:2}] = selectLocalMaximaVer3(vf,va,st.nThres,st.M);
case 'gridMaximum'
pr = Correlation(st,'gridValue',0:st.K-1,varargin{3:end});
[ma,ind] = max(abs(pr));
varargout = {(ind-1)/st.K,ma};
end
function tab = SampleFourierCorrelation1DVer1(a,K)
M = length(a);
m1 = FirstIndex(M);
a = a(:).';
tab = fft([a(-m1+1:end),zeros(1,K-M),a(1:-m1)]);
function out = centreFreq(f)
out = f + ceil(-0.5-f);
function P = whatP(T,B,epsilon)
P = ceil(acsch(epsilon)/(pi*(1-B*T)));
function d = dist(x,y)
if x<y
d = min([y-x,x+1-y]);
else
d = min([x-y,y+1-x]);
end
function v = sincMask(x,y,M)
v = zeros(size(x));
for k = 1:length(x)
x1 = max([dist(x(k),y)-1/(2*M),0]);
v(k) = 1/max([M*sin(pi*x1),1]);
end
function v = slBoundVer2(x,vam,vfe,nThres,M1)
vam = vam(:);
vfe = vfe(:);
sll = sincMask(vfe,x,M1).'*vam;
v = max([sll+nThres,1.1*sll]);
function varargout = selectLocalMaximaVer3(vf,va,nThres,M)
va1 = abs(va(:));
[pr,ord] = sort(va1);
ord = flipud(ord);
vf1 = vf(ord);
va1 = va1(ord);
n1 = length(va1);
if va1(1) < nThres
varargout = {[],[]};
return
end
va2 = va1(1);
va2Abs = va1(1);
vf2 = vf1(1);
for n=2:n1
if va1(n)<nThres
continue
end
if va1(n) < slBoundVer2(vf1(n),va2Abs,vf2,nThres,M)
continue
end
vf2(end+1) = vf1(n);
va2(end+1) = va1(n);
va2Abs(end+1) = abs(va(ord(n)));
end
varargout = {vf2,va2};
function out = Barycentric(varargin)
%Performs barycentric interpolation using the method in [Bary].
if ~isstruct(varargin{1})
st = [];
[st.n0,st.a,st.B,st.T,st.P] = deal(varargin{:});
st.N = length(st.a);
st.alpha = baryWeights(st.T,st.B,st.P);
st.alpha = st.alpha([1:st.P,st.P+2:end,st.P+1]);
st.vt = (-st.P:st.P)*st.T;
st.vt = st.vt([1:st.P,st.P+2:end,st.P+1]);
out = st;
else
[st,v] = deal(varargin{1:2});
if nargin == 3
nd = varargin{3};
else
nd = 0;
end
v = v(:);
Nv = numel(v);
n = floor(v/st.T+0.5);
u = v - n*st.T;
n1 = n - st.n0;
pr = [-st.P:-1,1:st.P,0];
da = st.a(mod(n1(:,ones(1,2*st.P+1)) + pr(ones(Nv,1),:),st.N)+1);
out = DerBarycentricInterp3Vec(st.alpha,da,st.vt,u,nd);
end
function w = baryWeights(T,B,P)
vt = (-P:P)*T;
g = APPulse(vt,1/T-B,T*(P+1));
gam = gamma(vt/T+P+1) .* gamma(-vt/T+P+1) .* g;
N = length(vt);
LD = ones(1,N);
for k = 1:N
LD(k) = prod(vt(k)-vt(1:k-1))* prod(vt(k)-vt(k+1:N));
end
w = gam ./ LD;
w = w / max(abs(w));
function v = APPulse(t,B,TSL)
v = real(sinc(B*sqrt(t.^2-TSL^2)))/real(sinc(1i*pi*B*TSL));
function out = DerBarycentricInterp3Vec(alpha,zL,t,tauL,nd)
NtauL = size(tauL,1);
LBtau = 200; %Block size for vectorized processing. Change this to adjust the memory
%usage.
NBlock = 1+floor(NtauL/LBtau);
t = t(:).';
Nt = length(t);
out = zeros(NtauL,nd+1);
for r = 0:NBlock-1
ind1 = 1+r*LBtau;
ind2 = min([(r+1)*LBtau,NtauL]);
Ntau = ind2-ind1+1;
z = zL(ind1:ind2,:);
tau = tauL(ind1:ind2);
vD = zeros(Ntau,1);
LF = [1,zeros(1,Nt-1)];
LF = LF(ones(Ntau,1),:);
for k = 1:Nt-1
vD = vD .* (tau-t(k))+alpha(k)*LF(:,k);
LF(:,k+1) = LF(:,k) .* (tau-t(k));
end
vD = vD .* (tau-t(Nt)) + alpha(Nt) * LF(:,Nt);
for kd = 0:nd
pr = z(:,end); z1 = z-pr(:,ones(1,Nt)); cN = zeros(Ntau,1);
for k = 1:Nt-1
cN = cN .* (tau-t(k)) + z1(:,k) .* alpha(k) .* LF(:,k);
end
cN = cN ./ vD;
ztau = z(:,end)+(tau-t(end)) .* cN;
out(ind1:ind2,kd+1) = ztau;
if kd < nd
pr1 = ztau;
pr1 = z(:,1:end-1) - pr1(:,ones(1,Nt-1));
pr2 = t(1:end-1);
pr2 = pr2(ones(Ntau,1),:)-tau(:,ones(1,Nt-1));
z = [pr1 ./ pr2, cN] * (kd+1);
end
end
end
function varargout = Cost(varargin)
%Computes the ML cost function differentials by calling "Correlation" and "RegGeoSum".
if ~isstruct(varargin{1})
st = [];
[st.EnergyRyy,st.icf] = deal(varargin{:});
varargout = {st};
return
end
st = varargin{1};
if nargin > 2
vf = varargin{3};
else
vf = [];
end
switch varargin{2}
case 'value'
vf = varargin{3};
Rss = RegGeoSum(st.icf.M,vf,vf,0,0);
Rsd = RegGeoSum(st.icf.M,vf,vf,0,1);
Rdd = RegGeoSum(st.icf.M,vf,vf,1,1);
pr = Correlation(st.icf,'value',vf,1);
Rsy = pr(:,1);
Rdy = pr(:,2);
pr = Rss\[Rsy,Rsd];
MRsy = pr(:,1);
MRsd = pr(:,2:end);
varargout = {st.EnergyRyy-RealTrace(Rsy',MRsy),-2*RealDiag(MRsy,Rdy'-MRsy'*Rsd),...
2*RealHadamard((Rdd'-Rsd'*MRsd).',MRsy*MRsy')};
case 'value0'
if length(vf) == 0
varargout = {st.EnergyRyy};
return
end
vf = varargin{3};
Rss = RegGeoSum(st.icf.M,vf,vf,0,0);
Rsy = Correlation(st.icf,'value',vf);
[varargout{1}] = st.EnergyRyy-RealTrace(Rsy',Rss\Rsy);
end
function v = RealHadamard(A,B)
v = real(A).*real(B)-imag(A).*imag(B);
function C = RealDiag(A,B)
% Same as real(diag(A*B)) but using less flops.
C = RealHadamard(A.',B);
if size(C,1) == 1
C = C.';
else
C = sum(C).';
end
function C = RealTrace(A,B)
% Same as real(trace(A*B)) but using less flops.
C=RealHadamard(A.',B);
C=sum(C(:));
function [listv,listLv,listCL,ErrRatio] = ...
LookForMinimum(L,v0,nThres,freqThres,freqStopThres,varargin)
%Executes a descent algorithm using a method similar to that in [ML].
if ~isempty(varargin)
nItMax = varargin{1};
else
nItMax = 20;
end
v0 = v0(:);
Nf = numel(v0);
listv = [v0,zeros(Nf,nItMax)];
listLv = [Cost(L,'value0',v0),zeros(1,nItMax)];
listCL = [NaN,zeros(1,nItMax)];
succ = 1;
nIt = 1;
while succ && nIt <= nItMax
[succ,v,Lv,CL] = NewtonIt(L,listv(:,nIt),freqThres);
if succ
nIt = nIt + 1;
listv(:,nIt) = v;
listLv(nIt) = Lv;
listCL(nIt) = CL;
else
break
end
if max(abs(listv(:,nIt)-listv(:,nIt-1))) < freqStopThres
break
end
end
ErrRatio = (Cost(L,'value0')-listLv(end))/listLv(end);
listv(:,nIt+1:end) = [];
listLv(nIt+1:end) = [];
listCL(nIt+1:end) = [];
function [v,vf] = validVector(vf,thres)
vf = sort(vf(:));
if length(vf)==1
v = 1;
return
end
m = vf(2:end)-vf(1:end-1) < thres;
if any(m)
v = 0;
vf([logical(0); m]) = [];
else
v = 1;
end
function [succ,v1,Lv1,nCutsLeft] = NewtonIt(L,v0,thres)
Lv1 = [];
nCutsLeft = [];
[vv,v1] = validVector(v0,thres);
if ~vv
succ = 0;
return
end
[Lv0,g0,h0] = Cost(L,'value',v0);
d = -h0\g0;
nCutsLeft = 20;
while nCutsLeft > 0
[vv,v1] = validVector(v0+d,thres);
if ~vv
succ = 0;
v1 = v0;
end
Lv1 = Cost(L,'value0',v0+d);
if Lv1 < Lv0
succ = 1;
v1 = v0 + d;
break;
end
nCutsLeft = nCutsLeft - 1;
d = d/2;
end
if nCutsLeft == 0
succ = 0;
end
function out = RegGeoSum(varargin)
%Computes the value of a geometric sum or its differentials using the exact formula. It
%takes care of the preventable singularity at zero.
switch nargin
case 2
[M,vf] = deal(varargin{:});
nd = 0;
msg = 'value';
case 3
[M,vf,nd] = deal(varargin{:});
msg = 'value';
case 5
[M,vf1,vf2,nd1,nd2] = deal(varargin{:});
msg = 'corr';
end
switch msg
case 'value'
out = RegGeoSum1(M,vf,nd);
case 'corr'
vf1 = vf1(:);
Nf1 = length(vf1);
vf2 = vf2(:);
Nf2 = length(vf2);
out = RegGeoSum1(M,-vf1(:,ones(1,Nf2))+vf2(:,ones(1,Nf1)).',nd1+nd2)*(-1)^nd1;
end
function v = RegGeoSum1(M,vf,nd)
m1 = -ceil(M/2);
ind = abs(vf)>eps^(1/3);
vf1 = vf(ind);
v = zeros(size(vf));
switch nd
case 0
v(ind) = -(exp(1i*2*pi*m1*vf1)-exp(1i*2*pi*(m1+M)*vf1))./(-1+exp(1i*2*pi*vf1));
v(~ind) = M;
case 1
v(ind) = (-2*1i*pi*( exp(1i*2*pi*vf1*(1 + m1))*(-1 + m1) - exp(1i*2*pi*vf1*m1)*m1 ...
-exp(1i*pi*2*vf1*(1 + M + m1))*(-1 + M + m1) + ...
exp(1i*2*pi*vf1*(M + m1))*(M + m1)))./(-1 + exp(1i*2*pi*vf1)).^2;
v(~ind) = pi*1i*M*(-1+M+2*m1);
case 2
v(ind) = (4*(exp(1i*2*pi*vf1*(2 + m1))*(-1 + m1)^2 + exp(1i*2*pi*vf1*m1)*m1^2 - ...
exp(1i*2*pi*vf1*(2 + M + m1))*(-1 + M + m1)^2 - ...
exp(1i*2*pi*vf1*(M + m1))*(M + m1)^2 + ...
exp(1i*2*pi*vf1*(1 + m1))*(1 - 2*(-1 + m1)*m1) + ...
exp(1i*2*pi*vf1*(1 + M + m1))* ...
(-1 + 2*M^2 + 2*(-1 + m1)*m1 + M*(-2 + 4*m1)))*pi^2)./ ...
(-1 + exp(1i*2*pi*vf1)).^3;
v(~ind) = -(2*pi^2/3)*M*(1+6*m1^2+6*m1*(-1+M)-3*M+2*M^2);
end
function st = Report(p,vfEst)
st = [];
st.Freqs = vfEst;
st.Amps = Correlation(p.cfd.icf,'ampEstimates',vfEst);
function FreqsDevBounds = CRBBound(vf,va,M,nVar)
%Computes the Cramer-Rao bound of the frequencies.
va = va(:);
vf = vf(:);
Raa = RegGeoSum(M,vf,vf,0,0);
Rad = RegGeoSum(M,vf,vf,0,1);
Rdd = RegGeoSum(M,vf,vf,1,1);
FreqsDevBounds = sqrt((nVar/2)*diag(inv(real((conj(va)*va.') .* (Rdd-Rad'*(Raa\Rad))))));
function st = MLEstimate(y,nVar,PFA)
p = Preprocessing(y,nVar,PFA);
st = MethodSerial1(p);
st.CRBBoundsOnFreqDevs = CRBBound(st.Freqs,st.Amps,length(y),nVar);
Jesus Selva ● June 15, 2011
● Coded in
Matlab function Demo
%Author: J. Selva.
%Date: June 2011.
close all
format short g
format compact
T = 1; %Sampling period in both dimensions.
B = 1/(1.223*T); %This is the two-sided bandwidth in both dimensions. It is a typical ...
%value in SAR image co-registration.
P = 18; %Kernel semi-length.
Q = 10; %Polynomial order in Farrow structure.
%The next call to "ChebApproxPolys" computes all Farrow coefficients.
t0 = -P*T-T/2; %The parameters t0, Nt, 2*P+1 specify the 2*P+1 intervals in which Knab's ...
%pulse is to be approximated by polynomials.
Nt = 2*P+1;
Dt = T;
IOut = [-T/2,T/2];
%Perform Chebyshev interpolation. The parameters T, B and P are passed to KnabAPPulse.
K = ChebApproxPolys('KnabAPPulse',t0,Nt,Dt,IOut,40,Q,T,B,P).';
%Evaluate the test image in a regular grid with spacing 1.
Lx = 100;
Ly = 100;
x0 = (0:Lx-1)';
x0 = x0(:,ones(1,Lx));
y0 = x0.';
x0 = x0(:);
y0 = y0(:);
disp('Computing image in a regular grid...')
tic;
Image = reshape(TestImage5(x0,y0,B),[Lx,Ly]);
disp(['(Timing: ',num2str(toc), ' sec.)'])
disp(' ')
clear x0 y0
%Specify NP points on the image at random.
NP = 10000;
x = P*T+rand(1,NP)*(Lx-2*P)*T;
y = P*T+rand(1,NP)*(Ly-2*P)*T;
%Evaluate image directly.
disp('Computing image at selected positions exactly...')
tic;
vRef = TestImage5(x,y,B);
disp(['(Timing: ',num2str(toc), ' sec.)'])
disp(' ')
%Interpolate the image using the 2D-Farrow structure.
disp('Computing image at selected positions using Farrow 2-D...')
tic;
v = Farrow2DInterp(Image,K,x,y);
disp(['(Timing: ',num2str(toc), ' sec.)'])
disp(' ')
%Show the error.
disp(['Maximum interpolation error: ', ...
num2str(20*log10(max(abs(vRef-v))/max(abs(vRef)))),' dB'])
disp(' ')
disp('Display image and interpolation error:')
disp(' ')
%Evaluate the test image in a regular grid with spacing 1.
Lx = 100;
Ly = Lx;
x1 = (P*T:0.5*T:(Lx-P)*T)';
x1 = x1(:,ones(1,length(x1)));
y1 = x1.';
x1 = x1(:);
y1 = y1(:);
tic;
disp('Computing image in over-sampled grid ...')
vRef = TestImage5(x1,y1,B);
disp(['(Timing: ',num2str(toc), ' sec.)'])
disp(' ')
tic;
disp('Computing signal in over-sampled grid using 2-D Farrow ...')
v = Farrow2DInterp(Image,K,x1,y1);
disp(['(Timing: ',num2str(toc), ' sec.)'])
disp(' ')
x1 = (P*T:0.5*T:(Lx-P)*T)';
y1 = x1.';
disp(['Maximum interpolation error: ', ...
num2str(20*log10(max(abs(vRef-v))/max(abs(vRef)))),' dB'])
vRef = reshape(vRef,[length(x1),length(y1)]);
v = reshape(v,[length(x1),length(y1)]);
mesh(x1,y1,abs(vRef))
xlabel('x coordinate')
ylabel('y coordinate')
zlabel('Image''s abs. value')
title('Test Image (abs. value)')
figure
mesh(x1,y1,20*log10(abs(v-vRef)/max(abs(vRef(:)))))
xlabel('x coordinate')
ylabel('y coordinate')
zlabel('Interp. error (dB)')
title('Interpolation error (dB)')
function P = ChebApproxPolys(Func,t0,Nt,Dt,IOut,NPoints,Q,varargin)
%This function constructs Chebyshev interpolation polynomials in a set of consecutive ...
%intervals.
%
% Func is the function to be interpolated.
% t0 is the first instant of the first interval.
% Nt is the number of intervals.
% IOut is the interval in which the output polynomials will be evaluated.
% Dt is the interval width. The intervals are [t0,t0+Dt], [t0+Dt,t0+2*Dt],....
% [t0+(Nt-1)*Dt,t0+Nt*Dt].
% NPoints is the number of points in which Func is evaluated, for each interval.
% Q is the output polynomial order.
% varargin contains extra parameters which are passed to varargin.
%Construct matrix that maps function values to polynomial coefficients.
%The algorithm is that in section 9.8 of
%
% W. H. Press et al., Numerical Recipes in C. Cambridge University
% Press, 1997.
M = MapPolyMatrix([-1,1],IOut,Q) * ChebCoefMatrix(Q) * ...
BlockChebInvMatrix(NPoints,Q);
% Get Chebyshev roots.
r = ChebRoots(NPoints);
r = r(:);
P = zeros(Q,Nt);
%Compute polynomial for each interval.
for k = 0:Nt-1,
v1 = feval(Func,MapIaToIb(r,[-1,1],t0+[k,k+1]*Dt),varargin{:});
P(:,k+1) = double(M * v1(:));
end
function M = MapPolyMatrix(IIn,IOut,Q)
%For a polynomial with Q coefficients which is evaluated in interval IIn, this function ...
%gives the matrix that transforms its coefficients, so that it can be evaluated in the ...
%associated points in interval IOut.
ayx = (IIn(1)-IIn(2))/(IOut(1)-IOut(2));
byx = IIn(1)-ayx*IOut(1);
if byx == 0
M = diag(ayx.^(0:Q-1));
else
M = zeros(Q);
for q = 0:Q-1
for r = 0:q
M(r+1,q+1) = nchoosek(q,r) * ayx^r * byx^(q-r);
end
end
end
function M = ChebCoefMatrix(N)
M = zeros(N);
M(1) = 1;
if N > 1
M(2,2) = 1;
end
for n = 3:N,
M(2:n,n) = 2*M(1:n-1,n-1);
M(1:n-2,n) = M(1:n-2,n) - M(1:n-2,n-2);
end
function M = BlockChebInvMatrix(N,Q)
r = ChebRoots(N);
r= r(:);
M = zeros(N,Q);
if Q >=1
M(:,1) = 1;
end
if Q >= 2
M(:,2) = r;
end
for q = 3:Q,
M(:,q) = 2 * r .* M(:,q-1) - M(:,q-2);
end
M = M.'/N;
M(2:end,:) = M(2:end,:) * 2;
function R = ChebRoots(N)
R = cos(pi*((1:N)-1/2)/N);
function y = MapIaToIb(x,Ia,Ib)
%Linear transformation that maps interval Ia onto interval Ib.
y = (Ib(1)+Ib(2))/2 + (x-(Ia(1)+Ia(2))/2) * ((Ib(2)-Ib(1))/(Ia(2)-Ia(1)));
function [g,BoundAP] = KnabAPPulse(t,T,B,P)
%Knab interpolation pulse in
%
% J. J. Knab, 'Interpolation of band-limited functions using the Approximate
% Prolate series', IEEE Transactions on Information Theory, vol. IT-25,
% no. 6, pp. 717-720, Nov 1979.
%
%Grid points -P*T:T:P*T.
%Interpolation interval: [-T/2,T/2].
g = sinc(t/T) .* sinc((1/T-B)*sqrt(t.^2-(P*T)^2)) / sinc(j*(1-B*T)*P);
if nargout == 2
BoundAP = 1/sinh(P*pi*(1-B*T));
end
function v = TestImage5(x,y,B)
%Test image formed by a sum of random sincs.
Seed = rand('state');
rand('state',362436069);
L = 100; %Image size in both dimensions.
Np = 5000;
p = rand(Np,2)*L;
A = exp(j*2*pi*rand(Np,1));
rand('state',Seed);
v = zeros(size(x));
for k = 1:length(x(:))
v(k) = A(:).' * (sinc(B*(x(k)-p(:,1))) .* sinc(B*(y(k)-p(:,2))));
end
function vx = Farrow2DInterp(Image,Kernel,x,y)
%Farrow interpolation in two dimensions using interlaced FFTs.
vxSize = size(x);
[LK,NK] = size(Kernel);
[Lx,Ly] = size(Image);
nx = round(x);
ny = round(y);
ux = x-nx;
uy = y-ny;
np = 1+nx+ny*Lx;
np = np(:);
clear x y nx ny
P = (LK-1)/2;
NFFTx = 2^ceil(log2(Lx+LK-1));
NFFTy = 2^ceil(log2(Ly+LK-1));
KernelFx = fft(Kernel,NFFTx);
KernelFy = fft(Kernel,NFFTy);
ImageFx = fft(Image,NFFTx);
vx = zeros(size(np))+j*zeros(size(np));
for kx = NK:-1:1,
ImageCx = ifft(ImageFx .* KernelFx(:,kx(ones(1,Ly))),NFFTx);
ImageCx = ImageCx(1+P:Lx+LK-1-P,:);
ImageCxFy = fft(ImageCx,NFFTy,2);
vy = zeros(size(np))+j*zeros(size(np));
for ky = NK:-1:1,
ImageCxCy = ifft(ImageCxFy .* KernelFy(:,ky(ones(1,Lx))).',NFFTy,2);
ImageCxCy = ImageCxCy(:,1+P:Ly+LK-1-P);
vy(:) = vy(:) .* uy(:) + ImageCxCy(np);
end
vx(:) = vx(:) .* ux(:) + vy(:);
end
vx = reshape(vx,vxSize);
Jesus Selva ● May 18, 2011
● ● 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
Jesus Selva ● April 18, 2011
● Coded in
Matlab function Demo(L)
%Author: J. Selva.
%Date: April 2011.
%This demo constructs trials of an undamped exponential plus noise, and then computes the
%maximum likelihood estimate of its frequency, amplitude, and phase. Only one FFT of size
%2^ceil(log2(st.M)+1) is computed, plus a few operations. The complexity is that given by
%the FFT. The algoritm is explained in
%
% http://arxiv.org/abs/1104.3069
%
format long g
format compact
M = 1024; %Number of samples.
gdB = 20; %SNR at matched filter output.
ARef = 1.3; %Amplitude.
fRef = 0.372; %Frequency in [-0.5,0.5[.
PhiRef = 0.43; %Phase (rad).
if nargin > 0
for k = 1:2:length(L)
eval([L{k} '=L{k+1};']);
end
end
%Define a random vector with 1024 elements.
SNRdB = gdB - 10*log10(M); %Signal-to-noise ratio (dB).
%Complex sinusoidal amplitude, frequency, and phase.
NVar = ARef^2*10.^(-SNRdB/10); %Noise variance.
m1 = -ceil(M/2);
while 1
disp(' ')
disp(['M = ', num2str(M), ', gdB = ', num2str(gdB), ' dB, ARef = ', num2str(ARef), ...
', fRef = ', num2str(fRef), ', PhiRef = ', num2str(PhiRef), ' rad'])
disp(' ')
%Set up a trial.
a = ARef*exp(1i*2*pi*fRef*(m1:m1+M-1).'+1i*PhiRef)+(randn(M,1)+1i*randn(M,1))*sqrt(NVar/2);
%Call the nufft function. The output is a struct with some variables for later use.
st = nufft(a);
%Compute the DFT peak.
[f1,A1,Phi1,L1,Deg,nIt] = MLEstimate(st);
%The output is the following,
%
% f1, A1, and Phi1 are the ML estimates of frequency, amplitude and phase, respectively.
%
% L1 is the ML cost function value at f1 and its first two derivatives. L1(1) should be
% large since it is the square value of the peak amplitude. L1(2) is the cost
% function derivative, and it should be small since we are at a critical
% point. Finally L1(2) should be a negative number with large amplitude since f1 is
% the abcissa of a maximum (second derivative).
%
% Deg: If Deg==1 then the Newton iteration was cut too many times.
% nIt: Number of Newton iterations.
disp('Exact and estimated frequency, amplitude, and phase:')
disp([fRef,ARef,PhiRef;f1,A1,Phi1])
disp(['Number of iterations: ' num2str(nIt)])
%
%Plot the spectrum close to the ML estimate
%
f = fRef+(-10:.05:10)/st.K;
sp = nufft(st,f);
plot(f,20*log10(abs(sp)/M))
hold on
stem(f1,20*log10(A1),'r')
hold off
grid on
set(gca,'YLim',[-40,10])
xlabel('f')
ylabel('Amplitude (dB)')
text(fRef+2/st.K,5,{['Freq. ' num2str(f1)],['Ampl. ' num2str(A1)],...
['Phase. (rad) ' num2str(Phi1)]})
%
R = input('Press Enter for a new trial, q to quit: ','s');
if any(strcmp(R,{'q','Q'}))
break
end
end
%------------------------------------------------------
%------------------------------------------------------
function [f1,A1,Phi1,L1,Deg,nIt] = MLEstimate(st)
%This function finds out the peak of a DFT. Its input is the struct produced by the nufft function.
%Find the spectral peak in the frequency grid.
[pr11,ind] = max(abs(st.aF).^2);
%Compute its frequency.
f0 = (ind-1)/st.K;
f0 = f0+ceil(-0.5-f0);
%L0 contains the differentials of the ML cost function of orders 0,1, 2 ...
%at f0.
L0 = LDiff(nufft(st,f0,2));
nItMax = 20;
err = 10.^-12;
nIt = 0;
%Call Iter (Newton iteration) until there is no change or the solution is degenerate.
while 1
[Deg,pr11,f1,L1] = Iter(st,f0,L0);
if Deg || abs(f1-f0) <= err || nIt >= nItMax
break
end
f0 = f1;
L0 = L1;
nIt = nIt+1;
end
pr = nufft(st,f1);
A1 = abs(pr)/st.M;
Phi1 = angle(pr);
function L = LDiff(c)
%Computes the differentials of the ML cost function from the nufft differentials.
L = [real(c(1)*conj(c(1))); 2*real(conj(c(1))*c(2)); ...
2*real(conj(c(2))*c(2)+conj(c(1))*c(3))];
function [Deg,nCut,f1,L1] = Iter(st,f0,L0)
%Performs the Newton iteration, cutting the Newton direction if necessary.
f1 = f0-L0(2)/L0(3); %Newton iteration.
f1 = f1+ceil(-0.5-f1); %Make it fall in [-0.5,0.5].
L1 = LDiff(nufft(st,f1,2)); %Compute the differentials at f1.
nCutMax = 20; %The difference f1-f0 will be reduced in the sequel if there is no ...
%improvement in the cost function
nCut = 0;
d = 1;
%In this loop the difference f1-f0 is halved a maximum of nCutMax times until there is ...
%an increase in the cost function
while L1(1) < L0(1) && nCut < nCutMax
d = d/2;
f1 = f0-d*L0(2)/L0(3);
f1 = f1+ceil(-0.5-f1);
L1 = LDiff(nufft(st,f1,2));
nCut = nCut+1;
end
Deg = nCut >= nCutMax;
%----------------------------------------------------------------
%----------------------------------------------------------------
%Non-uniform FFT function. This code is discussed in
% http://www.dsprelated.com/showcode/159.php
function varargout = nufft(varargin)
%Author: J. Selva
%Date: April 2011.
%
%nufft computes the FFT at arbitrary frequencies using barycentric interpolation.
%
%For the interpolation method, see
%
%J. Selva, 'Design of barycentric interpolator for uniform and nonuniform sampling
% grids', IEEE Trans. on Signal Processing, vol. 58, n. 3, pp. 1618-1627, March 2010.
%
%Usage:
%
% -First call nufft with the vector of samples as argument,
%
% st = nufft(a); %a is the vector of samples.
%
% The output is an struct. The field st.aF is the DFT of vector a, sampled with spacing
% 1/st.K.
%
% The DFT is defined using
%
% A_k = \sum_{n=m_1}^{m_1+M-1} a_n e^{-j 2 \pi n f}
%
% where m_1=-ceil(M/2) and M the number of samples.
%
% -Then call nufft with sintax
%
% out = nufft(st,f,nd,method)
%
% where
%
% f: Vector of frequencies where the DFT is evaluated. Its elements must follow
% abs(f(k))<=0.5
%
% nd: Derivative order. nufft computes derivatives up to order nd. At the output
% out(p,q) is the (q-1)-order derivative of the DFT at frequency f(p). The first
% column of out contains the zero-order derivatives, i.e, the values of the DFT at
% frequencies in vector f.
% method: If omitted, it is 'baryVec'. Four methods have been implemented:
%
% +'baryLoop': The output is interpolated using the barycentric method and a loop
% implementation.
%
% +'baryVec': The output is interpolated using the barycentric method and a
% vectorized implementation.
%
% +'directLoop': The output is computed using the DFT sum directly and a loop
% implementation.
%
% +'directVec': The output is computed using the DFT sum directly and a vectorized
% implementation.
if ~isstruct(varargin{1})
st = [];
st.a = varargin{1};
st.a = st.a(:);
st.M = length(st.a);
st.m1 = -ceil(st.M/2);
st.K = 2^ceil(log2(st.M)+1);
st.aF = fft([st.a(-st.m1+1:st.M); zeros(st.K-st.M,1); st.a(1:-st.m1)]);
errb = 10.^-13; %This is the interpolation accuracy. You can set any larger value, and ...
%this would reduce st.P. The number of interpolation samples is 2*st.P+1.
st.T = 1/st.K;
st.B = -2*st.m1;
st.P = ceil(acsch(errb)/(pi*(1-st.B*st.T)));
st.vt = MidElementToEnd((-st.P:st.P)*st.T);
st.alpha = MidElementToEnd(baryWeights(st.T,st.B,st.P));
varargout = {st};
return
end
[st,f] = deal(varargin{1:2});
nd = 0;
if nargin > 2
nd = varargin{3};
end
method = 'baryVec';
if nargin > 3
method = varargin{4};
end
Nf = length(f);
out = zeros(Nf,nd+1);
switch method
case 'baryLoop' %Interpolated method using loops
for k = 1:length(f)
x = f(k);
n = floor(x/st.T+0.5);
u = x - n * st.T;
da = MidElementToEnd(st.aF(1+mod(n-st.P:n+st.P,st.K)).');
out(k,:) = DerBarycentricInterp3(st.alpha,da,st.vt,u,nd);
end
case 'baryVec' %Vectorized interpolated method
f = f(:);
Nf = length(f);
n = floor(f/st.T+0.5);
u = f - n * st.T;
pr = [-st.P:-1 , 1:st.P , 0];
ms = st.aF(1+mod(n(:,ones(1,2*st.P+1)) + pr(ones(Nf,1),:),st.K));
if length(f) == 1
ms = ms.';
end
out = DerBarycentricInterp3Vec(st.alpha,ms,st.vt,u,nd);
case 'directLoop' % Direct method using loops
for k = 1:length(f)
out(k,1) = 0;
for r = st.m1:st.m1+st.M-1
out(k,1) = out(k,1)+exp(-1i*2*pi*r*f(k))*st.a(r-st.m1+1);
end
for kd = 1:nd
out(k,kd+1) = 0;
for r = st.m1:st.m1+st.M-1
out(k,kd+1) = out(k,kd+1) + ...
((-1i*2*pi*r).^kd .* exp(-1i*2*pi*r*f(k)))*st.a(r-st.m1+1);
end
end
end
case 'directVec' %Vectorized direct method
for k = 1:length(f)
out(k,1) = exp(-1i*2*pi*(st.m1:st.m1+st.M-1)*f(k))*st.a;
for kd = 1:nd
out(k,kd+1) = ...
((-1i*2*pi*(st.m1:st.m1+st.M-1)).^kd .* ...
exp(-1i*2*pi*(st.m1:st.m1+st.M-1)*f(k)))*st.a;
end
end
end
varargout = {out};
function v = MidElementToEnd(v)
ind = ceil(length(v)/2);
v = [v(1:ind-1),v(ind+1:end),v(ind)];
function v = APPulse(t,B,TSL)
v = real(sinc(B*sqrt(t.^2-TSL^2)))/real(sinc(1i*pi*B*TSL));
function w = baryWeights(T,B,P)
vt = (-P:P)*T;
g = APPulse(vt,1/T-B,T*(P+1));
gam = gamma(vt/T+P+1) .* gamma(-vt/T+P+1) .* g;
N = length(vt);
LD = ones(1,N);
for k = 1:N
LD(k) = prod(vt(k)-vt(1:k-1))* prod(vt(k)-vt(k+1:N));
end
w = gam ./ LD;
w = w / max(abs(w));
function out = DerBarycentricInterp3(alpha,s,t,tau,nd)
vD = 0;
Nt = length(t);
LF = [1,zeros(1,Nt-1)];
out = zeros(1,nd+1);
for k = 1:Nt-1
vD = vD*(tau-t(k))+alpha(k)*LF(k);
LF(k+1) = LF(k)*(tau-t(k));
end
vD = vD * (tau-t(Nt)) + alpha(Nt) * LF(Nt);
z = s;
for kd = 0:nd
z1 = z-z(end); cN = 0;
for k = 1:Nt-1
cN = cN * (tau-t(k))+z1(k)*alpha(k)*LF(k);
end
cN = cN/vD;
ztau = z(end)+(tau-t(end))*cN;
out(kd+1) = ztau;
if kd < nd
z = [ (z(1:end-1)-ztau)./(t(1:end-1)-tau) , cN ] * (kd+1);
end
end
function out = DerBarycentricInterp3Vec(alpha,zL,t,tauL,nd)
NtauL = size(tauL,1);
LBtau = 200; %Block size for vectorized processing. Change this to adjust the memory
%usage.
NBlock = 1+floor(NtauL/LBtau);
Nt = length(t);
out = zeros(NtauL,nd+1);
for r = 0:NBlock-1
ind1 = 1+r*LBtau;
ind2 = min([(r+1)*LBtau,NtauL]);
Ntau = ind2-ind1+1;
z = zL(ind1:ind2,:);
tau = tauL(ind1:ind2);
vD = zeros(Ntau,1);
LF = [1,zeros(1,Nt-1)];
LF = LF(ones(Ntau,1),:);
for k = 1:Nt-1
vD = vD .* (tau-t(k))+alpha(k)*LF(:,k);
LF(:,k+1) = LF(:,k) .* (tau-t(k));
end
vD = vD .* (tau-t(Nt)) + alpha(Nt) * LF(:,Nt);
for kd = 0:nd
pr = z(:,end); z1 = z-pr(:,ones(1,Nt)); cN = zeros(Ntau,1);
for k = 1:Nt-1
cN = cN .* (tau-t(k)) + z1(:,k) .* alpha(k) .* LF(:,k);
end
cN = cN ./ vD;
ztau = z(:,end)+(tau-t(end)) .* cN;
out(ind1:ind2,kd+1) = ztau;
if kd < nd
pr1 = ztau;
pr1 = z(:,1:end-1) - pr1(:,ones(1,Nt-1));
pr2 = t(1:end-1);
pr2 = pr2(ones(Ntau,1),:)-tau(:,ones(1,Nt-1));
z = [pr1 ./ pr2, cN] * (kd+1);
end
end
end
