Download this notebook by clicking on the download icon or find it in the repository

Fin Call Index

Copyright (c) 2021, MBARI

• Distributed under the terms of the GPL License
• Maintainer: ryjo@mbari.org
• Author: John Ryan ryjo@mbari.org

Fin whale song¶

---

Fin whales produce rhythmic structured sequences of sound, typically in a series of pulses described as 'song'. This tutorial describes use of the Pacific Ocean Sound Recordings archive to examine temporal patterns of occurrence of fin whale song.

If you use this data set, please cite our project.

Data Overview¶

---

Recording site¶

The recording site is located on the continental slope of the eastern North Pacific, within Monterey Bay National Marine Sanctuary. The region is known to be important foraging habitat for the regional whale populations.

Hydrophone calibration¶

For the low-frequency (2 kHz) data, calibration data are not frequency dependent; a single low-frequency calibration value is used. Its value depends on time of data collection, as two hydrophones have been deployed sequentially at the same site. Before 14 June 2017, the calibration value is -168.8 dB re V / uPa (measured at 26 Hz). After this date the value is -177.9 dB re V / uPa (measured at 250 Hz). See also:

• https://bitbucket.org/mbari/pacific-sound/src/master/MBARI_MARS_Hydrophone_Deployment01.json
• https://bitbucket.org/mbari/pacific-sound/src/master/MBARI_MARS_Hydrophone_Deployment02.json

The first hydrophone exhibited calibration drift, while the second (deployed 13 June 2017 and currently operational) has not. This observation is consistent with differences in the technologies of the two instruments. However, for this application the calibration drift of the first hydrophone is not problematic because the CI is computed as a signal to noise ratio. Therefore, time-series analysis of CI can reliably span the full archive.

Data files and archive organization¶

The decimated audio data are in daily WAV files in an s3 bucket named pacific-sound-2khz, grouped by year and month. Buckets are stored as objects, so the data are not physically stored in folders or directories as you may be famaliar with, but you can think of it conceptually as follows:

pacific-sound-2khz
      |
      ----2020
        |
        |----01
        ...
        |----12

Install required dependencies¶

First, let's install the required software dependencies.

If you are using this notebook in a cloud environment, select a Python3 compatible kernel and run this next section. This only needs to be done once for the duration of this notebook.

If you are working on local computer, you can skip this next cell. Change your kernel to pacific-sound-notebooks, which you installed according to the instructions in the README - this has all the dependencies that are needed.

In [1]:

!pip install -q boto3 --quiet
!pip install -q soundfile --quiet
!pip install -q scipy --quiet
!pip install -q numpy --quiet
!pip install -q matplotlib --quiet

!pip install -q boto3 --quiet !pip install -q soundfile --quiet !pip install -q scipy --quiet !pip install -q numpy --quiet !pip install -q matplotlib --quiet

     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 139.3/139.3 kB 1.6 MB/s eta 0:00:00
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 11.9/11.9 MB 26.8 MB/s eta 0:00:00
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 82.1/82.1 kB 8.7 MB/s eta 0:00:00

Import all packages¶

In [2]:

import boto3, botocore
from botocore import UNSIGNED
from botocore.client import Config
from six.moves.urllib.request import urlopen
import io
import scipy
from scipy import signal
import numpy as np
import soundfile as sf
import matplotlib.pyplot as plt

import boto3, botocore from botocore import UNSIGNED from botocore.client import Config from six.moves.urllib.request import urlopen import io import scipy from scipy import signal import numpy as np import soundfile as sf import matplotlib.pyplot as plt

Data Access¶

---

This section covers file listing, metadata retrieval, and data loading.

List files¶

Files are organized by year and month; list all of the files available for one month of one year.

In [3]:

s3 = boto3.client('s3',
    aws_access_key_id='',
    aws_secret_access_key='',
    config=Config(signature_version=UNSIGNED))

s3 = boto3.client('s3', aws_access_key_id='', aws_secret_access_key='', config=Config(signature_version=UNSIGNED))

In [4]:

year = 2023
month = 12
bucket = 'pacific-sound-2khz'

for obj in s3.list_objects_v2(Bucket=bucket, Prefix=f'{year:04d}/{month:02d}')['Contents']:
    print(obj['Key'])

year = 2023 month = 12 bucket = 'pacific-sound-2khz' for obj in s3.list_objects_v2(Bucket=bucket, Prefix=f'{year:04d}/{month:02d}')['Contents']: print(obj['Key'])

2023/12/MARS-20231201T000000Z-2kHz.wav
2023/12/MARS-20231202T000000Z-2kHz.wav
2023/12/MARS-20231203T000000Z-2kHz.wav
2023/12/MARS-20231204T000000Z-2kHz.wav
2023/12/MARS-20231205T000000Z-2kHz.wav
2023/12/MARS-20231206T000000Z-2kHz.wav
2023/12/MARS-20231207T000000Z-2kHz.wav
2023/12/MARS-20231208T000000Z-2kHz.wav
2023/12/MARS-20231209T000000Z-2kHz.wav
2023/12/MARS-20231210T000000Z-2kHz.wav
2023/12/MARS-20231211T000000Z-2kHz.wav
2023/12/MARS-20231212T000000Z-2kHz.wav
2023/12/MARS-20231213T000000Z-2kHz.wav
2023/12/MARS-20231214T000000Z-2kHz.wav
2023/12/MARS-20231215T000000Z-2kHz.wav
2023/12/MARS-20231216T000000Z-2kHz.wav
2023/12/MARS-20231217T000000Z-2kHz.wav
2023/12/MARS-20231218T000000Z-2kHz.wav
2023/12/MARS-20231219T000000Z-2kHz.wav
2023/12/MARS-20231220T000000Z-2kHz.wav
2023/12/MARS-20231221T000000Z-2kHz.wav
2023/12/MARS-20231222T000000Z-2kHz.wav
2023/12/MARS-20231223T000000Z-2kHz.wav
2023/12/MARS-20231224T000000Z-2kHz.wav
2023/12/MARS-20231225T000000Z-2kHz.wav
2023/12/MARS-20231226T000000Z-2kHz.wav
2023/12/MARS-20231227T000000Z-2kHz.wav
2023/12/MARS-20231228T000000Z-2kHz.wav
2023/12/MARS-20231229T000000Z-2kHz.wav
2023/12/MARS-20231230T000000Z-2kHz.wav
2023/12/MARS-20231231T000000Z-2kHz.wav

Retrieve metadata¶

Read and show metadata for a single daily file.

In [5]:

filename = 'MARS-20231201T000000Z-2kHz.wav'
key = f'{year:04d}/{month:02d}/{filename}'

url = f'https://{bucket}.s3.amazonaws.com/{key}'

sf.info(io.BytesIO(urlopen(url).read(1_000)), verbose=True)

filename = 'MARS-20231201T000000Z-2kHz.wav' key = f'{year:04d}/{month:02d}/{filename}' url = f'https://{bucket}.s3.amazonaws.com/{key}' sf.info(io.BytesIO(urlopen(url).read(1_000)), verbose=True)

Out[5]:

<_io.BytesIO object at 0x7a4c6ebcad40>
samplerate: 2000 Hz
channels: 1
duration: 222 samples
format: WAV (Microsoft) [WAV]
subtype: Signed 24 bit PCM [PCM_24]
endian: FILE
sections: 1
frames: 222
extra_info: """
    Length : 1000
    RIFF : 518400324 (should be 992)
    WAVE
    fmt  : 16
      Format        : 0x1 => WAVE_FORMAT_PCM
      Channels      : 1
      Sample Rate   : 2000
      Block Align   : 3
      Bit Width     : 24
      Bytes/sec     : 6000
    LIST : 280
      INFO
        INAM : MBARI ocean audio data, start 20231201T000000 UTC
        ICMT : If you use these data, please cite https://doi.org/10.1109/OCEANS.2016.7761363. Recording metadata can be found at https://bitbucket.org/mbari/pacific-sound/src/master/MBARI_MARS_Hydrophone_Deployment02.json.
    data : 518400000 (should be 668)
    End
    """

Load data¶

Read a single daily file.

In [6]:

# read full-day of data
print(f'Reading from {url}')
v, sample_rate = sf.read(io.BytesIO(urlopen(url).read()),dtype='float32')
v = v*3   # convert scaled voltage to volts
nsec = (v.size)/sample_rate # number of seconds in vector
print(f'Read {nsec} seconds of data')

# read full-day of data print(f'Reading from {url}') v, sample_rate = sf.read(io.BytesIO(urlopen(url).read()),dtype='float32') v = v*3 # convert scaled voltage to volts nsec = (v.size)/sample_rate # number of seconds in vector print(f'Read {nsec} seconds of data')

Reading from https://pacific-sound-2khz.s3.amazonaws.com/2023/12/MARS-20231201T000000Z-2kHz.wav
Read 86400.0 seconds of data

A view of fin whale song¶

---

To understand the method of quantifying song occurrence using an energy metric, it is useful to first consider the attributes of fin whale song. Analysis approaches include (1) detecting, classifying, and counting calls, and (2) quantifying the energy within the frequency band of the call, relative to that at background frequencies. The first approach becomes difficult during periods when the whales chorus because the presence of overlapping calls thwarts distinction of individual calls. The second approach can be applied consistently regardless of whether or not whale calls overlap, and it is effective for quantifying the integrated signal received from many calling whales.

In [7]:

# Compute spectrogram
w = scipy.signal.get_window('hann',sample_rate)
f, t, psd = scipy.signal.spectrogram(v, sample_rate,nperseg=sample_rate,noverlap=0,window=w,nfft=sample_rate)
sens = -168.8  # hydrophone sensitivity at 26 Hz
psd = 10*np.log10(psd) - sens
print(f':: psd.shape           = {psd.shape}')
print(f':: f.size              = {f.size}')
print(f':: t.size              = {t.size}')

# Compute spectrogram w = scipy.signal.get_window('hann',sample_rate) f, t, psd = scipy.signal.spectrogram(v, sample_rate,nperseg=sample_rate,noverlap=0,window=w,nfft=sample_rate) sens = -168.8 # hydrophone sensitivity at 26 Hz psd = 10*np.log10(psd) - sens print(f':: psd.shape = {psd.shape}') print(f':: f.size = {f.size}') print(f':: t.size = {t.size}')

:: psd.shape           = (1001, 86400)
:: f.size              = 1001
:: t.size              = 86400

<ipython-input-7-9921a84d1056>:5: RuntimeWarning: divide by zero encountered in log10
  psd = 10*np.log10(psd) - sens

In [8]:

# View 15 minutes
start_hour = 15
start_sec = int(start_hour * 3600 + 1)
end_sec = start_sec+900-1
psd_subset = psd[:,start_sec:end_sec]
plt.figure(dpi=200, figsize = [9,3])
plt.imshow(psd_subset,aspect='auto',origin='lower',vmin=45,vmax=95)
plt.colorbar()
plt.ylim(10,50)
plt.xlabel('Second of hour 15')
plt.ylabel('Frequency (Hz)')
plt.title('Spectrum level (dB re 1 $\mu$Pa$^2$/Hz)')
plt.annotate("fin whale pulses",(400,30),color='w')

# View 15 minutes start_hour = 15 start_sec = int(start_hour * 3600 + 1) end_sec = start_sec+900-1 psd_subset = psd[:,start_sec:end_sec] plt.figure(dpi=200, figsize = [9,3]) plt.imshow(psd_subset,aspect='auto',origin='lower',vmin=45,vmax=95) plt.colorbar() plt.ylim(10,50) plt.xlabel('Second of hour 15') plt.ylabel('Frequency (Hz)') plt.title('Spectrum level (dB re 1 $\mu$Pa$^2$/Hz)') plt.annotate("fin whale pulses",(400,30),color='w')

Out[8]:

Text(400, 30, 'fin whale pulses')

Fin whale pulse calls occurred throughout the time period shown. The stronger signals with a regular cadence are likely from a single whale, which appears to have paused for a surface breathing interval between ~600 and 700 seconds. These calls are part of a greater chorus evident in the background, from multiple fin whales calling from different locations (distances) relative to the hydrophone.

Call index¶

---

To consider the peak and background frequencies used in the calculation of the call index, let's examine the average spectrum levels for period shown above.

In [9]:

m = np.mean(psd_subset,axis=1)
plt.figure(dpi=200, figsize = [7,3])
plt.plot(m,f,'k')
plt.ylim(10,38)
plt.xlim(62,80)
plt.plot([62, 80],[20, 20],'b--')
plt.plot([62, 80],[21, 21],'b--')
plt.plot([62, 80],[12, 12],'r--')
plt.plot([62, 80],[34, 34],'r--')
plt.xlabel('Mean spectrum level (dB re 1 $\mu$Pa$^2$/Hz)')
plt.ylabel('Frequency (Hz)')

m = np.mean(psd_subset,axis=1) plt.figure(dpi=200, figsize = [7,3]) plt.plot(m,f,'k') plt.ylim(10,38) plt.xlim(62,80) plt.plot([62, 80],[20, 20],'b--') plt.plot([62, 80],[21, 21],'b--') plt.plot([62, 80],[12, 12],'r--') plt.plot([62, 80],[34, 34],'r--') plt.xlabel('Mean spectrum level (dB re 1 $\mu$Pa$^2$/Hz)') plt.ylabel('Frequency (Hz)')

Out[9]:

Text(0, 0.5, 'Frequency (Hz)')

Energy within the peak caused by fin whale pulses (frequencies indicated by dashed blue lines) is much higher than the background (frequencies indicated by dashed red lines). The call index (CI) is simply the ratio of signal (average spectrum level at the peak frequencies) to noise (average spectrum level at the background frequencies). A value of 1 effectively means no signal (no difference between peak and background). To examine response of the index to both background whale calls and strong individual calls, let's compute CI at 1-second resolution.

In [10]:

# find the frequencies of the peak and average spectrum levels
p1 = psd_subset[f==20]; p2 = psd_subset[f==21];
pk = np.squeeze(np.array([p1,p2]))
pk = np.mean(pk,axis=0); pk.shape
# plt.plot(pk)
# find the frequencies of the background and average
b1 = psd_subset[f==12]; b2 = psd_subset[f==34];
bg = np.squeeze(np.array([b1,b2]))
bg = np.mean(bg,axis=0); bg.shape
# plt.plot(pk/bg)
# CI
CI = pk/bg;
plt.figure(dpi=200, figsize = [9,3])
plt.plot([0, 900],[1, 1],'m--')
plt.plot(CI)
plt.xlabel('Seconds')
plt.ylabel('CI')

# find the frequencies of the peak and average spectrum levels p1 = psd_subset[f==20]; p2 = psd_subset[f==21]; pk = np.squeeze(np.array([p1,p2])) pk = np.mean(pk,axis=0); pk.shape # plt.plot(pk) # find the frequencies of the background and average b1 = psd_subset[f==12]; b2 = psd_subset[f==34]; bg = np.squeeze(np.array([b1,b2])) bg = np.mean(bg,axis=0); bg.shape # plt.plot(pk/bg) # CI CI = pk/bg; plt.figure(dpi=200, figsize = [9,3]) plt.plot([0, 900],[1, 1],'m--') plt.plot(CI) plt.xlabel('Seconds') plt.ylabel('CI')

Out[10]:

Text(0, 0.5, 'CI')

This high-resolution time series emphasizes the signal of the loudest singer (closest to the hydrophone) with a consistent interpulse interval for many of the highest CI events. It also shows oscillations in the background total sound energy received from calling fin whales.

Time-series analysis¶

---

Let's produce daily statistics for the same month as the brief period examined above. For time-series analysis, additional considerations for the methods arise. The first is that the recording time-series contains gaps ranging from minutes to days. Each daily 2 kHz wav file is fully populated (contains 1 day of data), and missing data periods are indicated by zero values. Therefore, gaps can be identified as contiguous segments of data having zero values. The second consideration is that computing power spectral density over periods longer than 1 second (as above) can be useful data reduction in time-series analysis. Therefore, below we will use a 60-second window for psd calculation, and we will screen missing data periods.

In [11]:

# Determine the number of data files available for the month
num_files = 0
for obj in s3.list_objects_v2(Bucket=bucket, Prefix=f'{year}/{month}')['Contents']:
    print(obj['Key'])
    num_files = num_files+1
print(f':: number of files to process = {num_files}')

# Determine the number of data files available for the month num_files = 0 for obj in s3.list_objects_v2(Bucket=bucket, Prefix=f'{year}/{month}')['Contents']: print(obj['Key']) num_files = num_files+1 print(f':: number of files to process = {num_files}')

2023/12/MARS-20231201T000000Z-2kHz.wav
2023/12/MARS-20231202T000000Z-2kHz.wav
2023/12/MARS-20231203T000000Z-2kHz.wav
2023/12/MARS-20231204T000000Z-2kHz.wav
2023/12/MARS-20231205T000000Z-2kHz.wav
2023/12/MARS-20231206T000000Z-2kHz.wav
2023/12/MARS-20231207T000000Z-2kHz.wav
2023/12/MARS-20231208T000000Z-2kHz.wav
2023/12/MARS-20231209T000000Z-2kHz.wav
2023/12/MARS-20231210T000000Z-2kHz.wav
2023/12/MARS-20231211T000000Z-2kHz.wav
2023/12/MARS-20231212T000000Z-2kHz.wav
2023/12/MARS-20231213T000000Z-2kHz.wav
2023/12/MARS-20231214T000000Z-2kHz.wav
2023/12/MARS-20231215T000000Z-2kHz.wav
2023/12/MARS-20231216T000000Z-2kHz.wav
2023/12/MARS-20231217T000000Z-2kHz.wav
2023/12/MARS-20231218T000000Z-2kHz.wav
2023/12/MARS-20231219T000000Z-2kHz.wav
2023/12/MARS-20231220T000000Z-2kHz.wav
2023/12/MARS-20231221T000000Z-2kHz.wav
2023/12/MARS-20231222T000000Z-2kHz.wav
2023/12/MARS-20231223T000000Z-2kHz.wav
2023/12/MARS-20231224T000000Z-2kHz.wav
2023/12/MARS-20231225T000000Z-2kHz.wav
2023/12/MARS-20231226T000000Z-2kHz.wav
2023/12/MARS-20231227T000000Z-2kHz.wav
2023/12/MARS-20231228T000000Z-2kHz.wav
2023/12/MARS-20231229T000000Z-2kHz.wav
2023/12/MARS-20231230T000000Z-2kHz.wav
2023/12/MARS-20231231T000000Z-2kHz.wav
:: number of files to process = 31

Single-day processing¶

First, define the standard daily processing steps.

In [12]:

# define standard segment processing
w = scipy.signal.get_window('hann',sample_rate)
total_sec = v.size / sample_rate # total number of seconds in vector
spa = 60  # seconds per average
num_segments = int(total_sec/spa)
print(f'{num_segments} segments of length {spa} seconds in {total_sec} seconds of audio')

# initialize empty spectrogram
nyquist_freq = int(sample_rate/2+1)
sg = np.empty((nyquist_freq, num_segments), float)
print(f':: individual spectrogram shape   =   {sg.shape}')

# process spectrogram, capturing number of zero values in data
for x in range(0,num_segments):
  cstart = x*spa*sample_rate; cend = (x+1)*spa*sample_rate
  cv = v[cstart:cend]
  f,psd = scipy.signal.welch(cv,fs=sample_rate,window=w,nfft=sample_rate)
  psd = 10*np.log10(psd) - sens
  sg[:,x] = psd

# plot spectrogram for example day
plt.figure(dpi=200, figsize = [9,3])
plt.imshow(sg,aspect='auto',origin='lower',vmin=40,vmax=100)
plt.yscale('log')
plt.ylim(10,1000)
plt.colorbar()
plt.xlabel('Minute of day')
plt.ylabel('Frequency (Hz)')
plt.title('Calibrated spectrum levels')

# define standard segment processing w = scipy.signal.get_window('hann',sample_rate) total_sec = v.size / sample_rate # total number of seconds in vector spa = 60 # seconds per average num_segments = int(total_sec/spa) print(f'{num_segments} segments of length {spa} seconds in {total_sec} seconds of audio') # initialize empty spectrogram nyquist_freq = int(sample_rate/2+1) sg = np.empty((nyquist_freq, num_segments), float) print(f':: individual spectrogram shape = {sg.shape}') # process spectrogram, capturing number of zero values in data for x in range(0,num_segments): cstart = x*spa*sample_rate; cend = (x+1)*spa*sample_rate cv = v[cstart:cend] f,psd = scipy.signal.welch(cv,fs=sample_rate,window=w,nfft=sample_rate) psd = 10*np.log10(psd) - sens sg[:,x] = psd # plot spectrogram for example day plt.figure(dpi=200, figsize = [9,3]) plt.imshow(sg,aspect='auto',origin='lower',vmin=40,vmax=100) plt.yscale('log') plt.ylim(10,1000) plt.colorbar() plt.xlabel('Minute of day') plt.ylabel('Frequency (Hz)') plt.title('Calibrated spectrum levels')

1440 segments of length 60 seconds in 86400.0 seconds of audio
:: individual spectrogram shape   =   (1001, 1440)

Out[12]:

Text(0.5, 1.0, 'Calibrated spectrum levels')

In the frequency band of peak energy from fin whale calls, 20-21 Hz, we can see calls present throughout the day, with greatest intensity centered near minute 900.

To consider day-to-day variability, let's process a month of data, representing mean CI for each day. We will still process in 1-minute segments, which allows us to identify periods of missing data and exclude them from analysis.

In [13]:

# Batch process the month of daily files
F = -1   # first row index will be 0 after increment
spa = 60  # seconds per average
num_segments = int(86400/spa)   # number of processing segments
# Initialize arrays to hold CI and zero count results
# dimensions: number of days x number of segments per day
CIm = np.zeros((num_files,num_segments))
zero_count = np.zeros((num_files,num_segments))

sample_rate = 2000   # data sample rate
w = scipy.signal.get_window('hann',sample_rate)   # 1-second window for 2 kHz data
nyquist_freq = int(sample_rate/2+1)   # number of output frequencies in spectrogram

for obj in s3.list_objects_v2(Bucket=bucket, Prefix=f'{year}/{month}')['Contents']:
  F = F+1
  # read file
  filename = obj['Key']
  url = f'https://{bucket}.s3.amazonaws.com/{filename}'
  print(f'Reading from {url}')
  v, fs = sf.read(io.BytesIO(urlopen(url).read()),dtype='float32')
  v = v*3   # convert scaled voltage to volts

  # initialize empty spectrogram
  sg = np.empty((nyquist_freq, num_segments), float)

  # process spectrogram, capturing number of zero values in data
  for x in range(0,num_segments):
    cstart = x*spa*fs; cend = (x+1)*spa*fs; cv = v[cstart:cend]
    f,psd = scipy.signal.welch(cv,fs,window=w,nfft=fs)
    psd = 10*np.log10(psd) - sens
    sg[:,x] = psd
    zero_count[F,x] = (cv == 0).sum(0)

  # CI
  # find the frequencies of the peak and average spectrum levels
  p1 = sg[f==20]; p2 = sg[f==21];
  pk = np.squeeze(np.array([p1,p2]))
  pk = np.mean(pk,axis=0);
  # find the frequencies of the background and average
  b1 = sg[f==12]; b2 = sg[f==34];
  bg = np.squeeze(np.array([b1,b2]))
  bg = np.mean(bg,axis=0);

  # CI
  CIm[F,:] = pk/bg;

# Batch process the month of daily files F = -1 # first row index will be 0 after increment spa = 60 # seconds per average num_segments = int(86400/spa) # number of processing segments # Initialize arrays to hold CI and zero count results # dimensions: number of days x number of segments per day CIm = np.zeros((num_files,num_segments)) zero_count = np.zeros((num_files,num_segments)) sample_rate = 2000 # data sample rate w = scipy.signal.get_window('hann',sample_rate) # 1-second window for 2 kHz data nyquist_freq = int(sample_rate/2+1) # number of output frequencies in spectrogram for obj in s3.list_objects_v2(Bucket=bucket, Prefix=f'{year}/{month}')['Contents']: F = F+1 # read file filename = obj['Key'] url = f'https://{bucket}.s3.amazonaws.com/{filename}' print(f'Reading from {url}') v, fs = sf.read(io.BytesIO(urlopen(url).read()),dtype='float32') v = v*3 # convert scaled voltage to volts # initialize empty spectrogram sg = np.empty((nyquist_freq, num_segments), float) # process spectrogram, capturing number of zero values in data for x in range(0,num_segments): cstart = x*spa*fs; cend = (x+1)*spa*fs; cv = v[cstart:cend] f,psd = scipy.signal.welch(cv,fs,window=w,nfft=fs) psd = 10*np.log10(psd) - sens sg[:,x] = psd zero_count[F,x] = (cv == 0).sum(0) # CI # find the frequencies of the peak and average spectrum levels p1 = sg[f==20]; p2 = sg[f==21]; pk = np.squeeze(np.array([p1,p2])) pk = np.mean(pk,axis=0); # find the frequencies of the background and average b1 = sg[f==12]; b2 = sg[f==34]; bg = np.squeeze(np.array([b1,b2])) bg = np.mean(bg,axis=0); # CI CIm[F,:] = pk/bg;

Reading from https://pacific-sound-2khz.s3.amazonaws.com/2023/12/MARS-20231201T000000Z-2kHz.wav
Reading from https://pacific-sound-2khz.s3.amazonaws.com/2023/12/MARS-20231202T000000Z-2kHz.wav
Reading from https://pacific-sound-2khz.s3.amazonaws.com/2023/12/MARS-20231203T000000Z-2kHz.wav
Reading from https://pacific-sound-2khz.s3.amazonaws.com/2023/12/MARS-20231204T000000Z-2kHz.wav
Reading from https://pacific-sound-2khz.s3.amazonaws.com/2023/12/MARS-20231205T000000Z-2kHz.wav
Reading from https://pacific-sound-2khz.s3.amazonaws.com/2023/12/MARS-20231206T000000Z-2kHz.wav
Reading from https://pacific-sound-2khz.s3.amazonaws.com/2023/12/MARS-20231207T000000Z-2kHz.wav
Reading from https://pacific-sound-2khz.s3.amazonaws.com/2023/12/MARS-20231208T000000Z-2kHz.wav
Reading from https://pacific-sound-2khz.s3.amazonaws.com/2023/12/MARS-20231209T000000Z-2kHz.wav
Reading from https://pacific-sound-2khz.s3.amazonaws.com/2023/12/MARS-20231210T000000Z-2kHz.wav
Reading from https://pacific-sound-2khz.s3.amazonaws.com/2023/12/MARS-20231211T000000Z-2kHz.wav
Reading from https://pacific-sound-2khz.s3.amazonaws.com/2023/12/MARS-20231212T000000Z-2kHz.wav
Reading from https://pacific-sound-2khz.s3.amazonaws.com/2023/12/MARS-20231213T000000Z-2kHz.wav
Reading from https://pacific-sound-2khz.s3.amazonaws.com/2023/12/MARS-20231214T000000Z-2kHz.wav
Reading from https://pacific-sound-2khz.s3.amazonaws.com/2023/12/MARS-20231215T000000Z-2kHz.wav
Reading from https://pacific-sound-2khz.s3.amazonaws.com/2023/12/MARS-20231216T000000Z-2kHz.wav
Reading from https://pacific-sound-2khz.s3.amazonaws.com/2023/12/MARS-20231217T000000Z-2kHz.wav
Reading from https://pacific-sound-2khz.s3.amazonaws.com/2023/12/MARS-20231218T000000Z-2kHz.wav
Reading from https://pacific-sound-2khz.s3.amazonaws.com/2023/12/MARS-20231219T000000Z-2kHz.wav
Reading from https://pacific-sound-2khz.s3.amazonaws.com/2023/12/MARS-20231220T000000Z-2kHz.wav
Reading from https://pacific-sound-2khz.s3.amazonaws.com/2023/12/MARS-20231221T000000Z-2kHz.wav
Reading from https://pacific-sound-2khz.s3.amazonaws.com/2023/12/MARS-20231222T000000Z-2kHz.wav
Reading from https://pacific-sound-2khz.s3.amazonaws.com/2023/12/MARS-20231223T000000Z-2kHz.wav
Reading from https://pacific-sound-2khz.s3.amazonaws.com/2023/12/MARS-20231224T000000Z-2kHz.wav
Reading from https://pacific-sound-2khz.s3.amazonaws.com/2023/12/MARS-20231225T000000Z-2kHz.wav
Reading from https://pacific-sound-2khz.s3.amazonaws.com/2023/12/MARS-20231226T000000Z-2kHz.wav
Reading from https://pacific-sound-2khz.s3.amazonaws.com/2023/12/MARS-20231227T000000Z-2kHz.wav
Reading from https://pacific-sound-2khz.s3.amazonaws.com/2023/12/MARS-20231228T000000Z-2kHz.wav
Reading from https://pacific-sound-2khz.s3.amazonaws.com/2023/12/MARS-20231229T000000Z-2kHz.wav
Reading from https://pacific-sound-2khz.s3.amazonaws.com/2023/12/MARS-20231230T000000Z-2kHz.wav
Reading from https://pacific-sound-2khz.s3.amazonaws.com/2023/12/MARS-20231231T000000Z-2kHz.wav

Viewing the arrays containing 1-minute CI and missing data counts.

In [14]:

# View CI
plt.figure(dpi=200, figsize = [9,3])
plt.imshow(CIm,aspect='auto',origin='lower')
plt.xlabel('Minute of day')
plt.ylabel('Day of month')
plt.title('Fin whale call index')
plt.colorbar()

# View CI plt.figure(dpi=200, figsize = [9,3]) plt.imshow(CIm,aspect='auto',origin='lower') plt.xlabel('Minute of day') plt.ylabel('Day of month') plt.title('Fin whale call index') plt.colorbar()

Out[14]:

<matplotlib.colorbar.Colorbar at 0x7a4c37b15c00>

In the above CI analysis result we see significant variations in fin whale call activity at a range of time scales. For example, calling was stronger during the first week, particularly during the latter half of the day (UTC). Poor sampling of any one-minute temporal bin yields untrustworthy results, so we also want to examine sampling coverage.

In [15]:

# View percent of each 1 minute sample having missing data
percent_missing = 100*zero_count/(spa*sample_rate)
plt.figure(dpi=200, figsize = [9,3])
plt.imshow(percent_missing,aspect='auto',origin='lower')
plt.xlabel('Minute of day')
plt.ylabel('Day of month')
plt.title('Percent of minute having missing data')
plt.colorbar()

# View percent of each 1 minute sample having missing data percent_missing = 100*zero_count/(spa*sample_rate) plt.figure(dpi=200, figsize = [9,3]) plt.imshow(percent_missing,aspect='auto',origin='lower') plt.xlabel('Minute of day') plt.ylabel('Day of month') plt.title('Percent of minute having missing data') plt.colorbar()

Out[15]:

<matplotlib.colorbar.Colorbar at 0x7a4c38583bb0>

This month was very well sampled, with nearly no missing data. Nearly every time bin for the month was completely sampled (0% missing), and the brief period having any missing data still had at least 70% coverage per time bin. These attributes do not motivate screening any of the CI results for quality control. However, analysis periods with greater amounts of missing data would motivate such screening.