CRLB Range Estimation for GNSS meta-signals

Hi,

I try to compute CRLB for range estimation but the results are not convincing. In the following code, I generate a meta-signal and apply the following equation to compute CRLB estimation for certain CN0 values.

From Steven.M.Kay - Fundamentals of Statistical Signal Processing Vol. 1 the CRLB for range estimation is:

where, F^2 and E are Gabor bandwidth and energy of the signal respectly.

I attach my code to obtain the simulation.

clear;
close all;
clc;

c       =   physconst('LightSpeed');

Tc      =   1e-3/1023;                  % Code CA time (~1us).
Ts      =   2*Tc;                       % Auto-correlation of CA chip duration.
Ti      =   1e-3;                       % Integration time (1ms).
fsc     =   1.5/Tc;                     % Sub-carrier frequency.
A       =   1;                          % Amplitude.
wsc     =   2*pi*fsc;                   % Angular frequency.
Npi     =   100;                        % #Samples.
Niter   =   1e3;                        % #Realizations.
Ti      =   1e-3;                       % Integration time.

% Autocorrelation of the CA chips.

if mod(Npi,2) == 0
    Npi = Npi + 1
end

tau     =   Tc*linspace(-1,1,Npi);
%tau     =   [tau(1:Npi/2) 0 tau(Npi/2 + 1:end)];    % Adjust time axis.
Rx      =   A - abs(tau)/Tc;
Rx2     =   Rx;

figure;
plot(tau,Rx);
title("Autocorrelation CA chip");

R_p     =   Rx.*exp(1i*wsc*tau);
R_n     =   Rx.*exp(-1i*wsc*tau);

figure;
plot(R_p(Npi/2 + 1:end), 'b');
hold on
plot(R_n(Npi/2 + 1:end), 'r');
legend("R_{+}(\tau) for 0\leq\tau\leqT_{c}","R_{-}(\tau) for 0\leq\tau\leqT_{c}");
xlabel('Real');
ylabel('Imag');

%%

R_ms    =   R_n + R_p;

figure;
plot(tau,real(R_ms));
title("Real[R_{MS}"+"(\tau)]");

figure;
plot(tau,abs(R_ms));
title("|R_{MS}"+"(\tau)|");

figure;
plot(tau,imag(R_ms));
title("Imag[R_{MS}]"+"(\tau)]");

% ========== CN0 =========

CN0     =   linspace(20,45,Npi);        % CN0 range values (dB-Hz).

% Spectral density: Fourier transform and change the axis so zero is at the center and have negative and positive
% frequencies (dual band).

Sx      =   fftshift(fft(R_ms));

Fs      =   1/(tau(2)-tau(1));
f       =   (Fs/2)*linspace(-1,1,Npi);

figure;
plot(f,abs(Sx));
title('x(t) spectral density');
xlabel('Frequency [Hz]');
ylabel('Power [W]');

% ========== CRLB ==========

num     =   trapz(f,((2*pi*f).^2).*(abs(Sx).^2));

E       =   trapz(tau,Rx.^2);
Px      =   E/Ts;

No      =   Px./db2pow(CN0);
Pno     =   No./Ti;

F2      =   num/trapz(f,abs(Sx).^2);

CRLB    =   1./((E./(Pno/2))*F2);

% ============================================================

tau_est     =   zeros(1,Niter);
tau_hat     =   zeros(1,Npi);
e           =   zeros(1,Npi);
delay       =   0;

Pn          =   No./Ti;

for ii=1:Npi
    sig2    =   Pn(ii);
    for jj=1:1e4
        w                 =   (sqrt(sig2)/2)*randn(1,Npi) + (sqrt(sig2)/2)*1i*randn(1,Npi);
        S                 =   R_ms + w;
        [~,tau_est(jj)]   =   max(S);
        tau_est(jj)       =   tau(tau_est(jj));
    end
    tau_hat(ii)     =   mean(tau_est);
    e(ii)           =   sqrt((tau_hat(ii)-delay)^2);        % Root Mean Square Error (RMSE).
end


figure;
plot(CN0,c*sqrt(CRLB));
hold on
plot(CN0,c*e);
title("CN0 comparison with CRLB.");
xlabel("CN0 [dBHz]");
ylabel("RMSE [m]");
legend('CRLB','Signal');
% set(gca, 'YLim', [0,20], 'YTick', 0:2:20,...
%  'YTickLabel', 0:2:20);

The CRLB estimation obtained is

The obtained values are too high.

I would greatly appreciate finding the possible fault.

Thanks in advance.

[Matthieu]

Did you find what was wrong ? I am interested in what would have solved your issue.