[Feature Req.] Decide on default filtering for real-space images.

[nikos-triantafyllou]

For tuesday.core.intrument_models.noise .

When adding noise in UV space, the transform back to real space for imaging is not obvious.

Below I show an image where different decisions can be made about the filtering.

Upper and lower rows are different projections of a cubic chunk of a lightcone (y-x and y-z respectively).

First column is the intrinsic signal form one simulation at z=8 for default parameters (RMSE of the signal is ~5) for 300Mpc side and 200 pixels.
The values of the first column are divided by the standard deviation of the signal. The rest of the columns are standardized.

The rest of the columns show different schemes for filtering.

• 2nd column: no filtering, stadard fft from UV to image space.
• 3rd column: no filtering, stadard fft from UV to image space but only considering baselines<1KM.
• 4th column: inverse variance weighting
• 5th column: wiener filtering. Note that this depend on the power spectrum of the signal, which is inherently not known in real life but results in an almost perfect visually revovery of the signal.

To reproduce:

add this code chunck in the end of the tuesday.core.instrument_models.noise.observe_lightcone function:

    if spatial_taper is not None:
        print(lc_uv_nu.shape)
        _, nx,  ny, _  = lc_uv_nu.shape
        window_fnc = taper2d(lightcone.shape[0], spatial_taper)[:, -ny:]
        window_fnc = np.fft.fftshift(
            window_fnc, axes=(0,)
        )  # shift the window to be in the right format for FFT
        lc_uv_nu *= window_fnc[None, ..., None]
        
    if apply_inverse_variance_weighting:
        with np.errstate(divide="ignore", invalid="ignore"):
            w = 1.0 / (thermal_rms_uv**2).value   # inverse variance weights
        w[np.isinf(w)] = 0.0                      # set unsampled cells to 0, w is in uv_nu space
        wsum = w.sum(axis=(0,1), keepdims=True)
        lc_uv_nu *= w.shape[0] * w.shape[1] * w[None, ...] / wsum[None, ...] # 2nd and 3rd dimensions of w and wsum will not match but numpy is broadcasting
        
    if apply_wiener:
        # Estimate signal power spectrum from the lightcone itself
        # lc_uv_nu is already in uv space, use it to estimate P_signal
        # Average over realizations and frequency axis to get a stable estimate
        p_signal = np.abs(np.fft.rfft2(lightcone.value, axes=(0,1)))**2  # shape (nx, ny_rfft)
        # Smooth it to get a stable estimate (optional but recommended)
        # p_signal = gaussian_filter(p_signal, sigma=2)  # from scipy.ndimage

        p_noise = (thermal_rms_uv**2).value  # your existing noise variance

        with np.errstate(divide="ignore", invalid="ignore"):
            w = p_signal / (p_signal + p_noise)   # Wiener filter, values in [0,1]

        w[np.isnan(w)] = 0.0
        w[thermal_rms_uv.value == 0] = 0.0  # zero where unsampled

        wsum = w.sum(axis=(0,1), keepdims=True)
        lc_uv_nu *= w.shape[0] * w.shape[1] * w[None, ...] / wsum[None, ...]

    noisy_lc_real = np.fft.irfft2(lc_uv_nu, s=(nx, nx), axes=(1, 2)).to(lightcone.unit)
    

    return noisy_lc_real, lightcone_redshifts```

[steven-murray]

This is a very nicely illustrated issue, thanks @nikos-triantafyllou.

Yes, you're absolutely right that there is not one unique way to go from UV to image space. The "best" way depends on what you're trying to do with the image in the end, and how you want your noise statistics to correlate.

Given that there is no "correct" way, I think the best thing we can do is to provide the user with options. I think the default should be to do as we are doing now (i.e. not applying any filtering beyond the uv coverage). But I think we should make the choice prominent in the docs.

I would probably make a separate function for each of the UV filters to make this a little more elegant.

The "<1km" filter is fine, but of course can already be achieved by simply passing in an Observatory that has put a filter on its baselines. NEvertheless, changing that filter requires regridding, whereas just applying a UV filter does not, so perhaps it's wise to keep that option in there.