Skip to content

Commit

Permalink
Ensure that noise subtraction is working properly
Browse files Browse the repository at this point in the history
  • Loading branch information
@jmcvey3
jmcvey3 committed Jan 12, 2024
1 parent 16b4b37 commit d0ddc77
Show file tree
Hide file tree
Showing 5 changed files with 72 additions and 24 deletions.
19 changes: 17 additions & 2 deletions dolfyn/adp/turbulence.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -623,7 +623,7 @@ def check_turbulence_cascade_slope(self, psd, freq_range=[0.2, 0.4]):

return m, b

def calc_dissipation_LT83(self, psd, U_mag, freq_range=[0.2, 0.4]):
def calc_dissipation_LT83(self, psd, U_mag, freq_range=[0.2, 0.4], noise=None):
"""Calculate the TKE dissipation rate from the velocity spectra.
Parameters
Expand All @@ -633,8 +633,11 @@ def calc_dissipation_LT83(self, psd, U_mag, freq_range=[0.2, 0.4]):
U_mag : xarray.DataArray (time)
The bin-averaged horizontal velocity (a.k.a. speed) from a single depth bin (range)
f_range : iterable(2)
The range over which to integrate/average the spectrum, in units
The range over which to integrate/average the spectrum, in units
of the psd frequency vector (Hz or rad/s)
noise : float or array-like
A vector of the noise levels of the velocity data with
the same first dimension as the velocity vector (time).
Returns
-------
Expand Down Expand Up @@ -669,6 +672,18 @@ def calc_dissipation_LT83(self, psd, U_mag, freq_range=[0.2, 0.4]):
raise Exception('PSD should be 2-dimensional (time, frequency)')
if len(U_mag.shape) != 1:
raise Exception('U_mag should be 1-dimensional (time)')
if noise is not None:
if np.shape(noise)[0] != np.shape(psd)[0]:
raise Exception(
'Noise should have same first dimension as PSD')
else:
noise = np.array(0)

# Noise subtraction from binner.TimeBinner.calc_psd_base
psd = psd.copy()
if noise is not None:
psd -= noise**2 / (self.fs / 2)
psd = psd.where(psd > 0, np.min(np.abs(psd)) / 100)

freq = psd.freq
idx = np.where((freq_range[0] < freq) & (freq < freq_range[1]))
Expand Down
23 changes: 20 additions & 3 deletions dolfyn/adv/turbulence.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -165,7 +165,7 @@ def calc_csd(self, veldat,
'time': time,
'coh_freq': coh_freq},
dims=['C', 'time', 'coh_freq'],
attrs={'units': units,
attrs={'units': units,
'n_fft_coh': n_fft,
'long_name': 'Cross Spectral Density'})
csd['coh_freq'].attrs['units'] = freq_units
Expand Down Expand Up @@ -296,7 +296,7 @@ def check_turbulence_cascade_slope(self, psd, freq_range=[6.28, 12.57]):

return m, b

def calc_epsilon_LT83(self, psd, U_mag, freq_range=[6.28, 12.57]):
def calc_epsilon_LT83(self, psd, U_mag, freq_range=[6.28, 12.57], noise=None):
"""Calculate the dissipation rate from the PSD
Parameters
Expand All @@ -308,6 +308,9 @@ def calc_epsilon_LT83(self, psd, U_mag, freq_range=[6.28, 12.57]):
freq_range : iterable(2) (default: [6.28, 12.57])
The range over which to integrate/average the spectrum, in units
of the psd frequency vector (Hz or rad/s)
noise : float or array-like
A vector of the noise levels of the velocity data with
the same first dimension as the velocity vector.
Returns
-------
Expand Down Expand Up @@ -339,8 +342,22 @@ def calc_epsilon_LT83(self, psd, U_mag, freq_range=[6.28, 12.57]):
"""

# Ensure time has been averaged
if len(psd.time)!=len(U_mag.time):
if len(psd.time) != len(U_mag.time):
raise Exception("`U_mag` should be from ensembled-averaged dataset")
if noise is not None:
if np.shape(noise)[0] != 3:
raise Exception(
'Noise should have same first dimension as velocity')
# Edit first dim to work with PSD
noise = noise.assign_coords({'dir': psd['S'].values}).rename({'dir': "S"})
else:
noise = np.array([0, 0, 0])[:, None, None]

# Noise subtraction from binner.TimeBinner.calc_psd_base
psd = psd.copy()
if noise is not None:
psd -= (noise**2 / (self.fs / 2))
psd = psd.where(psd > 0, np.min(np.abs(psd)) / 100)

freq = psd.freq
idx = np.where((freq_range[0] < freq) & (freq < freq_range[1]))
Expand Down
2 changes: 1 addition & 1 deletion dolfyn/binned.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -401,7 +401,7 @@ def calc_psd_base(self, dat, fs=None, window='hann', noise=0,
for slc in slice1d_along_axis(dat.shape, -1):
out[slc] = psd(dat[slc], n_fft, fs,
window=window, step=step)
if noise != 0:
if np.any(noise):
out -= noise**2 / (fs/2)
# Make sure all values of the PSD are >0 (but still small):
out[out < 0] = np.min(np.abs(out)) / 100
Expand Down
18 changes: 16 additions & 2 deletions dolfyn/tests/test_analysis.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -114,14 +114,16 @@ def test_adv_turbulence(make_data=False):

tdat['stress_detrend'] = bnr.calc_stress(dat['vel'])
tdat['stress_demean'] = bnr.calc_stress(dat['vel'], detrend=False)
tdat['csd'] = bnr.calc_csd(
dat['vel'], freq_units='rad', window='hamm', n_fft_coh=10)
tdat['csd'] = bnr.calc_csd(dat['vel'], freq_units='rad', window='hamm', n_fft_coh=10)
tdat['LT83'] = bnr.calc_epsilon_LT83(tdat['psd'], tdat.velds.U_mag)
tdat['noise'] = bnr.calc_doppler_noise(tdat['psd'], pct_fN=0.8)
tdat['LT83_noise'] = bnr.calc_epsilon_LT83(tdat['psd'], tdat.velds.U_mag, noise=tdat['noise'])
tdat['SF'] = bnr.calc_epsilon_SF(dat['vel'][0], tdat.velds.U_mag)
tdat['TE01'] = bnr.calc_epsilon_TE01(dat, tdat)
tdat['L'] = bnr.calc_L_int(acov, tdat.velds.U_mag)
slope_check = bnr.check_turbulence_cascade_slope(
tdat['psd'][-1].mean('time'), freq_range=[10, 100])
tdat['psd_noise'] = bnr.calc_psd(dat['vel'], freq_units='rad', noise=[0.06, 0.04, 0.01])

if make_data:
save(tdat, 'vector_data01_bin.nc')
Expand All @@ -144,6 +146,7 @@ def test_adcp_turbulence(make_data=False):
tdat['tau2'] = bnr.calc_shear2(tdat.vel)
tdat['I'] = tdat.velds.I
tdat['ti'] = bnr.calc_ti(dat.velds.U_mag, detrend=False)
#dat.velds.rotate2('beam')
tdat['psd'] = bnr.calc_psd(dat['vel'].isel(
dir=2, range=len(dat.range)//2), freq_units='Hz')
tdat['noise'] = bnr.calc_doppler_noise(tdat['psd'], pct_fN=0.8)
Expand All @@ -157,16 +160,27 @@ def test_adcp_turbulence(make_data=False):
tdat['wpwp'] = bnr.calc_tke(dat['vel_b5'], noise=tdat['noise'])
tdat['dissipation_rate_LT83'] = bnr.calc_dissipation_LT83(
tdat['psd'], tdat.velds.U_mag.isel(range=len(dat.range)//2), freq_range=[0.2, 0.4])
tdat['dissipation_rate_LT83_noise'] = bnr.calc_dissipation_LT83(
tdat['psd'], tdat.velds.U_mag.isel(range=len(dat.range)//2), freq_range=[0.2, 0.4], noise=tdat['noise'])
tdat['dissipation_rate_SF'], tdat['noise_SF'], tdat['D_SF'] = bnr.calc_dissipation_SF(
dat.vel.isel(dir=2), r_range=[1, 5])
tdat['friction_vel'] = bnr.calc_ustar_fit(
tdat, upwp_=tdat['stress_vec5'].sel(tau='upwp_'), z_inds=slice(1, 5), H=50)
slope_check = bnr.check_turbulence_cascade_slope(
tdat['psd'].mean('time'), freq_range=[0.4, 4])
tdat['psd_noise'] = bnr.calc_psd(dat['vel'].isel(
dir=2, range=len(dat.range)//2), freq_units='Hz', noise=0.01)

with pytest.raises(Exception):
bnr.calc_psd(dat['vel'], freq_units='Hz', noise=0.01)

if make_data:
save(tdat, 'Sig1000_tidal_bin.nc')
return

with pytest.raises(Exception):
bnr.calc_psd(dat['vel'], freq_units='Hz', noise=0.01)
with pytest.raises(Exception):
bnr.calc_psd(dat['vel'][0], freq_units='Hz', noise=0.01)
assert np.round(slope_check[0].values, 4), -1.0682
assert_allclose(tdat, load('Sig1000_tidal_bin.nc'), atol=1e-6)
34 changes: 18 additions & 16 deletions dolfyn/velocity.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -963,7 +963,8 @@ def calc_ti(self, U_mag, noise=0, thresh=0, detrend=False):
coords=coords,
dims=dims,
attrs={'units': '% [0,1]',
'long_name': 'Turbulence Intensity'})
'long_name': 'Turbulence Intensity',
'comment': f'TI was corrected from a noise level of {noise} m/s'})

def calc_tke(self, veldat, noise=None, detrend=True):
"""Calculate the turbulent kinetic energy (TKE) (variances
Expand Down Expand Up @@ -1055,16 +1056,15 @@ def calc_psd(self, veldat,
veldat : xr.DataArray
The raw velocity data (of dims 'dir' and 'time').
freq_units : string
Frequency units of the returned spectra in either Hz or rad/s
Frequency units of the returned spectra in either Hz or rad/s
(`f` or :math:`\\omega`)
fs : float (optional)
The sample rate (default: from the binner).
window : string or array
Specify the window function.
Options: 1, None, 'hann', 'hamm'
noise : float or array-like
A vector of the noise levels of the velocity data with
the same first dimension as the velocity vector.
noise : numeric
The noise level in the same units as velocity
n_bin : int (optional)
The bin-size (default: from the binner).
n_fft : int (optional)
Expand Down Expand Up @@ -1107,14 +1107,16 @@ def calc_psd(self, veldat,
).astype('float32')

# Spectra, if input is full velocity or a single array
if len(vel.shape) == 2:
assert vel.shape[0] == 3, "Function can only handle 1D or 3D arrays." \
" If ADCP data, please select a specific depth bin."
if (noise is not None) and (np.shape(noise)[0] != 3):
raise Exception(
'Noise should have same first dimension as velocity')
if len(vel.shape) >= 2:
if vel.shape[0] == 3:
raise ValueError("Function can only handle 1D or 3D arrays."
" If ADCP data, please select a specific depth bin.")
if noise is not None:
if np.size(noise) != 3:
raise ValueError('Noise is expected to be an array of 3 scalars')
else:
noise = np.array([0, 0, 0])
noise = 0

out = np.empty(self._outshape_fft(vel[:3].shape, n_fft=n_fft, n_bin=n_bin),
dtype=np.float32)
for idx in range(3):
Expand All @@ -1131,11 +1133,11 @@ def calc_psd(self, veldat,
'freq': freq}
dims = ['S', 'time', 'freq']
else:
if (noise is not None) and (len(np.shape(noise)) > 1):
raise Exception(
'Noise should have same first dimension as velocity')
if noise is not None:
if np.size(noise) > 1:
raise ValueError('Noise is expected to be a scalar')
else:
noise = np.array(0)
noise = 0
out = self.calc_psd_base(vel,
fs=fs,
noise=noise,
Expand Down

0 comments on commit d0ddc77

Please sign in to comment.