how to improve the noise estimation in PSD for non-coherent sampled data?

Question

I am testing an ADC and the signal generator's frequency cannot be programmed to be exact so I am always doing non-coherent sampling (the noise/phase noise is ok, but low in precision). When there is non-coherent sampling, the noise estimation is off a lot with window.

To model this problem , assume an analog signal is sampled by fs, fin = (k+fdelta)/n*fs, where k is an positive integer, fdelta is a value between 0 to 1, n is 2^x and x is a positive integer, as fdelta changes, the estimation of the signal stays pretty accurate but the noise goes off by a lot, like:

Here is this experiment matlab script, where the "goodness of noise estimation" and "goodness of signal estimation" is defined. coswin function is from here . Not sure if I am doing it correctly (maybe I am using the wrong window?)

clear;
format long
% Number of N
n = 32768;
% bin offset
bin_offset = 0:0.1:1;
% fin bin
korig = 271;
% sampling frequency
fs = 1.5e6;
% amplitude
A = 1;
% standard deviation of noise
std_n = A/2^12;
for i = 1:1:length(bin_offset)
k=korig+bin_offset(i);    
% Generate the sinusoid 
t = 0:1/fs:(n-1)/fs;
fin = (k/n)*fs;
y = A*cos(2*pi*fin*t)+std_n*randn(1,length(t));
[sig_pwr,n_pwr] = get_signal_and_noise(y);
% goodness of signal amplitude estimated 
sig_est = (sqrt(sig_pwr)*sqrt(2)/A);
% goodness of noise std estimated
nd_est = n_pwr/(std_n^2);

nn(i) = nd_est;
ss(i) = sig_est;

end
subplot(2,1,1)
plot(bin_offset,nn)
title('goodness of noise estimation')
subplot(2,1,2)
plot(bin_offset,ss)
title('goodness of amplitude estimation')

function [sig_pwr,n_pwr] = get_signal_and_noise(x)
    L=length(x);
    x = x - mean(x);
    Window = coswin(L,4);
    sigL = 3;
    sigR = 3;
    % coherent gain 
    coherent_gain = sum(Window)/ L;
    % equivalent noise bandwidth  
    crct_eqnbw_win = L*sum(Window.^2)/(sum(Window))^2;
    % windowed data
    [sw_r,sw_c] = size(Window);
    if(sw_r>sw_c)
        Window = Window';
    end
    xw = x.*Window;
    % fft
    fft_ret_dsb = 1/L*fft(xw);
    % psd
    psd_dsb = 1/coherent_gain^2/crct_eqnbw_win*(fft_ret_dsb.*conj(fft_ret_dsb));
    % double sided band to single sided band
    ssb_m = floor(L/2)+1;
    psd_ssb = psd_dsb(1:ssb_m);
    psd_ssb(2:end-1) = 2 * psd_ssb(2:end-1) ;
    % signal power and noise power
    [val,idx]=max(psd_ssb);
    sig_pwr = sum(psd_ssb(idx-sigL:idx+sigR));
    n_pwr = sum(psd_ssb(1:idx-sigL-1)) + sum(psd_ssb(idx+sigR+1:end));
end

UPDATE on 2/21/2022:

I found this discrepancy is related to k as well as sigR and sigL. As long as I change sigR and sigL to 4 and makes k large like (> 2000 over n =32768), this is relieved enough to meet my measurement needs. So for k<2000, another noise estimation procedure is derived. But I will still keep this problem open for an elegant solution.

UPDATE on 2/24/2022:

To better explain the discrepancy gets larger as k is lowered, insert the code in get_signal_and_noise to plot the PSD.

NBW = 1/L;
freq = 0:NBW:(L-1)*NBW;
freq_ssb = freq(1:ssb_m);
figure
semilogx(freq_ssb,10*log10(psd_ssb))

when k=100, bin_offset(i)=0.3

when k=1000, bin_offset(i)=0.3

excluding the bin power in vicinity nearby DC could also improve the estimation of noise, which it seems the "low frequency artifcat" is responsible for this discrepancy.

Answer 1

As far as compensating for coherent and non-coherent gain of the window for noise estimation, I think the following line that the OP has should be further multiplied by L/fs; what is originally:

psd_dsb = 1/coherent_gain^2/crct_eqnbw_win*(fft_ret_dsb.*conj(fft_ret_dsb));

should be:

psd_dsb = 1/coherent_gain^2/crct_eqnbw_win* L/fs * (fft_ret_dsb.*conj(fft_ret_dsb));

This would properly express the power density in a per Hz basis and account for the otherwise double-counting from the window's resolution bandwidth as given by crct_eqnbw_win, which is the resolution bandwidth in number of bins. Each bin of the FFT once squared is representing the power in its own bins and the adjacent bins, so dividing by crct_eqnbw_win then provides the power per bin, and then multiplying by L/fs then gives power per unit Hz.

The window used (Blackman-Nuttall) is reasonable if the noise being measured is expected to be higher than the sidelobes of the window. To confirm (as I have preference toward using the Kaiser Window), I plotted the Kernel (Fourier Transform) for a 128 length Blackman-Nuttall below where we see the loss in coherent gain (as the OP computes) and the excellent sidelobe attenuation, possibly surpassing the noise floor being measured that would be within the resolution bandwidth of this window. Further the OP properly avoids any scalloping loss issues in the estimate of the signal power by including all samples over a region in vicinity of the tone being measured.

The OP has further clarified in the comments that the purpose is to measure the signal and noise power in the presence of amplified thermal noise (and if an ADC is involved, I assume the thermal noise has been amplified to be above the quantization noise floor).

The issue lies in the assumption that the sidelobes of the window are below the level of the noise in the samples used to estimate noise. If sigL and sigR are two low (and this will vary based on total window length) then signal power will be included in the noise estimate.

Here we see that for shorter windows the adjacent immediate sidelobe rejection increases. An improved strategy would be to use two different regions for tone estimation and noise estimation. The residual power beyond 3 bins for any tone would be negligible in all cases, so a span of $\pm3$ for signal estimation would be appropriate and would minimize noise in those samples from affecting the signal estimate. However for noise estimation I would recommend a higher number such as $\pm5$ or even $\pm8$ as as the noise being measured is high enough to be above the noise floor of the window beyond this range. We lose the improvement of the estimate based on the ratio of fewer samples (For white noise specifically, the noise estimate improves by $\sqrt{N}$ where $N$ is the total number of samples used. I would imagine for all cases the loss by increasing the exclusion zone to $\pm8$ would be negligible for noise densities above -100 dB/bin).

Additional Details:

The OP mentioned in the comments continued confusion as to why shifting $k$ would make this work. The reason is the noise level chosen is just on the threshold of sensitivity of the window. The plot below specifically shows the window performance for total number of samples $N= 32768$.

The issue is the level of noise being used for measurement, compared to the "noise floor" of this window. The OP is using a standard deviation $\sigma = 1/2^{12}$, and modelled as a Gaussian white noise source (each sample is independent and Gaussian distributed as provided by randn). This total noise power is distributed evenly across the sampled spectrum, so the noise in each bin (without any windowing other than a rectangular window associated with the number of samples used) would be in dB: $20log_{10}(1/2^{12}) - 10log_{10}(N) = -117 $ dB.

The equivalent noise bandwidth of this window (as the OP has properly computed) is 2 bins, which means the measured noise in each bin (without the scaling the OP has provided to compensate for this for a "true-power" measurement) will report the noise equivalent to two bins, or -114 dB. Zooming in on the Kernel below, we see all the areas where the sidelobes of a tone representing signal would fall above this line and thus change a noise measurement if those bins were included. Further, even if the sidelobe is below the line, it will contribute non-coherently to the noise and effect the result (for example a sidelobe 10 dB lower would raise the noise estimate, for that bin, 0.4 dB).

By modifying $k$, the OP is shifting the location of the tone lower in frequency and starting to see the additional interaction of the negative and positive frequency components, further degrading the leakage components in areas where noise is being measured. This effect limits the low frequency range for testing of any real tone, but improving the window under lower noise conditions can help push this lower, as well as increasing the exclusion zone for the noise estimation.

Note that the plot above was generated with a zero padded FFT which serves to interpolate more samples in frequency, and thus shows what we can expect in one plot when the tone is shifted relative to bin center (as the OP had done with fdelta).

The solution is to use the Kaiser window (or Slepian) for an improved algorithm since for that we can dynamically trade the main lobe width for sidelobe rejection and therefore have windowing options available for even higher SNR conditions, and then to modify the exclusion zone to be much wider than the inclusion zone for measuring "signal". As we see in the plots above, the (nice) Blackmann-Nuttall window limits the application to noise densities with around -100 dB/bin (within the window equivalent noise bandwidth), while the OP is using a noise level and window that would result in -114 dB/bin measurements.

I further detail the proper compensation for windowing for noise measurements at this post (but as I mentioned it appears the OP is doing this properly).