calc_wn is biased

[ykyohei]

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 2*n degrees of freedom, thus correction factor will be 2*n/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.