Gaussian noise returns garbage below PSD fmin

@nikhil.sarin and I discovered that below the minimum frequency in the PSD file passed to gaussian_noise, the function returns garbage. The safest fix is to just cut these extra frequencies so that the minimum frequency returned corresponds to the PSD fmin.

To be clear, this is only an issue if we keep the arbitrarily low frequency noise. Even if we do, I think this shouldn't have any effect on pe. It'll just change the noise evidence.

I'm not sure if we currently do this, but we should be more consistent with our use of fmin and fmax. This is a good reminder that we should set the data below fmin to be zero.

I'm happy to deal with this.

I was just going to add a cut inside the function so that it only returns data for frequencies greater than the lowest frequency in the PSD text file. It might not affect PE, but if we're no longer maintaining MGWT separately, then this is an easy fix to avoid other issues (i.e. with making frames).

Actually I think the better way to do this is to allow the user to specify what fmin they want, make this an argument of the function, and return data all the way down to zero but so that it plateaus below fmin. This will work well if there's currently a method for specifying fmin in the likelihood (I don't know if there is, but if not there should be). The ASD calculated from time domain gaussian noise data gets funky if you set f<fmin=0 or cut it out altogether. Thoughts?

I would prefer cutting it out below fmin instead of user specifying fmin. Don't have a strong opinion on this though

What would be the fmin that you'd cut below? For example, I made a time series after cutting the output of gaussian_noise at fmin=10Hz, then from this time series I calculated the ASD, and the result still has spectral content below 10Hz and it biases the neighboring frequencies. This seems dangerous to me. The same thing happened when I set the values of f<fmin to 0. First graph is zeroing, second is cutting

hmmm alright, zeroing might be the way to go then. But it looks like its doing a pretty poor job interpolating outside the noise curve so that needs to be checked. Instead of cutting the output can you set the extrapolated part of the noise curve to be 0 instead?

It's not actually interpolating anything here. This is the ASD calculated from the time domain data. So even when I've zeroed the raw freq domain data from gaussian_noise once I nifft into the time domain and calculate the ASD from the time domain data, there is still spectral content at f<fmin when it should either not exist (cutting, second graph) or be set to zero (first graph). Which is why I'm advocating that the safest thing is to let the data plateau for frequencies less than some user specified frequency. Which is essentially what the function currently does, except instead of taking fmin as an argument, it just uses the lowest frequency in the PSD file. But you'll get biased evidences if the fmin in your data doesn't match the fmin in the PE likelihood, so that's why I think fmin should be an argument of the function

Should definitely be interpolating...

Ah I mean to make these plots I am not interpolating. The function is interpolating, sure. But these are the ASDs estimated from time-domain data, not the raw frequency domain data from MGWT.

Yeah so I think this problem might be mitigated significantly if in the function we interpolate with a different spline/or different fill value, rather than the interp1d default

This problem isn't coming from the function though. The function output is fine and does what it's supposed to do. I think this is coming from spectral leakage from the FFT or from windowing or some other signal processing thing

Surely an FFT should not have any spectral leaking? Idk

changed milestone to %0.1

I've set this as 0.1 as it seems important. What is the status, is it possible to solve this soon or should we tag as 'future'? Sorry I don't really understand the problem hence the (probably ignorant) questions.

It's not really an issue for PE, but if we're not maintaining MonashGWTools independently anymore, then this refers more to the gaussian_noise function's other applications. @colm.talbot wanted to set the strain for all frequencies below fmin to 0, while the default is to set them all to h(fmin). I don't really have a preference between the two, but I'd like fmin to be an argument of the function, so the user knows what fmin they're getting automatically instead of having to implement it by hand afterwards. We never really reached a consensus...

Solution: change extrapolation to NANs rather than max of psd. Need to make sure this doesn't break the kode. Do this by merge, please!!

changed milestone to %Future

I just pushed a commit to https://git.ligo.org/Monash/tupak/merge_requests/79.

NANs were causing issues so now:

• the interpolated PSD has fill_value=np.inf outside the specified frequency range, meaning that the PSD is infinite outside the range.
• get_noise_realisation returns 0 outside the frequency range the PSD is known for.

closed

unassigned @sylvia.biscoveanu