{
    "nbformat": 4,
    "nbformat_minor": 5,
    "metadata": {
        "kernelspec": {
            "display_name": "Python",
            "language": "python",
            "name": "python"
        },
        "language_info": {
            "name": "python"
        }
    },
    "cells": [
        {
            "id": "0e549749",
            "cell_type": "markdown",
            "source": "<div class='alert alert-warning'>\n\nSciPy's interactive examples with Jupyterlite are experimental and may not always work as expected. Execution of cells containing imports may result in large downloads (up to 60MB of content for the first import from SciPy). Load times when importing from SciPy may take roughly 10-20 seconds. If you notice any problems, feel free to open an [issue](https://github.com/scipy/scipy/issues/new/choose).\n\n</div>",
            "metadata": {}
        },
        {
            "id": "01cd1068",
            "cell_type": "markdown",
            "source": "The following example plots the cross power spectral density of two signals with\nsome common features:\n",
            "metadata": {}
        },
        {
            "id": "259f22ef",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "import numpy as np\nfrom scipy import signal\nimport matplotlib.pyplot as plt\nrng = np.random.default_rng()\n\n# Generate two test signals with some common features:\nN, fs = 100_000, 10e3  # number of samples and sampling frequency\namp, freq = 20, 1234.0  # amplitude and frequency of utilized sine signal\nnoise_power = 0.001 * fs / 2\ntime = np.arange(N) / fs\nb, a = signal.butter(2, 0.25, 'low')\nx = rng.normal(scale=np.sqrt(noise_power), size=time.shape)\ny = signal.lfilter(b, a, x)\nx += amp*np.sin(2*np.pi*freq*time)\ny += rng.normal(scale=0.1*np.sqrt(noise_power), size=time.shape)\n\n# Compute and plot the magnitude of the cross spectral density:\nnperseg, noverlap, win = 1024, 512, 'hann'\nf, Pxy = signal.csd(x, y, fs, win, nperseg, noverlap)\nfig0, ax0 = plt.subplots(tight_layout=True)\nax0.set_title(f\"CSD ({win.title()}-window, {nperseg=}, {noverlap=})\")\nax0.set(xlabel=\"Frequency $f$ in kHz\", ylabel=\"CSD Magnitude in V²/Hz\")\nax0.semilogy(f/1e3, np.abs(Pxy))\nax0.grid()\nplt.show()",
            "outputs": []
        },
        {
            "id": "34657e40",
            "cell_type": "markdown",
            "source": "The cross spectral density is calculated by taking the average over the time slices\nof a spectrogram:\n",
            "metadata": {}
        },
        {
            "id": "7abb8bbb",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "SFT = signal.ShortTimeFFT.from_window('hann', fs, nperseg, noverlap,\n                                      scale_to='psd', fft_mode='onesided2X',\n                                      phase_shift=None)\nSxy1 = SFT.spectrogram(y, x, detr='constant', k_offset=nperseg//2,\n                       p0=0, p1=(N-noverlap) // SFT.hop)\nPxy1 = Sxy1.mean(axis=-1)\nnp.allclose(Pxy, Pxy1)  # same result as with csd()",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "True"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "6e2d5085",
            "cell_type": "markdown",
            "source": "As discussed in the Notes section, the results of using an approach analogous to\nthe code snippet above and the `csd` function may deviate due to implementation\ndetails.\n\nNote that the code snippet above can be easily adapted to determine other\nstatistical properties than the mean value.\n",
            "metadata": {}
        }
    ]
}