raab/GP2023_chirp_detection
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

added notebook

Browse Source
This commit is contained in:
weygoldt 2023-04-13 15:14:10 +02:00
parent 95a256b517
commit c4f05d9891
No known key found for this signature in database
2 changed files with 397 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

389
chirp_instantaneous_freq/chirp_exploration.ipynb Normal file
View File

@ -0,0 +1,389 @@
{
"cells": [
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"# Why can the instantaneous frequency of a band-pass filtered chirp recording go down ...\n",
"... if a chirp is an up-modulation of the frequency? \n",
"\n",
"This is the question we try to answer in this notebook"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"QApplication: invalid style override passed, ignoring it.\n",
" Available styles: Windows, Fusion\n"
]
}
],
"source": [
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"import thunderfish.fakefish as ff \n",
"from filters import instantaneous_frequency, bandpass_filter\n",
"%matplotlib qt\n",
"\n",
"# parameters that stay the same\n",
"samplerate = 20000\n",
"duration = 0.2\n",
"chirp_freq = 5\n",
"smooth = 3"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"qt.qpa.wayland: Wayland does not support QWindow::requestActivate()\n"
]
}
],
"source": [
"def make_chirp(eodf, size, width, kurtosis, contrast, phase0):\n",
"\n",
" chirp_trace, amp = ff.chirps(\n",
" eodf = eodf,\n",
" samplerate = samplerate,\n",
" duration = duration,\n",
" chirp_freq = chirp_freq,\n",
" chirp_size = size,\n",
" chirp_width = width,\n",
" chirp_kurtosis = kurtosis,\n",
" chirp_contrast = contrast,\n",
" )\n",
"\n",
" chirp = ff.wavefish_eods(\n",
" fish = 'Alepto',\n",
" frequency = chirp_trace,\n",
" samplerate = samplerate,\n",
" duration = duration,\n",
" phase0 = phase0,\n",
" noise_std = 0,\n",
" )\n",
"\n",
" chirp *= amp\n",
"\n",
" return chirp_trace, chirp\n",
"\n",
"def filtered_chirp(eodf, size, width, kurtosis, contrast, phase0):\n",
"\n",
" time = np.arange(0, duration, 1/samplerate)\n",
" chirp_trace, chirp = make_chirp(\n",
" eodf = eodf, \n",
" size = size, \n",
" width = width, \n",
" kurtosis = kurtosis, \n",
" contrast = contrast, \n",
" phase0 = phase0,\n",
" )\n",
"\n",
" # apply filters\n",
" narrow_filtered = bandpass_filter(chirp, samplerate, eodf-10, eodf+10)\n",
" narrow_freqtime, narrow_freq = instantaneous_frequency(narrow_filtered, samplerate, smooth)\n",
" broad_filtered = bandpass_filter(chirp, samplerate, eodf-300, eodf+300)\n",
" broad_freqtime, broad_freq = instantaneous_frequency(broad_filtered, samplerate, smooth)\n",
"\n",
" original = (time, chirp_trace, chirp)\n",
" broad = (broad_freqtime, broad_freq, broad_filtered)\n",
" narrow = (narrow_freqtime, narrow_freq, narrow_filtered)\n",
"\n",
" return original, broad, narrow\n",
"\n",
"def plot(original, broad, narrow, axs):\n",
"\n",
" axs[0].plot(original[0], original[1], label='chirp trace')\n",
" axs[0].plot(broad[0], broad[1], label='broad filtered')\n",
" axs[0].plot(narrow[0], narrow[1], label='narrow filtered')\n",
" axs[1].plot(original[0], original[2], label='unfiltered')\n",
" axs[1].plot(original[0], broad[2], label='broad filtered')\n",
" axs[1].plot(original[0], narrow[2], label='narrow filtered')\n",
"\n",
"original, broad, narrow = filtered_chirp(600, 100, 0.02, 1, 0.1, 0)\n",
"fig, axs = plt.subplots(2, 1, figsize=(10, 5), sharex=True)\n",
"plot(original, broad, narrow, axs)\n",
"fig.align_labels()\n",
"plt.show()"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## Chirp size\n",
"now that we have established an easy way to simulate and plot the chirps, lets change the chirp size and see how the narrow-filtered instantaneous frequency changes."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"size 10 Hz; Integral 0.117\n",
"size 30 Hz; Integral 0.35\n",
"size 50 Hz; Integral 0.584\n",
"size 70 Hz; Integral 0.818\n",
"size 90 Hz; Integral 1.051\n",
"size 110 Hz; Integral 1.285\n",
"size 130 Hz; Integral 1.518\n",
"size 150 Hz; Integral 1.752\n",
"size 170 Hz; Integral 1.986\n",
"size 190 Hz; Integral 2.219\n",
"size 210 Hz; Integral 2.453\n",
"size 230 Hz; Integral 2.687\n",
"size 250 Hz; Integral 2.92\n",
"size 270 Hz; Integral 3.154\n",
"size 290 Hz; Integral 3.387\n",
"size 310 Hz; Integral 3.621\n",
"size 330 Hz; Integral 3.855\n",
"size 350 Hz; Integral 4.088\n",
"size 370 Hz; Integral 4.322\n",
"size 390 Hz; Integral 4.555\n",
"size 410 Hz; Integral 4.789\n",
"size 430 Hz; Integral 5.023\n",
"size 450 Hz; Integral 5.256\n",
"size 470 Hz; Integral 5.49\n",
"size 490 Hz; Integral 5.724\n",
"size 510 Hz; Integral 5.957\n",
"size 530 Hz; Integral 6.191\n",
"size 550 Hz; Integral 6.424\n",
"size 570 Hz; Integral 6.658\n",
"size 590 Hz; Integral 6.892\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"qt.qpa.wayland: Wayland does not support QWindow::requestActivate()\n"
]
}
],
"source": [
"sizes = np.arange(10, 600, 20)\n",
"fig, axs = plt.subplots(2, len(sizes), figsize=(10, 5), sharex=True, sharey='row')\n",
"integrals = []\n",
"\n",
"for i, size in enumerate(sizes):\n",
" original, broad, narrow = filtered_chirp(600, size, 0.02, 1, 0.1, 0)\n",
"\n",
" integral = np.sum(original[1]-600)/(20000)\n",
" integrals.append(integral)\n",
"\n",
" plot(original, broad, narrow, axs[:, i])\n",
" axs[:, i][0].set_xlim(0.06, 0.14)\n",
" axs[0, i].set_title(np.round(integral, 3))\n",
" print(f'size {size} Hz; Integral {np.round(integral,3)}')\n",
" \n",
"fig.legend(handles=axs[0,0].get_lines(), loc='upper center', ncol=3)\n",
"axs[0,0].set_ylabel('frequency [Hz]')\n",
"axs[1,0].set_ylabel('amplitude [a.u.]')\n",
"fig.supxlabel('time [s]')\n",
"fig.align_labels()\n",
"plt.show()"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## Chirp width"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"widths = np.arange(0.02, 0.08, 0.005)\n",
"fig, axs = plt.subplots(2, len(widths), figsize=(10, 5), sharex=True, sharey='row')\n",
"integrals = []\n",
"\n",
"for i, width in enumerate(widths):\n",
" if i > 9:\n",
" break\n",
"\n",
" original, broad, narrow = filtered_chirp(600, 100, width, 1, 0.1, 0)\n",
"\n",
" integral = np.sum(original[1]-600)/(20000)\n",
"\n",
" plot(original, broad, narrow, axs[:, i])\n",
" axs[:, i][0].set_xlim(0.06, 0.14)\n",
" axs[0, i].set_title(f'width {np.round(width, 2)} s')\n",
" print(f'width {width} s; Integral {np.round(integral, 3)}')\n",
" \n",
"fig.legend(handles=axs[0,0].get_lines(), loc='upper center', ncol=3)\n",
"axs[0,0].set_ylabel('frequency [Hz]')\n",
"axs[1,0].set_ylabel('amplitude [a.u.]')\n",
"fig.supxlabel('time [s]')\n",
"fig.align_labels()\n",
"plt.show()"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## Chirp kurtosis"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"kurtosiss = np.arange(0, 20, 1.6)\n",
"fig, axs = plt.subplots(2, len(kurtosiss), figsize=(10, 5), sharex=True, sharey='row')\n",
"integrals = []\n",
"\n",
"for i, kurtosis in enumerate(kurtosiss):\n",
"\n",
" original, broad, narrow = filtered_chirp(600, 100, 0.02, kurtosis, 0.1, 0)\n",
"\n",
" integral = np.sum(original[1]-600)/(20000)\n",
"\n",
" plot(original, broad, narrow, axs[:, i])\n",
" axs[:, i][0].set_xlim(0.06, 0.14)\n",
" axs[0, i].set_title(f'kurt {np.round(kurtosis, 2)}')\n",
" print(f'kurt {kurtosis}; Integral {np.round(integral, 3)}')\n",
" \n",
"fig.legend(handles=axs[0,0].get_lines(), loc='upper center', ncol=3)\n",
"axs[0,0].set_ylabel('frequency [Hz]')\n",
"axs[1,0].set_ylabel('amplitude [a.u.]')\n",
"fig.supxlabel('time [s]')\n",
"fig.align_labels()\n",
"plt.show()"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## Chirp contrast"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"contrasts = np.arange(0.0, 1.1, 0.1)\n",
"fig, axs = plt.subplots(2, len(sizes), figsize=(10, 5), sharex=True, sharey='row')\n",
"integrals = []\n",
"\n",
"for i, contrast in enumerate(contrasts):\n",
" if i > 9:\n",
" break\n",
" original, broad, narrow = filtered_chirp(600, 100, 0.02, 1, contrast, 0)\n",
"\n",
" integral = np.trapz(original[2], original[0])\n",
" integrals.append(integral)\n",
"\n",
" plot(original, broad, narrow, axs[:, i])\n",
" axs[:, i][0].set_xlim(0.06, 0.14)\n",
" axs[0, i].set_title(f'contr {np.round(contrast, 2)}')\n",
" \n",
"fig.legend(handles=axs[0,0].get_lines(), loc='upper center', ncol=3)\n",
"axs[0,0].set_ylabel('frequency [Hz]')\n",
"axs[1,0].set_ylabel('amplitude [a.u.]')\n",
"fig.supxlabel('time [s]')\n",
"fig.align_labels()\n",
"plt.show()"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## Chirp phase "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"phases = np.arange(0.0, 2 * np.pi, 0.2)\n",
"fig, axs = plt.subplots(2, len(sizes), figsize=(10, 5), sharex=True, sharey='row')\n",
"integrals = []\n",
"for i, phase in enumerate(phases):\n",
" if i > 9:\n",
" break\n",
"\n",
" original, broad, narrow = filtered_chirp(600, 100, 0.02, 1, 0.1, phase)\n",
"\n",
" integral = np.trapz(original[2], original[0])\n",
" integrals.append(integral)\n",
"\n",
" plot(original, broad, narrow, axs[:, i])\n",
" axs[:, i][0].set_xlim(0.06, 0.14)\n",
" axs[0, i].set_title(f'phase {np.round(phase, 2)}')\n",
"\n",
" \n",
"fig.legend(handles=axs[0,0].get_lines(), loc='upper center', ncol=3)\n",
"axs[0,0].set_ylabel('frequency [Hz]')\n",
"axs[1,0].set_ylabel('amplitude [a.u.]')\n",
"fig.supxlabel('time [s]')\n",
"fig.align_labels()\n",
"plt.show()"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"These experiments show, that the narrow filtered instantaneous freuqency only switches its sign, when the integral of the instantaneous frequency (that was used to make the signal)\n",
"changes. Specifically, when the instantaneous frequency is 0.57, 1.57, 2.57 etc., the sign swithes. "
]
}
],
"metadata": {
"kernelspec": {
"display_name": "chirpdetection",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.2"
},
"orig_nbformat": 4
},
"nbformat": 4,
"nbformat_minor": 2
}

14
chirp_instantaneous_freq/test_parameters.py
View File

@ -1,7 +1,7 @@
import numpy as np
import matplotlib.pyplot as plt
import numpy as np
from filters import bandpass_filter, inst_freq, instantaneous_frequency
from fish_signal import chirps, wavefish_eods
from filters import bandpass_filter, instantaneous_frequency, inst_freq
from IPython import embed
@ -28,13 +28,14 @@ def extract_dict(dict, index):
return {key: value[index] for key, value in dict.items()}
def main(test1, test2, resolution=10):
def test(test1, test2, resolution=10):
assert test1 in [
"width",
"size",
"kurtosis",
"contrast",
], "Test1 not recognized"
assert test2 in [
"width",
"size",
@ -139,10 +140,11 @@ def main(test1, test2, resolution=10):
iter0 += 1
fig, ax = plt.subplots()
ax.imshow(distances, cmap="jet")
plt.show()
def main():
test("contrast", "kurtosis")
if __name__ == "__main__":
main("width", "size")
main()