simonsobs / sotodlib

Simons Observatory: Time-Ordered Data processing library.
MIT License
16 stars 19 forks source link

calc_wn is biased #995

Open ykyohei opened 1 month ago

ykyohei commented 1 month ago

fft_ops.calc_wn takes median of psd as a estimator of white noise level (wn), and this is biased from the average of psd.

The power spectrum density is usually calculated by the fft_ops.calc_psd using scipy.welch, and the amount of bias can be different between different observations because we allow different nperseg for different observations.

There are at least two ways to obtain the unbiased estimate of wn.

  1. Take an average of the psd in the frequency range without peaks.
  2. Take a median of the psd to be robust against peaks and then debias. For example. we set noverlap = 0 in scipy welch and debias depending on the number of segments. The average PSD of non-overlapping n segments (scipy.welch psd with noverlap = 0) follows chi square distribution with 2n degrees of freedom, thus correction factor will be `2n/scipy.stats.chi2.ppf(0.5, 2*n)`. The maximum correction factor for psd is ~1.44 when n=1. The current fft_ops.calc_psd allows overlap between segments and it's more complicated to obtain the correction factor.

Controlling the bias of wn is important for correct understanding of the detector noise performance, and for correct inverse variance weighting for mapmaking.

Figure 25 (3)