# FFT-based phase estimation better than CRLB in MATLAB simulation

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:

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)');
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?

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.

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.

• 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. Dec 31, 2023 at 14:09
• @DanBoschen Yes, there is lots else to be improved with the way this was approached... but I decided to start somewhere. :-)
– Peter K.
Dec 31, 2023 at 17:06
• 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! Jan 1 at 14:16
• @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!
– Peter K.
Jan 1 at 16:00
• @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. Jan 1 at 20:35