\n",
" \n",
" * Distributed under the terms of the GPL License\n",
" * Maintainer: yzhang@mbari.org\n",
" * Authors: Yanwu Zhang yzhang@mbari.org, Paul McGill mcgill@mbari.org, Danelle Cline (dcline@mbari.org), John Ryan ryjo@mbari.org"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "O-yXYo5vEEu8"
},
"source": [
"## Decimation of MARS hydrophone data\n",
"\n",
"An extensive (6+ years and growing) archive of sound recordings from a deep-sea location [along the eastern margin of the North Pacific Ocean](https://www.mbari.org/at-sea/cabled-observatory/) has been made available through AWS Open data. Temporal coverage of the recording archive has been 95% since project inception in July 2015. The original recordings have a sample rate of 256 kHz. Many research topics can be effectively studied using data with a lower sample rate, and this Open Data project includes daily files decimated to 2 kHz and 16 kHz. The purpose of this notebook is to illustrate a method of optimizing decimation processing, which is applicable to any desired sample rate. The demonstration uses Python, but the algorithms can be implemented in other languages.\n",
"\n",
"This method enables design of the optimal windowed-sinc anti-aliasing low-pass filter (using a certain window) that meets the desired passband and stopband requirements based on signal retention and exclusion needs. The best combination of the sinc function’s cutoff frequency and the main-lobe bandwidth of the window function generates the shortest qualifying filter.\n",
"\n",
"If you use this data set, please **[cite our project](https://ieeexplore.ieee.org/document/7761363).**\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "26Vtb9VwEEu9"
},
"source": [
"## Data Overview\n",
"The full-resolution audio data are in [WAV](https://en.wikipedia.org/wiki/WAV) format in s3 buckets named pacific-sound-256khz-yyyy, where yyyy is 2015 or later. Buckets are stored as objects, so the data aren't physically stored in folders or directories as you may be familiar with, but you can think of it conceptually as follows:\n",
"\n",
"```\n",
"pacific-sound-256khz-2021\n",
" |\n",
" individual 10-minute files\n",
"```\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "tc33lkaFEEu9"
},
"source": [
"## Install required dependencies\n",
"\n",
"First, let's install the required software dependencies. \n",
"\n",
"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.\n",
"\n",
"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](https://github.com/mbari-org/pacific-sound-notebooks/) - this has all the dependencies that are needed. "
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"id": "-8PsxkDAEEu-",
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "9f095ebf-8445-4310-e0c4-3cd4833a58ae"
},
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"\u001b[K |████████████████████████████████| 132 kB 4.7 MB/s \n",
"\u001b[K |████████████████████████████████| 79 kB 5.3 MB/s \n",
"\u001b[K |████████████████████████████████| 9.2 MB 40.3 MB/s \n",
"\u001b[K |████████████████████████████████| 140 kB 20.1 MB/s \n",
"\u001b[31mERROR: pip's dependency resolver does not currently take into account all the packages that are installed. This behaviour is the source of the following dependency conflicts.\n",
"requests 2.23.0 requires urllib3!=1.25.0,!=1.25.1,<1.26,>=1.21.1, but you have urllib3 1.26.12 which is incompatible.\u001b[0m\n",
"\u001b[?25h"
]
}
],
"source": [
"!pip install -q boto3 --quiet\n",
"!pip install -q soundfile --quiet\n",
"!pip install -q scipy --quiet\n",
"!pip install -q numpy --quiet\n",
"!pip install -q matplotlib --quiet"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "S5Cn498hEEu-"
},
"source": [
"### Import all packages"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"id": "BFfThmwIEEu-"
},
"outputs": [],
"source": [
"import boto3\n",
"from botocore import UNSIGNED\n",
"from botocore.client import Config\n",
"from six.moves.urllib.request import urlopen\n",
"import io\n",
"import math\n",
"import scipy\n",
"from scipy import signal, interpolate\n",
"import numpy as np\n",
"import soundfile as sf\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "F5KeMFHaEEu_"
},
"source": [
"## List the contents of a monthly directory"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"id": "1PnE5fX8EEu_"
},
"outputs": [],
"source": [
"s3 = boto3.client('s3',\n",
" aws_access_key_id='',\n",
" aws_secret_access_key='', \n",
" config=Config(signature_version=UNSIGNED))"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"id": "gewUkRQbEEu_",
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "6a29511b-2497-4eee-e663-002d1aafee04"
},
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"01/MARS_20210101_000424.wav\n",
"01/MARS_20210101_001424.wav\n",
"01/MARS_20210101_002424.wav\n",
"01/MARS_20210101_003424.wav\n",
"01/MARS_20210101_004424.wav\n",
"01/MARS_20210101_005424.wav\n",
"01/MARS_20210101_010424.wav\n",
"01/MARS_20210101_011424.wav\n",
"01/MARS_20210101_012424.wav\n",
"01/MARS_20210101_013424.wav\n",
"01/MARS_20210101_014424.wav\n",
"01/MARS_20210101_015424.wav\n",
"01/MARS_20210101_020424.wav\n",
"01/MARS_20210101_021424.wav\n",
"01/MARS_20210101_022424.wav\n",
"01/MARS_20210101_023424.wav\n",
"01/MARS_20210101_024424.wav\n",
"01/MARS_20210101_025424.wav\n",
"01/MARS_20210101_030424.wav\n",
"01/MARS_20210101_031424.wav\n",
"01/MARS_20210101_032424.wav\n",
"01/MARS_20210101_033424.wav\n"
]
}
],
"source": [
"bucket = 'pacific-sound-256khz-2021'\n",
"\n",
"for i, obj in enumerate(s3.list_objects_v2(Bucket=bucket)['Contents']):\n",
" print(obj['Key'])\n",
" if i > 20:\n",
" break"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "-WSJMl9QEEvA"
},
"source": [
"## Read metadata from a file"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"id": "jc3SWLX5EEvA",
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "3e420614-3eb2-40ff-e7b1-66d6be4fe396"
},
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"Reading metadata from https://pacific-sound-256khz-2021.s3.amazonaws.com/09/MARS_20210915_070829.wav\n"
]
},
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<_io.BytesIO object at 0x7fcdb059ea10>\n",
"samplerate: 256000 Hz\n",
"channels: 1\n",
"duration: 218 samples\n",
"format: WAV (Microsoft) [WAV]\n",
"subtype: Signed 24 bit PCM [PCM_24]\n",
"endian: FILE\n",
"sections: 1\n",
"frames: 218\n",
"extra_info: \"\"\"\n",
" Length : 1000\n",
" RIFF : 460800336 (should be 992)\n",
" WAVE\n",
" LIST : 292\n",
" INFO\n",
" IART : icListen HF #1689\n",
" IPRD : RB9-ETH R4\n",
" ICRD : 2021-09-15T07:08:29+00\n",
" ISFT : Lucy V4.3.6\n",
" INAM : MARS_20210915_070829.wav\n",
" ICMT : 3.000000 V pk, -177 dBV re 1uPa, 44.0 % RH, 6.0 deg C, 8388608 = Max Count\n",
" fmt : 16\n",
" Format : 0x1 => WAVE_FORMAT_PCM\n",
" Channels : 1\n",
" Sample Rate : 256000\n",
" Block Align : 3\n",
" Bit Width : 24\n",
" Bytes/sec : 768000\n",
" data : 460800000 (should be 656)\n",
" End\n",
" \"\"\""
]
},
"metadata": {},
"execution_count": 5
}
],
"source": [
"bucket = 'pacific-sound-256khz-2021'\n",
"filename = '09/MARS_20210915_070829.wav'\n",
"url = f'https://{bucket}.s3.amazonaws.com/{filename}'\n",
"print(f'Reading metadata from {url}')\n",
"sf.info(io.BytesIO(urlopen(url).read(1000)), verbose=True) "
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "3wCE630CEEvA"
},
"source": [
"## Read data from a file"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"id": "ddx_VLpJEEvA",
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "5ddd0ddd-d9b1-45ec-e969-e6b1fcd9795c"
},
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"Reading data from https://pacific-sound-256khz-2021.s3.amazonaws.com/09/MARS_20210915_070829.wav\n"
]
}
],
"source": [
"print(f'Reading data from {url}')\n",
"x, sample_rate = sf.read(io.BytesIO(urlopen(url).read()),dtype='float32')"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "BJG91K_aEEvB"
},
"source": [
"## Produce a spectrogram"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "X_g0Ou26EEvB"
},
"source": [
"### Calibrated Spectrum Levels"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "guJX5TBEEEvB"
},
"source": [
"### Calibration metadata\n",
"Frequency-dependent hydrophone sensitivity data are defined in the following files, one for each deployment:\n",
"* https://bitbucket.org/mbari/pacific-sound/src/master/MBARI_MARS_Hydrophone_Deployment01.json\n",
"* https://bitbucket.org/mbari/pacific-sound/src/master/MBARI_MARS_Hydrophone_Deployment02.json"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "A0m2-J7CEEvB"
},
"source": [
"### Compute spectrogram\n",
"\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"id": "BUynuK_OEEvB",
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "2f807fef-5f6f-4916-ce37-813353cd0f51"
},
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"600 segments of length 1 seconds in 600.0 seconds of audio\n"
]
}
],
"source": [
"# convert scaled voltage to volts\n",
"v = x*3 \n",
"nsec = (v.size)/sample_rate # number of seconds in vector\n",
"spa = 1 # seconds per average\n",
"nseg = int(nsec/spa)\n",
"print(f'{nseg} segments of length {spa} seconds in {nsec} seconds of audio')"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"id": "FdpXH4xUEEvB",
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "6daa1252-1c56-449f-85b6-fe7ef380cc28"
},
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"262144\n"
]
},
{
"output_type": "execute_result",
"data": {
"text/plain": [
"(131073, 600)"
]
},
"metadata": {},
"execution_count": 8
}
],
"source": [
"lenfft_input = 2**int(np.ceil(np.log2(sample_rate)))\n",
"print(lenfft_input)\n",
"#\n",
"# initialize empty LTSA\n",
"nfreq = int(lenfft_input/2+1)\n",
"sg_input = np.empty((nfreq, nseg), float)\n",
"sg_input.shape"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"id": "NlBpRSwtEEvC",
"colab": {
"base_uri": "https://localhost:8080/"
},
"outputId": "335908c3-9562-448a-e97f-39e1eab87ad2"
},
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"Comparing power of the input signal computed in the time domain versus that computed in the frequency domain:\n",
"0.00012095761\n",
"0.000121922813689169\n"
]
}
],
"source": [
"# get window\n",
"w_input = scipy.signal.get_window('blackman',sample_rate)\n",
"window_correction = np.mean(np.square(w_input))\n",
"#\n",
"numDataPoints_input = int(sample_rate*spa)\n",
"#\n",
"Ind_keep = np.arange(0, int(lenfft_input/2)+1)\n",
"f_input = (Ind_keep/lenfft_input) * sample_rate\n",
"#\n",
"# Calculate spectrogram\n",
"for x in range(0,nseg):\n",
" cstart = x*spa*sample_rate\n",
" cend = (x+1)*spa*sample_rate\n",
"# f,psd = scipy.signal.welch(v[cstart:cend], fs=sample_rate, window=w_input, nfft=sample_rate)\n",
" psd_input = np.square(np.absolute(np.fft.fft(np.multiply(v[cstart:cend],w_input), n=lenfft_input)))/(numDataPoints_input*window_correction)/sample_rate\n",
" psd_input_log10 = 10*np.log10(psd_input[Ind_keep])\n",
" sg_input[:,x] = psd_input_log10\n",
" if (x == 0):\n",
" psd_input_check = psd_input\n",
" print(\"Comparing power of the input signal computed in the time domain versus that computed in the frequency domain:\")\n",
" print(np.mean(np.square(v[cstart:cend])))\n",
" print(sum(psd_input_check)*(sample_rate/lenfft_input))"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "06SkIb43EEvC"
},
"source": [
"### Apply calibration\n",
"Frequency-dependent hydrophone sensitivity data are reported in the json files identified above. This example file is from the second hydrophone deployment, for which the calibration data are manually entered below. Note that the lowest measured value, at 250 Hz, is assumed to cover lower frequencies and repeated as a value at 0 Hz to allow interpolation to the spectrogram output frequencies across the full frequency range. \n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"id": "X9J_0gfgEEvC",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 1000
},
"outputId": "a1b8a469-264a-4aab-f664-a5d341f6b874"
},
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
""
]
},
"metadata": {},
"execution_count": 10
},
{
"output_type": "display_data",
"data": {
"text/plain": [
"