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
| 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?
There was a problem hiding this comment.
It's intentional to reduce the number of nans in the resulting values. Makes the test a little more robust
| 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.
There was a problem hiding this comment.
It is yes, I'm trying to make sure all of the options are covered in testing.
| 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.
There was a problem hiding this comment.
Yes, I should do this
| 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?
There was a problem hiding this comment.
Whoa that's a lot. I just gave the most common examples, but we can link to that page
|
Alright @simmsa I mentioned |
Thank you! Nothing else from me. |
# MHKiT-Python v1.1.0 ## Additions * Acoustics: Add millicdecade and WISPR instrument support * Added millidecade spectral conversion * Added a WISPR hydrophone reader * Added a voltage-based `export_audio` resampling option * Refactored band-averaging to avoid losing information at frequency-band boundaries * #447 * Author: @jmcvey3 * Reviewer: @simmsa * DOLfYN: Add [Nortek Aquadopp](https://www.nortekgroup.com/oceanography/aquadopp-series) ADCP support * Added DOLfYN support for reading Aquadopp instruments * Cleaned Nortek parsing code * Simplified handling of the non-cabled ADV orientation flag * #434 * Author: @jmcvey3 * Reviewer: @akeeste * Examples: Add ADCP waves example: `example/adcp_waves_example.ipynb` * Added an example notebook showing how to ingest and analyze wave measurements from a dual-profile Nortek Signature 250 deployment at PacWave. * #430 * Author: @jmcvey3 * Reviewer: @akeeste ## Improvements ### Acoustics - Improved flexibility and robustness of the Acoustics module - Added a configurable FFT length for sound pressure PSDs - Renamed bin/windowing attributes for clarity - Fixed an incorrectly signed gain correction - #433 - Author: @jmcvey3 - Reviewer: @akeeste ### DOLfYN - Refactored PSD calculations to use `scipy.signal.welch` - Replaced DOLfYN's custom Welch-like PSD implementation (built each segment from bin slices, could not overlap FFT segments) with `scipy.signal.welch` - Removed the bin-based segment structure and added a `step` argument to control overlap - Dropped float32 casts in some functions in favor of float64 - Removed `fft.py` and renamed `tools/misc.py` to `tools.py` - Breaking Changes: - NaNs are no longer tolerated in PSD calculations - DOLfYN defaults to 0% overlap to preserve existing array shapes - Acoustics defaults to 50% overlap - PSD output time dimension renamed to `time_psd` - #452 - Author: @jmcvey3 - Reviewer: @akeeste, @simmsa - Improved handling of "averaged" profiles - Fixed handling of Nortek Signature dual-profile ADCP data by defaulting to "_avg" velocity variables when untagged ones are absent. - #430 - Author: @jmcvey3 - Reviewer: @akeeste - Critical PSD bugfix - Fixed a bug where individual FFTs received a 50% overlap twice, which corrupted the first and last spectrum of a timeseries. - #430 - Author: @jmcvey3 - Reviewer: @akeeste ### Examples - Added histograms to ADCP example - #448 - Author: @browniea - Reviewer: @akeeste ### River/IO - Fixed Qhull interpolation and D3D coordinate-system errors - #448 - Author: @browniea - Reviewer: @akeeste - Fixes: #442, #444 - Delft3D module updates - Added new Delft3D coordinate names - Added a new grid-convergence-index calculation function - Added support for xarray/netCDF4 input in the D3D module - #428 - Author: @browniea - Reviewer: @akeeste ### Wave - NDBC Directional Wave Units - Fixed NDBC directional wave spectrum output to return degrees instead of radians - Updated polar plots so 0 deg is at the top and increases clockwise - #437 - Author: @jmcvey3 - Reviewer: @akeeste - Fixes: #427 ### Wave/Hindcast - Added a `hindcast_guard` exception-handling decorator (`hindcast_exceptions.py`) that surfaces a clear error on HSDS request failures, distinguishing the known NLR HSDS outage (#450) from other failures - #449 - Author: @simmsa - Reviewer: @akeeste ## Maintenance - Added Python 3.13 support - #445 - Author: @simmsa - Reviewer: @akeeste - Fixes: #441 - Added [pandas 3](https://pandas.pydata.org/community/blog/pandas-3.0.html) Support - Updated the pandas dependency to allow pandas 3.x - Added a compatibility shim for NDBC missing-value handling related to pandas 3 object to String dtype api changes - Fixed a deprecated period alias in tests - #443 - Author: @simmsa - Reviewer: @akeeste - Fixes: #440 - Updated GitHub Actions CI, expanded installation/developer documentation, refreshed dev environment, and trimmed dependencies - Refactored optional dependencies - Standardized conda/conda-forge environment builds - Scoped black linting to changed files - #436 - Author: @simmsa - Reviewer: @akeeste - Update `rex` dependency to target pypi package to [`NLR-rex[hsds]>=0.5.0`](https://pypi.org/project/NLR-rex/) - #449 - Author: @simmsa - Reviewer: @akeeste
…#452) 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
# MHKiT-Python v1.1.0 ## Additions * Acoustics: Add millicdecade and WISPR instrument support * Added millidecade spectral conversion * Added a WISPR hydrophone reader * Added a voltage-based `export_audio` resampling option * Refactored band-averaging to avoid losing information at frequency-band boundaries * MHKiT-Software#447 * Author: @jmcvey3 * Reviewer: @simmsa * DOLfYN: Add [Nortek Aquadopp](https://www.nortekgroup.com/oceanography/aquadopp-series) ADCP support * Added DOLfYN support for reading Aquadopp instruments * Cleaned Nortek parsing code * Simplified handling of the non-cabled ADV orientation flag * MHKiT-Software#434 * Author: @jmcvey3 * Reviewer: @akeeste * Examples: Add ADCP waves example: `example/adcp_waves_example.ipynb` * Added an example notebook showing how to ingest and analyze wave measurements from a dual-profile Nortek Signature 250 deployment at PacWave. * MHKiT-Software#430 * Author: @jmcvey3 * Reviewer: @akeeste ## Improvements ### Acoustics - Improved flexibility and robustness of the Acoustics module - Added a configurable FFT length for sound pressure PSDs - Renamed bin/windowing attributes for clarity - Fixed an incorrectly signed gain correction - MHKiT-Software#433 - Author: @jmcvey3 - Reviewer: @akeeste ### DOLfYN - Refactored PSD calculations to use `scipy.signal.welch` - Replaced DOLfYN's custom Welch-like PSD implementation (built each segment from bin slices, could not overlap FFT segments) with `scipy.signal.welch` - Removed the bin-based segment structure and added a `step` argument to control overlap - Dropped float32 casts in some functions in favor of float64 - Removed `fft.py` and renamed `tools/misc.py` to `tools.py` - Breaking Changes: - NaNs are no longer tolerated in PSD calculations - DOLfYN defaults to 0% overlap to preserve existing array shapes - Acoustics defaults to 50% overlap - PSD output time dimension renamed to `time_psd` - MHKiT-Software#452 - Author: @jmcvey3 - Reviewer: @akeeste, @simmsa - Improved handling of "averaged" profiles - Fixed handling of Nortek Signature dual-profile ADCP data by defaulting to "_avg" velocity variables when untagged ones are absent. - MHKiT-Software#430 - Author: @jmcvey3 - Reviewer: @akeeste - Critical PSD bugfix - Fixed a bug where individual FFTs received a 50% overlap twice, which corrupted the first and last spectrum of a timeseries. - MHKiT-Software#430 - Author: @jmcvey3 - Reviewer: @akeeste ### Examples - Added histograms to ADCP example - MHKiT-Software#448 - Author: @browniea - Reviewer: @akeeste ### River/IO - Fixed Qhull interpolation and D3D coordinate-system errors - MHKiT-Software#448 - Author: @browniea - Reviewer: @akeeste - Fixes: MHKiT-Software#442, MHKiT-Software#444 - Delft3D module updates - Added new Delft3D coordinate names - Added a new grid-convergence-index calculation function - Added support for xarray/netCDF4 input in the D3D module - MHKiT-Software#428 - Author: @browniea - Reviewer: @akeeste ### Wave - NDBC Directional Wave Units - Fixed NDBC directional wave spectrum output to return degrees instead of radians - Updated polar plots so 0 deg is at the top and increases clockwise - MHKiT-Software#437 - Author: @jmcvey3 - Reviewer: @akeeste - Fixes: MHKiT-Software#427 ### Wave/Hindcast - Added a `hindcast_guard` exception-handling decorator (`hindcast_exceptions.py`) that surfaces a clear error on HSDS request failures, distinguishing the known NLR HSDS outage (MHKiT-Software#450) from other failures - MHKiT-Software#449 - Author: @simmsa - Reviewer: @akeeste ## Maintenance - Added Python 3.13 support - MHKiT-Software#445 - Author: @simmsa - Reviewer: @akeeste - Fixes: MHKiT-Software#441 - Added [pandas 3](https://pandas.pydata.org/community/blog/pandas-3.0.html) Support - Updated the pandas dependency to allow pandas 3.x - Added a compatibility shim for NDBC missing-value handling related to pandas 3 object to String dtype api changes - Fixed a deprecated period alias in tests - MHKiT-Software#443 - Author: @simmsa - Reviewer: @akeeste - Fixes: MHKiT-Software#440 - Updated GitHub Actions CI, expanded installation/developer documentation, refreshed dev environment, and trimmed dependencies - Refactored optional dependencies - Standardized conda/conda-forge environment builds - Scoped black linting to changed files - MHKiT-Software#436 - Author: @simmsa - Reviewer: @akeeste - Update `rex` dependency to target pypi package to [`NLR-rex[hsds]>=0.5.0`](https://pypi.org/project/NLR-rex/) - MHKiT-Software#449 - Author: @simmsa - Reviewer: @akeeste
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
# MHKiT-Python v1.1.0 ## Additions * Acoustics: Add millicdecade and WISPR instrument support * Added millidecade spectral conversion * Added a WISPR hydrophone reader * Added a voltage-based `export_audio` resampling option * Refactored band-averaging to avoid losing information at frequency-band boundaries * #447 * Author: @jmcvey3 * Reviewer: @simmsa * DOLfYN: Add [Nortek Aquadopp](https://www.nortekgroup.com/oceanography/aquadopp-series) ADCP support * Added DOLfYN support for reading Aquadopp instruments * Cleaned Nortek parsing code * Simplified handling of the non-cabled ADV orientation flag * #434 * Author: @jmcvey3 * Reviewer: @akeeste * Examples: Add ADCP waves example: `example/adcp_waves_example.ipynb` * Added an example notebook showing how to ingest and analyze wave measurements from a dual-profile Nortek Signature 250 deployment at PacWave. * #430 * Author: @jmcvey3 * Reviewer: @akeeste ## Improvements ### Acoustics - Improved flexibility and robustness of the Acoustics module - Added a configurable FFT length for sound pressure PSDs - Renamed bin/windowing attributes for clarity - Fixed an incorrectly signed gain correction - #433 - Author: @jmcvey3 - Reviewer: @akeeste ### DOLfYN - Refactored PSD calculations to use `scipy.signal.welch` - Replaced DOLfYN's custom Welch-like PSD implementation (built each segment from bin slices, could not overlap FFT segments) with `scipy.signal.welch` - Removed the bin-based segment structure and added a `step` argument to control overlap - Dropped float32 casts in some functions in favor of float64 - Removed `fft.py` and renamed `tools/misc.py` to `tools.py` - Breaking Changes: - NaNs are no longer tolerated in PSD calculations - DOLfYN defaults to 0% overlap to preserve existing array shapes - Acoustics defaults to 50% overlap - PSD output time dimension renamed to `time_psd` - #452 - Author: @jmcvey3 - Reviewer: @akeeste, @simmsa - Improved handling of "averaged" profiles - Fixed handling of Nortek Signature dual-profile ADCP data by defaulting to "_avg" velocity variables when untagged ones are absent. - #430 - Author: @jmcvey3 - Reviewer: @akeeste - Critical PSD bugfix - Fixed a bug where individual FFTs received a 50% overlap twice, which corrupted the first and last spectrum of a timeseries. - #430 - Author: @jmcvey3 - Reviewer: @akeeste ### Examples - Added histograms to ADCP example - #448 - Author: @browniea - Reviewer: @akeeste ### River/IO - Fixed Qhull interpolation and D3D coordinate-system errors - #448 - Author: @browniea - Reviewer: @akeeste - Fixes: #442, #444 - Delft3D module updates - Added new Delft3D coordinate names - Added a new grid-convergence-index calculation function - Added support for xarray/netCDF4 input in the D3D module - #428 - Author: @browniea - Reviewer: @akeeste ### Wave - NDBC Directional Wave Units - Fixed NDBC directional wave spectrum output to return degrees instead of radians - Updated polar plots so 0 deg is at the top and increases clockwise - #437 - Author: @jmcvey3 - Reviewer: @akeeste - Fixes: #427 ### Wave/Hindcast - Added a `hindcast_guard` exception-handling decorator (`hindcast_exceptions.py`) that surfaces a clear error on HSDS request failures, distinguishing the known NLR HSDS outage (#450) from other failures - #449 - Author: @simmsa - Reviewer: @akeeste ## Maintenance - Added Python 3.13 support - #445 - Author: @simmsa - Reviewer: @akeeste - Fixes: #441 - Added [pandas 3](https://pandas.pydata.org/community/blog/pandas-3.0.html) Support - Updated the pandas dependency to allow pandas 3.x - Added a compatibility shim for NDBC missing-value handling related to pandas 3 object to String dtype api changes - Fixed a deprecated period alias in tests - #443 - Author: @simmsa - Reviewer: @akeeste - Fixes: #440 - Updated GitHub Actions CI, expanded installation/developer documentation, refreshed dev environment, and trimmed dependencies - Refactored optional dependencies - Standardized conda/conda-forge environment builds - Scoped black linting to changed files - #436 - Author: @simmsa - Reviewer: @akeeste - Update `rex` dependency to target pypi package to [`NLR-rex[hsds]>=0.5.0`](https://pypi.org/project/NLR-rex/) - #449 - Author: @simmsa - Reviewer: @akeeste
# MHKiT-Python v1.1.0 ## Additions * Acoustics: Add millicdecade and WISPR instrument support * Added millidecade spectral conversion * Added a WISPR hydrophone reader * Added a voltage-based `export_audio` resampling option * Refactored band-averaging to avoid losing information at frequency-band boundaries * #447 * Author: @jmcvey3 * Reviewer: @simmsa * DOLfYN: Add [Nortek Aquadopp](https://www.nortekgroup.com/oceanography/aquadopp-series) ADCP support * Added DOLfYN support for reading Aquadopp instruments * Cleaned Nortek parsing code * Simplified handling of the non-cabled ADV orientation flag * #434 * Author: @jmcvey3 * Reviewer: @akeeste * Examples: Add ADCP waves example: `example/adcp_waves_example.ipynb` * Added an example notebook showing how to ingest and analyze wave measurements from a dual-profile Nortek Signature 250 deployment at PacWave. * #430 * Author: @jmcvey3 * Reviewer: @akeeste ## Improvements ### Acoustics - Improved flexibility and robustness of the Acoustics module - Added a configurable FFT length for sound pressure PSDs - Renamed bin/windowing attributes for clarity - Fixed an incorrectly signed gain correction - #433 - Author: @jmcvey3 - Reviewer: @akeeste ### DOLfYN - Refactored PSD calculations to use `scipy.signal.welch` - Replaced DOLfYN's custom Welch-like PSD implementation (built each segment from bin slices, could not overlap FFT segments) with `scipy.signal.welch` - Removed the bin-based segment structure and added a `step` argument to control overlap - Dropped float32 casts in some functions in favor of float64 - Removed `fft.py` and renamed `tools/misc.py` to `tools.py` - Breaking Changes: - NaNs are no longer tolerated in PSD calculations - DOLfYN defaults to 0% overlap to preserve existing array shapes - Acoustics defaults to 50% overlap - PSD output time dimension renamed to `time_psd` - #452 - Author: @jmcvey3 - Reviewer: @akeeste, @simmsa - Improved handling of "averaged" profiles - Fixed handling of Nortek Signature dual-profile ADCP data by defaulting to "_avg" velocity variables when untagged ones are absent. - #430 - Author: @jmcvey3 - Reviewer: @akeeste - Critical PSD bugfix - Fixed a bug where individual FFTs received a 50% overlap twice, which corrupted the first and last spectrum of a timeseries. - #430 - Author: @jmcvey3 - Reviewer: @akeeste ### Examples - Added histograms to ADCP example - #448 - Author: @browniea - Reviewer: @akeeste ### River/IO - Fixed Qhull interpolation and D3D coordinate-system errors - #448 - Author: @browniea - Reviewer: @akeeste - Fixes: #442, #444 - Delft3D module updates - Added new Delft3D coordinate names - Added a new grid-convergence-index calculation function - Added support for xarray/netCDF4 input in the D3D module - #428 - Author: @browniea - Reviewer: @akeeste ### Wave - NDBC Directional Wave Units - Fixed NDBC directional wave spectrum output to return degrees instead of radians - Updated polar plots so 0 deg is at the top and increases clockwise - #437 - Author: @jmcvey3 - Reviewer: @akeeste - Fixes: #427 ### Wave/Hindcast - Added a `hindcast_guard` exception-handling decorator (`hindcast_exceptions.py`) that surfaces a clear error on HSDS request failures, distinguishing the known NLR HSDS outage (#450) from other failures - #449 - Author: @simmsa - Reviewer: @akeeste ## Maintenance - Added Python 3.13 support - #445 - Author: @simmsa - Reviewer: @akeeste - Fixes: #441 - Added [pandas 3](https://pandas.pydata.org/community/blog/pandas-3.0.html) Support - Updated the pandas dependency to allow pandas 3.x - Added a compatibility shim for NDBC missing-value handling related to pandas 3 object to String dtype api changes - Fixed a deprecated period alias in tests - #443 - Author: @simmsa - Reviewer: @akeeste - Fixes: #440 - Updated GitHub Actions CI, expanded installation/developer documentation, refreshed dev environment, and trimmed dependencies - Refactored optional dependencies - Standardized conda/conda-forge environment builds - Scoped black linting to changed files - #436 - Author: @simmsa - Reviewer: @akeeste - Update `rex` dependency to target pypi package to [`NLR-rex[hsds]>=0.5.0`](https://pypi.org/project/NLR-rex/) - #449 - Author: @simmsa - Reviewer: @akeeste
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