From c4f05d989123ae02da3eb406b433d40a5a950d5a Mon Sep 17 00:00:00 2001 From: weygoldt <88969563+weygoldt@users.noreply.github.com> Date: Thu, 13 Apr 2023 15:14:10 +0200 Subject: [PATCH] added notebook --- .../chirp_exploration.ipynb | 389 ++++++++++++++++++ chirp_instantaneous_freq/test_parameters.py | 14 +- 2 files changed, 397 insertions(+), 6 deletions(-) create mode 100644 chirp_instantaneous_freq/chirp_exploration.ipynb diff --git a/chirp_instantaneous_freq/chirp_exploration.ipynb b/chirp_instantaneous_freq/chirp_exploration.ipynb new file mode 100644 index 0000000..8c27cf3 --- /dev/null +++ b/chirp_instantaneous_freq/chirp_exploration.ipynb @@ -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 +} diff --git a/chirp_instantaneous_freq/test_parameters.py b/chirp_instantaneous_freq/test_parameters.py index bad1e45..7348496 100644 --- a/chirp_instantaneous_freq/test_parameters.py +++ b/chirp_instantaneous_freq/test_parameters.py @@ -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()