390 lines
12 KiB
Plaintext
390 lines
12 KiB
Plaintext
{
|
|
"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
|
|
}
|