4
$\begingroup$

I am trying to understand the concept of CRLB for a simple sinusoid in AWGN with a simulation in MATLAB.

According to Kay (p. 57) the CRLB for phase estimation is:

$var\{\phi\}=\frac{2(2N-1)}{SNR\cdot N(N+1)}$.

I generated a MATALB script that simulates 200 sinusoid signals with different amplitudes and fixed noise power. The phase of the signals is calculated through a FFT. Then, the standard deviation is calculated for every amplitude and compared to the CRLB.

Unfortunately the results show, that the estimated phase has a lower standard deviation than the CRLB:

enter image description here

I can't figure out why and what the problem with my code is. Here is the code I use:

%% Init
clear all;close all;clc;

f1 = 500e3;%signal frequency
phi = 58*pi/180;%signal phase
T = 1e-3;%signal duration

fs = 4e6;%sampling frequency
t = 0:1/fs:T-1/fs;%time vector
nfft = 2*2^nextpow2(fs*T);%fft points
f = 0:fs/nfft:fs/2-fs/nfft;%frequency vector
Samples = round(fs*T);%number of samples
NumSigs = 200;%number of signals
Win = hanning(Samples);%window function
WinMat = repmat(Win, 1, NumSigs);%window matrix for fft
DataTime = zeros(Samples,NumSigs);%time samples of signals

sigmaN = 0.02;%rms noise amplitude
A = logspace(-2,0,15);%signal amplitudes
nAmps = length(A);%number of signal amplitudes

for a = 1:nAmps%for all signal amplitudes
    %% Signal generation
    for i = 1:NumSigs
        DataTime(:,i) = A(a)*cos(2*pi*f1.*t + phi) + sigmaN*randn(Samples,1)';%generate signals   
    end

    %% FFT
    DataTimeWin = WinMat.*DataTime;%apply window
    DataSpec = (1/sum(Win))*fft(DataTimeWin, nfft, 1);%fft
    DataPow = (2*abs(DataSpec)/sqrt(2)).^2  / (50*1e-3);%power in mW, 50 Ohm system
    DataPhase = atan2(imag(DataSpec),real(DataSpec));%phase

%     figure('units','normalized','outerposition',[0 0 1 1]);%plot
%     subplot(2,1,1)
%     plot(10*log10(DataPow(1:nfft/2,:)));%plot all signals
%     grid on; grid minor; 
%     subplot(2,1,2)
%     plot(DataPhase(1:nfft/2,:));%plot phase values
%     grid on; grid minor;  

    %% Signal Peak Search
    for j = 1:NumSigs
        [~, locs(j)] = findpeaks(DataPow(:,j),'SortStr', 'descend', 'NPeaks', 1);%find signal peak
        PeakPhase(j) = DataPhase(locs(j),j);%get phase value
    end

    %% CRLB frequency
    SNR(a) = 10*log10( (A(a)/sqrt(2))^2 ) - 10*log10( (sigmaN^2) );%calculate SNR

    SIGMAp(a) = std(PeakPhase);%standard deviation of phase

    CRLBp(a) = (2*(2*Samples-1))/(10^(SNR(a)/10)*Samples*(Samples+1));%CRLB for phase

end

%% Plot
figure('units','normalized','outerposition',[0 0 1 1]);
plot(SNR,10*log10(SIGMAp),SNR,10*log10(sqrt(CRLBp)),'--');
grid on; grid minor;
xlabel('SNR (dB)');
ylabel('Standard deviation (dBrad)');
title('Phase Est.');
legend('SD. Phase', 'CRLB Phase');

Since I use the formula of Kay for CRLB calculation and the obvious way to calculate the SNR, I assume there is something wrong with how I estimate the standard deviation of the phase. Can anyone give me a hint as to what I am doing wrong?

$\endgroup$

1 Answer 1

5
$\begingroup$

Frequency estimation is fraught with peril. I believe that one issue is that the maximum likelihood frequency estimator gets much better accuracy than just choosing the FFT bin maximizer as in the code.

You can see this from the second line of (3.41) in Kay: $$ {\tt var}(\hat{f}_0) \ge \frac{12}{(2\pi)^2 \eta N(N^2-1)} $$ whereas the FFT variance of the frequency estimate only goes as $O\left(\frac{1}{N^2}\right)$.

Because that estimator is poor, the bin being chosen for the phase is wrong, so of course the variance of the phase is out of whack with the bound.

One way to try to get around it using the existing code is to change

nfft = 2*2^nextpow2(fs*T);%fft points

to

nfft = 2^18;

so that there are many closer bins than in the original.

If I do that one change to the code then I get the following diagram.

Revised plot.

While it doesn't completely solve the problem, it does show that over part of the range the bound is actually a lower bound.

If a better frequency estimator were used, then the statistics should align.

Why does the poorer frequency estimator yield "better" phase estimation? This is likely because the better phase estimator is biased: the frequency bins are closer to the actual true frequency than the real peak of the noisy spectrum.

A stopped clock is right twice a day.

$\endgroup$
6
  • 2
    $\begingroup$ It's particularly suspicious to me that the estimate is very close to 3 dB better than the bound, such as we would get when we inadvertently take the real part of a complex noise process. $\endgroup$ Dec 31, 2023 at 14:09
  • 1
    $\begingroup$ @DanBoschen Yes, there is lots else to be improved with the way this was approached... but I decided to start somewhere. :-) $\endgroup$
    – Peter K.
    Dec 31, 2023 at 17:06
  • 1
    $\begingroup$ Oh it was just a quick observation and I hope you didn’t take it as a criticism to your approach; I had +1’d the insightful answer. Happy New Year Peter! $\endgroup$ Jan 1 at 14:16
  • $\begingroup$ @DanBoschen I took your comment as you intended it. It made me realize there are several other issues too. Happy New Year to you too! $\endgroup$
    – Peter K.
    Jan 1 at 16:00
  • $\begingroup$ @PeterK. Thank you for your answer. Very helpful! You mentioned that there is more room for improvement. Can you please give me a hint about what I could have done differently and how to do that? I tried to make the example as simple as possible to get an understanding of this topic. $\endgroup$
    – HaLoe
    Jan 1 at 20:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.