Fix ssp correction

[krisvanneste]

Claudio,

Here is my pull request to fix the interpolation issue when residual correction is applied. As you suggested, I created the branch from the fix_source_residuals branch, but now it says that it can't automatically merge... Let me know if I need to do it differently.

I ran the fix with my test case and it appears to work, but you should be aware that in the spectrum returned by the station_correction function, the linear-spaced and log-spaced versions may not have the same length and frequency range. The linear one is truncated to the common frequency range with the station residual, while the logarithmic one keeps the original length, with values outside the common frequency range set to NaN. I'm not sure if this could cause any problems further in the code, but I did not notice anything.

[krisvanneste]

Note that I won't be in the office tomorrow. We can continue on Thursday or Friday.

[claudiodsf]

Dear Kris, thanks for the PR.

There were indeed many conflicts (my suggestion to start from the previous branch was bad).
To fix that, I started a new branch from current upstream/main and "cherry-picked" the three commits, then renamed the branch to fix_ssp_correction and force-pushed.
Please force-pull on your side.

I can now check the PR contents... 😉

[claudiodsf]

I added a couple of commits to move common interpolation code to a function and to remove the unnecessary return spec_st (spec_st is altered in place)

[krisvanneste]

I force-pulled the branch.

[claudiodsf]

There is still a problem in my tests when computing spectral weights from S/N: signal spectrum and noise spectrum sometimes do not have the same length.

Today's is holiday in France and tomorrow I'll take the day off.
I'll be back working on this PR on Monday 😉

[krisvanneste]

On my side it still works.
Are the spectral weights computed before or after the residual correction?

[claudiodsf]

I think the problem is here:

specnoise is computed independently of spec, based on the assumption that both are computed from a time window that has the same length, and therefore they must have the same frequencies.

This assumption is no longer valid if the frequency range spec is modified during the correction.

A possible solution would be to pass the frequency range as an extra optional argument when calling _build_spectrum() for spec_noise, and interpolate spec_noise to this new frequency range.

So something like:

specnoise = _build_spectrum(config, trace_noise, interpolate_freqs=spec.freqs)

[krisvanneste]

OK, I will have a look as well. Enjoy your long weekend!

[krisvanneste]

I understand the problem now. I think it would be better not to change the frequency range of the signal spectrum, and set NaN values outside the frequency range of the residual for the linear-spaced spectrum as well, thus also keeping linear-spaced and log-spaced versions in sync. We then need to check if the NaN values do not give any problem when calculating spectral SNR. I will improve the ssp_correction.station_correction function.

[krisvanneste]

Claudio, I changed the station correction function in a way that the original frequency range is not modified. Can you check if this solves your problem with noise weighting?

[claudiodsf]

Hi Kris, thanks for this additional work.

I had to make some further modifications:

• replace NaNs in the weight spectrum with small values: otherwise the logspaced spectrum becomes all NaNs
• replace np.min() and np.max() with np.nanmin() and np.nanmax() in stacked spectral plotting

Now it seems to work, but I need to test it more extensively.

Could you also check on your side?

[krisvanneste]

OK

[claudiodsf]

One more commit related to NaN data 😉

[krisvanneste]

My results look identical.
Note that you can't use scipy.interpolate.interp1d with arrays containing NaN values. I had to ignore these to interpolate the log-spaced spectrum in station_correction.

[claudiodsf]

My results look identical.

Great!

Note that you can't use scipy.interpolate.interp1d with arrays containing NaN values. I had to ignore these to interpolate the log-spaced spectrum in station_correction.

Ok, so maybe I can use a similar approach for weight computation, instead of replacing NaN with 1e-9...

[krisvanneste]

Ok, so maybe I can use a similar approach for weight computation, instead of replacing NaN with 1e-9...

I did it like this:

nan_idxs = np.isnan(spec_corr.data_mag)
f = interp1d(spec_corr.freq[~nan_idxs], spec_corr.data_mag[~nan_idxs],
             bounds_error=False)
spec_corr.data_mag_logspaced = f(spec_corr.freq_logspaced)

This is possible because we know the NaN values only occur at the start and/or end. The interpolated spectrum will also be NaN where the original one is NaN, because we don't allow extrapolation.