Refactor DOLfYN PSD Scripts to use Scipy.Signal.Welch#452
Conversation
…y vector weighting for calculating band-averaged SPLs
…leanup deprecation warning strings
…ed, fix missing step input, set acoustics default to use 50% overlap
…ons to interpolate to the same timestamps (this does intro nan's, unlike previously when assuming timestamps were the same)
…d time size is the same as the other variables
jmcvey3
added
bug
|
Apologies for a mess of repeat commits from a previous PR, and I'll update the tests when I have another moment |
…and clean up power_spectral_density function
| dat["vel_b5"].isel(range_b5=len(dat["range_b5"]) // 2), | ||
| freq_units="Hz", | ||
| noise=0.01, | ||
| noise=0.0001, |
There was a problem hiding this comment.
Previous test used 0.01. Is it intentional that this change uses 0.001, an extra decimal?
| tdat["stress_demean"] = bnr.reynolds_stress(dat["vel"], detrend=False) | ||
| tdat["csd"] = bnr.cross_spectral_density( | ||
| dat["vel"], freq_units="rad", window="hamming", n_fft_coh=10 | ||
| dat["vel"], freq_units="rad", window="hann", n_fft_coh=10 |
There was a problem hiding this comment.
Is the change from hamming to hann window intentional? Mostly want to make sure this change is not a typo.
| n_bin=None, | ||
| n_fft=None, | ||
| step=None, | ||
| pct_overlap=0.5, |
There was a problem hiding this comment.
Can we put in a check in this method and also the _csd_base method to guard against obvious out of spec user inputs to the pct_overlap parameter, i.e. 50.0.
if pct_overlap < 0 or pct_overlap > 1:
raise ValueError(f"pct_overlap must be between 0 and 1, got {pct_overlap}.")
The specification is spelled out in the docstring, but I think we should try to prevent any out of spec values from passing into scipy through this line:
noverlap=int(pct_overlap * n_fft)
If the user does input an out of spec value it is likely that signal.welch will produce an error if noverlap > nperseg:
But the error will likely be confusing. I think it would be more correct catch this on the MHKiT side.
| Type of window to apply to each FFT segment. | ||
| Options ('boxcar', 'hann', 'hamming', 'blackman', 'bartlett'). | ||
| Default: 'hann'. | ||
| noise : numeric or array |
There was a problem hiding this comment.
Is it relevant to link to the available window functions in scipy here:
https://docs.scipy.org/doc/scipy/reference/signal.windows.html#module-scipy.signal.windows
Or are we only supporting a subset of windowing functions?
3 participants
Big refactor I'd been intending to finish for a while.
This refactor solves a DOLfYN structural problem between its binning architecture and the traditional sliding-window architecture of Welch's algorithm. Because DOLfYN was using slices of bins to compute each PSD segment, it was impossible to ensure data overlap between one bin and the next, meaning the industry-standard 50% overlap cannot be accomplished. This PR removes the bin structure from the PSD computation and ensures overlap using a "step" argument, which is also now built into the bin structure to properly average dimensions for the PSD.
In doing this, DOLfYN's custom PSD code, which was a custom implementation of Welch's algorithm without the ability to overlap FFT segments, has been replaced with scipy.signal.welch. This has one breaking change, in that NaN's are no longer tolerated in the PSD calculation.
To avoid breaking changes, the DOLfYN code will default to an overlap of 0%. This ensures array shapes will remain the same. On the other hand, since it's new, the Acoustics module will default to 50%.
The second breaking change is that the PSD code now renames the input time dimension to "time_psd". In this way, if a different overlap is used, and the PSD is saved into dataset with an averaged "time" dimension, the PSD will not follow the xarray default and go to nan.
Finally, I removed the float32 datatypes in certain functions since sometimes we do need float64 precision. Removing custom code means fft.py no longer exists, and tools/misc.py was renamed to tools.py