Calibration hybrid milidecades

[JPalmerK]

Hey there,
Not sure if this issue fits better in pyhydrophone or in pypam but I've been testing the hybrid milidecade band analysis a csv file. I've ensured that the data are being properly read into the hydrophone object (get_freq_cal within the hydrophone object, line 154) but the hybridmilidecades do not seem to change with or without the frequency dependent sensitivity included. Any help would be most appreciated. Thanks!

# Create a custom hydrophone class with the model number, serial number, and
# name (Dummy)
model = "alasdfasdfae"
name = "asdfsddd"
serial_number = 8675309
drifter_noCal = pyhy.RTSys(
    name=name, 
    model=model, 
    sensitivity=-172.8, 
    serial_number=serial_number,
    preamp_gain= 10,
    Vpp =2,
    mode = 'broadband')

drifter_Cal = pyhy.RTSys(
    name=name, 
    model=model, 
    sensitivity=-172.8, 
    serial_number=serial_number,
    preamp_gain= 10,
    Vpp =2,
    mode = 'broadband',
    calibration_file = calFileLoc)


# First, decide band to study. The top frequency should not be higher than the nyquist frequency (sampling rate/2)
band = [0, 4000]

# Then, set the nfft to double the sampling rate. If you want to pass None to your band, that is also an option, but
# then you need to know the sampling frequency to choose the nfft.
nfft = band[1] * 2  # or nfft = 8000

# Set the band to 1 minute
binsize = 60.0

include_dirs = False
zipped_files = False
dc_subtract = True

# Create an asa for the frequency specific calibration and non-frequency specific calibration 
asa_cal = pypam.ASA(
    hydrophone=drifter_Cal,
    folder_path="C:\\Users\\kaity\\Documents\\GitHub\\SPACIOUS-Propagation-Modes\\ExampleData",
    binsize=binsize,
    nfft=nfft,
    timezone="UTC",
    include_dirs=include_dirs,
    zipped=zipped_files,
    dc_subtract=dc_subtract)

asa_nocal = pypam.ASA(
    hydrophone=drifter_noCal,
    folder_path="C:\\Users\\kaity\\Documents\\GitHub\\SPACIOUS-Propagation-Modes\\ExampleData",
    binsize=binsize,
    nfft=nfft,
    timezone="UTC",
    include_dirs=include_dirs,
    zipped=zipped_files,
    dc_subtract=dc_subtract)


# Compute the hybrid millidecade bands

milli_psd_cal = asa_cal.hybrid_millidecade_bands(
    db=True, method="density", band=band, percentiles=None)
aa_cal = 20*np.log10(milli_psd_cal.band_density.values)


milli_psd_nocal = asa_nocal.hybrid_millidecade_bands(
    db=True, method="density", band=band, percentiles=None)
aa_nocal = 20*np.log10(milli_psd_nocal.band_density.values)

# Ach! they are the same.
difval = aa_cal-aa_nocal

[ryjombari]

Hi there. It looks like both your noCal and Cal cases specify a flat sensitivity (a single number for sensitivity), which would explain no difference in results. To use a frequency dependent calibration, you would specify the file containing that calibration data, and you would not specify a single number for sensitivity. Does your hydrophone have frequency dependent sensitivity or flat?

[cparcerisas]

Hi Kait!
Currently pypam does not automatically implement the frequency-dependent calibration directly (but it is in the TODOs).
How it is implemented atm is that only the flat sensitivity is used automatically in any method used from pypam.
To be able to use freq-dependent calibration, at the moment the way to go is this one (continuing from the end of your code).

(you can see the original example here: https://pyhydrophone.readthedocs.io/en/latest/_auto_examples/plot_freq_calibration.html#sphx-glr-auto-examples-plot-freq-calibration-py)

frequency_increment = drifter_Cal.freq_cal_inc(milli_psd_cal_corrected.frequency_bins)
milli_psd_cal_corrected = milli_psd_cal_corrected['millidecade_bands'] + frequency_increment['inc_value']

Let me know if this works for you!

Clea

[JPalmerK]

Thanks so much @cparcerisas! I'll give it a go now an let you you know how it works!

Update: I had do to some tweaking to the example because xarrays can't broadcast. This code also includes generic hydrophone structure ('custom') and I've made a local mortification to get the calibration from an nc file. Code and results are below.

import matplotlib as mpl
import matplotlib.pyplot as plt

# 1. Reset everything to the default Matplotlib style:
mpl.rcdefaults()     # → clears out any custom rcParams/style sheets

# 2. Make sure usetex is off (in case something later turns it on again):
mpl.rcParams['text.usetex'] = False

# 3. (Optionally) choose the default style explicitly:
plt.style.use('default')


model = "RandomPhone"
name = "Bob"
serial_number = 67416073
calFileLoc = '\\ExampleData\\sg679_CalCurCEAS_Aug2024_sensitivity_2025-05-19.nc'


drifter_Cal = pyhy.custom(
    name=name, 
    model=model, 
    sensitivity=0, 
    serial_number=serial_number,
    preamp_gain= 0,
    Vpp =2,
    #val="sensitivity",
    calibration_file = calFileLoc)

![Image](https://github.com/user-attachments/assets/84057727-5f93-45e7-bd87-2a5365ebc1dd)

# Get the end-to-end calibration values
# drifter_noCal.end_to_end_calibration()
drifter_Cal.end_to_end_calibration()

# Set up the acoustic study
# First, decide band to study. The top frequency should not be higher than the nyquist frequency (sampling rate/2)
band = [0, 90000]

# Then, set the nfft to double the sampling rate. 
nfft = band[1] * 2  # or nfft = 8000

# Set the band to 1 minute
binsize = 60.0

# Select features of 
include_dirs = False
zipped_files = False
dc_subtract = True
asa_cal = pypam.ASA(
    hydrophone=drifter_Cal,
    folder_path="\\ExampleData\\GliderData\\",
    binsize=binsize,
    nfft=nfft,
    timezone="UTC",
    include_dirs=include_dirs,
    zipped=zipped_files,
    dc_subtract=dc_subtract,
    calibration=-1)


# Compute the hybrid millidecade bands
# You can choose 'density' or 'spectrum' as a method
milli_psd_cal = asa_cal.hybrid_millidecade_bands(
    db=True, method="density", band=band, percentiles=None)

# Pull the frequency values from the hybrid milidecade bands
frequencies = milli_psd_cal.frequency.values
frequency_increment = drifter_Cal.freq_cal_inc(frequencies)

frequencies = milli_psd_cal.frequency_bins.values
frequency_increment = drifter_Cal.freq_cal_inc(frequencies)


# Plot the incement -giggle test
plt.plot(frequency_increment['frequency'][1:], -frequency_increment['inc_value'][1:])
plt.xlabel('Frequency (Hz)')
plt.ylabel('dB')
plt.title('End-to-End Sensitivity (hyd, preamp, voltage)')

# 1) grab the PSD DataArray (shape = (9, 2752))
psd = milli_psd_cal["millidecade_bands"]

# 2) subtract the 1-D increment (length 2752). xarray will broadcast it along the first axis:
calibrated = frequency_increment["inc_value"].values+psd  

# 3) put it back into a copy of the original Dataset:
milli_psd_cal_dB_upa = milli_psd_cal.copy()
milli_psd_cal_dB_upa["millidecade_bands"] = calibrated


# Plot the spectrum mean with the standard deviation
pypam.plots.plot_spectrum_median(
    milli_psd_cal, data_var="millidecade_bands", frequency_coord="frequency_bins")

# Plot the spectrum mean with the standard deviation
pypam.plots.plot_spectrum_median(
    milli_psd_cal_dB_upa, data_var="millidecade_bands", frequency_coord="frequency_bins")