{
    "nbformat": 4,
    "nbformat_minor": 5,
    "metadata": {
        "kernelspec": {
            "display_name": "Python",
            "language": "python",
            "name": "python"
        },
        "language_info": {
            "name": "python"
        }
    },
    "cells": [
        {
            "id": "74a95e8a",
            "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": "d1057f21",
            "cell_type": "markdown",
            "source": "The following example depicts a signal being up-sampled from 20 samples to 100\nsamples. The ringing at the beginning of the up-sampled signal is due to\ninterpreting the signal being periodic. The red square in the plot illustrates that\nperiodictiy by showing the first sample of the next cycle of the signal.\n",
            "metadata": {}
        },
        {
            "id": "ff61ff79",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "import numpy as np\nimport matplotlib.pyplot as plt\nfrom scipy.signal import resample\n\nn0, n1 = 20, 100  # number of samples\nt0 = np.linspace(0, 10, n0, endpoint=False)  # input time stamps\nx0 = np.cos(-t0**2/6)  # input signal\n\nx1 = resample(x0, n1)  # resampled signal\nt1 = np.linspace(0, 10, n1, endpoint=False)  # timestamps of x1\n\nfig0, ax0 = plt.subplots(1, 1, tight_layout=True)\nax0.set_title(f\"Resampling $x(t)$ from {n0} samples to {n1} samples\")\nax0.set(xlabel=\"Time $t$\", ylabel=\"Amplitude $x(t)$\")\nax0.plot(t1, x1, '.-', alpha=.5, label=f\"Resampled\")\nax0.plot(t0, x0, 'o-', alpha=.5, label=\"Original\")\nax0.plot(10, x0[0], 'rs', alpha=.5, label=\"Next Cycle\")\nax0.legend(loc='best')\nax0.grid(True)\nplt.show()",
            "outputs": []
        },
        {
            "id": "a5f2d69f",
            "cell_type": "markdown",
            "source": "The following example compares this function with a naive `~scipy.fft.rfft` /\n`~scipy.fft.irfft` combination: An input signal with a sampling interval of one\nsecond is upsampled by a factor of eight. The first figure depicts an odd number of\ninput samples whereas the second figure an even number. The upper subplots show the\nsignals over time: The input samples are marked by large green dots, the upsampled\nsignals by a continuous and a dashed line. The lower subplots show the magnitude\nspectrum: The FFT values of the input are depicted by large green dots, which lie\nin the frequency interval [-0.5, 0.5] Hz, whereas the frequency interval of the\nupsampled signal is [-4, 4] Hz. The continuous green line depicts the upsampled\nspectrum without antialiasing filter, which is a periodic continuation of the input\nspectrum. The blue x's and orange dots depict the FFT values of the signal created\nby the naive approach as well as this function's result.\n",
            "metadata": {}
        },
        {
            "id": "821abca8",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "import matplotlib.pyplot as plt\nimport numpy as np\nfrom scipy.fft import fftshift, fftfreq, fft, rfft, irfft\nfrom scipy.signal import resample, resample_poly\n\nfac, T0, T1 = 8, 1, 1/8  # upsampling factor and sampling intervals\nfor n0 in (15, 16):  # number of samples of input signal\n    n1 = fac * n0  # number of samples of upsampled signal\n    t0, t1 = T0 * np.arange(n0), T1 * np.arange(n1)  # time stamps\n    x0 = np.zeros(n0)  # input signal has two non-zero sample values\n    x0[n0//2], x0[n0//2+1] = n0 // 2, -(n0 // 2)\n\n    x1n = irfft(rfft(x0), n=n1) * n1 / n0  # naive resampling\n    x1r = resample(x0, n1)  # resample signal\n\n    # Determine magnitude spectrum:\n    x0_up = np.zeros_like(x1r)  # upsampling without antialiasing filter\n    x0_up[::n1 // n0] = x0\n    X0, X0_up = (fftshift(fft(x_)) / n0 for x_ in (x0, x0_up))\n    XX1 = (fftshift(fft(x_)) / n1 for x_ in (x1n, x1r))\n    f0, f1 = fftshift(fftfreq(n0, T0)), fftshift(fftfreq(n1, T1))  # frequencies\n    df = f0[1] - f0[0]  # frequency resolution\n\n    fig, (ax0, ax1) = plt.subplots(2, 1, layout='constrained', figsize=(5, 4))\n    ax0.set_title(rf\"Upsampling ${fac}\\times$ from {n0} to {n1} samples\")\n    ax0.set(xlabel=\"Time $t$ in seconds\", ylabel=\"Amplitude $x(t)$\", \n            xlim=(0, n1*T1))\n    ax0.step(t0, x0, 'C2o-', where='post', alpha=.3, linewidth=2, \n             label=\"$x_0(t)$ / $X_0(f)$\")\n    for x_, l_ in zip((x1n, x1r), ('C0--', 'C1-')):\n        ax0.plot(t1, x_, l_, alpha=.5, label=None)\n    ax0.grid()\n    ax1.set(xlabel=rf\"Frequency $f$ in hertz ($\\Delta f = {df*1e3:.1f}\\,$mHz)\", \n            ylabel=\"Magnitude $|X(f)|$\", xlim=(-0.7, 0.7))\n    ax1.axvspan(0.5/T0, f1[-1], color='gray', alpha=.2)\n    ax1.axvspan(f1[0], -0.5/T0, color='gray', alpha=.2)\n    ax1.plot(f1, abs(X0_up), 'C2-', f0, abs(X0),  'C2o', alpha=.3, linewidth=2)\n    for X_, n_, l_ in zip(XX1, (\"naive\", \"resample\"), ('C0x--', 'C1.-')): \n        ax1.plot(f1, abs(X_), l_, alpha=.5, label=n_)\n    ax1.grid()\n    fig.legend(loc='outside lower center', ncols=4)    \nplt.show()",
            "outputs": []
        },
        {
            "id": "79b20388",
            "cell_type": "markdown",
            "source": "The first figure shows that upsampling an odd number of samples produces identical\nresults. The second figure illustrates that the signal produced with the naive\napproach (dashed blue line) from an even number of samples does not touch all\noriginal samples. This deviation is due to `resample` correctly treating unpaired\nfrequency bins. I.e., the input `x1` has a bin pair ±0.5 Hz, whereas the output has\nonly one unpaired bin at -0.5 Hz, which demands rescaling of that bin pair.\nGenerally, special treatment is required if ``n_x != num`` and ``min(n_x, num)`` is\neven. If the bin values at `±m` are zero, obviously, no special treatment is\nneeded. Consult the source code of `resample` for details.\n\nThe final example shows how to utilize `resample_poly` to speed up the\ndown-sampling: The input signal a non-zero value at $t=0$ and is downsampled\nfrom 19937 to 128 samples. Since 19937 is prime, the FFT is expected to be slow. To\nspeed matters up, `resample_poly` is used to downsample first by a factor of ``n0\n// n1 = 155`` and then pass the result to `resample`. Two parameterization of \n`resample_poly` are used: Passing ``padtype='wrap'`` treats the input as being\nperiodic wheras the default parametrization performs zero-padding. The upper\nsubplot shows the resulting signals over time whereas the lower subplot depicts the\nresulting one-sided magnitude spectra.\n",
            "metadata": {}
        },
        {
            "id": "de3a2879",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "import matplotlib.pyplot as plt\nimport numpy as np\nfrom scipy.fft import rfftfreq, rfft\nfrom scipy.signal import resample, resample_poly\n\nn0 = 19937 # number of input samples - prime\nn1 = 128  # number of output samples - fast FFT length\nT0, T1 = 1/n0, 1/n1  # sampling intervals\nt0, t1 = np.arange(n0)*T0, np.arange(n1)*T1  # time stamps\n\nx0 = np.zeros(n0)  # Input has one non-zero sample\nx0[0] = n0\n\nx1r = resample(x0, n1)  # slow due to n0 being prime\n# This is faster:\nx1p = resample(resample_poly(x0, 1, n0 // n1, padtype='wrap'), n1)  # periodic \nx2p = resample(resample_poly(x0, 1, n0 // n1), n1)  # with zero-padding\n\nX0 = rfft(x0) / n0 \nX1r, X1p, X2p = rfft(x1r) / n1, rfft(x1p) / n1, rfft(x2p) / n1\nf0, f1 = rfftfreq(n0, T0), rfftfreq(n1, T1)\n\nfig, (ax0, ax1) = plt.subplots(2, 1, layout='constrained', figsize=(5, 4))\nax0.set_title(f\"Dowsampled Impulse response (from {n0} to {n1} samples)\")\nax0.set(xlabel=\"Time $t$ in seconds\", ylabel=\"Amplitude $x(t)$\", xlim=(-T1, 1)) \nfor x_ in (x1r, x1p, x2p):\n    ax0.plot(t1, x_, alpha=.5)\nax0.grid()\nax1.set(xlabel=rf\"Frequency $f$ in hertz ($\\Delta f = {f1[1]}\\,$Hz)\", \n        ylabel=\"Magnitude $|X(f)|$\", xlim=(0, 0.55/T1))\nax1.axvspan(0.5/T1, f0[-1], color='gray', alpha=.2)\nax1.plot(f1, abs(X1r), 'C0.-', alpha=.5, label=\"resample\")\nax1.plot(f1, abs(X1p), 'C1.-', alpha=.5, label=\"resample_poly(padtype='wrap')\")\nax1.plot(f1, abs(X2p), 'C2x-', alpha=.5, label=\"resample_poly\")\nax1.grid()\nfig.legend(loc='outside lower center', ncols=2)\nplt.show()    ",
            "outputs": []
        },
        {
            "id": "298cebc9",
            "cell_type": "markdown",
            "source": "The plots show that the results of the \"pure\" `resample` and the usage of the\ndefault parameters of `resample_poly` agree well.  The periodic padding of\n`resample_poly` (``padtype='wrap'``) on the other hand produces significant\ndeviations. This is caused by the disconiuity at the beginning of the signal, for\nwhich the default filter of `resample_poly` is not suited well. This example\nillustrates that for some use cases, adpating the `resample_poly` parameters may\nbe beneficial. `resample` has a big advantage in this regard: It uses the ideal\nantialiasing filter with the maximum bandwidth by default.\n\nNote that the doubled spectral magnitude at the Nyqist frequency of 64 Hz is due the\neven number of ``n1=128`` output samples, which requires a special treatment as \ndiscussed in the previous example. \n",
            "metadata": {}
        }
    ]
}