{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Power spectral density (PSD)\n\nPlotting power spectral density (PSD) using `~.Axes.psd`.\n\nThe PSD is a common plot in the field of signal processing. NumPy has\nmany useful libraries for computing a PSD. Below we demo a few examples\nof how this can be accomplished and visualized with Matplotlib.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nimport numpy as np\n\nimport matplotlib.mlab as mlab\n\n# Fixing random state for reproducibility\nnp.random.seed(19680801)\n\ndt = 0.01\nt = np.arange(0, 10, dt)\nnse = np.random.randn(len(t))\nr = np.exp(-t / 0.05)\n\ncnse = np.convolve(nse, r) * dt\ncnse = cnse[:len(t)]\ns = 0.1 * np.sin(2 * np.pi * t) + cnse\n\nfig, (ax0, ax1) = plt.subplots(2, 1, layout='constrained')\nax0.plot(t, s)\nax0.set_xlabel('Time (s)')\nax0.set_ylabel('Signal')\nax1.psd(s, NFFT=512, Fs=1 / dt)\n\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Compare this with the equivalent Matlab code to accomplish the same thing::\n\n    dt = 0.01;\n    t = [0:dt:10];\n    nse = randn(size(t));\n    r = exp(-t/0.05);\n    cnse = conv(nse, r)*dt;\n    cnse = cnse(1:length(t));\n    s = 0.1*sin(2*pi*t) + cnse;\n\n    subplot(211)\n    plot(t, s)\n    subplot(212)\n    psd(s, 512, 1/dt)\n\nBelow we'll show a slightly more complex example that demonstrates\nhow padding affects the resulting PSD.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "dt = np.pi / 100.\nfs = 1. / dt\nt = np.arange(0, 8, dt)\ny = 10. * np.sin(2 * np.pi * 4 * t) + 5. * np.sin(2 * np.pi * 4.25 * t)\ny = y + np.random.randn(*t.shape)\n\n# Plot the raw time series\nfig, axs = plt.subplot_mosaic([\n    ['signal', 'signal', 'signal'],\n    ['zero padding', 'block size', 'overlap'],\n], layout='constrained')\n\naxs['signal'].plot(t, y)\naxs['signal'].set_xlabel('Time (s)')\naxs['signal'].set_ylabel('Signal')\n\n# Plot the PSD with different amounts of zero padding. This uses the entire\n# time series at once\naxs['zero padding'].psd(y, NFFT=len(t), pad_to=len(t), Fs=fs)\naxs['zero padding'].psd(y, NFFT=len(t), pad_to=len(t) * 2, Fs=fs)\naxs['zero padding'].psd(y, NFFT=len(t), pad_to=len(t) * 4, Fs=fs)\n\n# Plot the PSD with different block sizes, Zero pad to the length of the\n# original data sequence.\naxs['block size'].psd(y, NFFT=len(t), pad_to=len(t), Fs=fs)\naxs['block size'].psd(y, NFFT=len(t) // 2, pad_to=len(t), Fs=fs)\naxs['block size'].psd(y, NFFT=len(t) // 4, pad_to=len(t), Fs=fs)\naxs['block size'].set_ylabel('')\n\n# Plot the PSD with different amounts of overlap between blocks\naxs['overlap'].psd(y, NFFT=len(t) // 2, pad_to=len(t), noverlap=0, Fs=fs)\naxs['overlap'].psd(y, NFFT=len(t) // 2, pad_to=len(t),\n                   noverlap=int(0.025 * len(t)), Fs=fs)\naxs['overlap'].psd(y, NFFT=len(t) // 2, pad_to=len(t),\n                   noverlap=int(0.1 * len(t)), Fs=fs)\naxs['overlap'].set_ylabel('')\naxs['overlap'].set_title('overlap')\n\nfor title, ax in axs.items():\n    if title == 'signal':\n        continue\n\n    ax.set_title(title)\n    ax.sharex(axs['zero padding'])\n    ax.sharey(axs['zero padding'])\n\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This is a ported version of a MATLAB example from the signal\nprocessing toolbox that showed some difference at one time between\nMatplotlib's and MATLAB's scaling of the PSD.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fs = 1000\nt = np.linspace(0, 0.3, 301)\nA = np.array([2, 8]).reshape(-1, 1)\nf = np.array([150, 140]).reshape(-1, 1)\nxn = (A * np.sin(2 * np.pi * f * t)).sum(axis=0)\nxn += 5 * np.random.randn(*t.shape)\n\nfig, (ax0, ax1) = plt.subplots(ncols=2, layout='constrained')\n\nyticks = np.arange(-50, 30, 10)\nyrange = (yticks[0], yticks[-1])\nxticks = np.arange(0, 550, 100)\n\nax0.psd(xn, NFFT=301, Fs=fs, window=mlab.window_none, pad_to=1024,\n        scale_by_freq=True)\nax0.set_title('Periodogram')\nax0.set_yticks(yticks)\nax0.set_xticks(xticks)\nax0.grid(True)\nax0.set_ylim(yrange)\n\nax1.psd(xn, NFFT=150, Fs=fs, window=mlab.window_none, pad_to=512, noverlap=75,\n        scale_by_freq=True)\nax1.set_title('Welch')\nax1.set_xticks(xticks)\nax1.set_yticks(yticks)\nax1.set_ylabel('')  # overwrite the y-label added by `psd`\nax1.grid(True)\nax1.set_ylim(yrange)\n\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This is a ported version of a MATLAB example from the signal\nprocessing toolbox that showed some difference at one time between\nMatplotlib's and MATLAB's scaling of the PSD.\n\nIt uses a complex signal so we can see that complex PSD's work properly.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "prng = np.random.RandomState(19680801)  # to ensure reproducibility\n\nfs = 1000\nt = np.linspace(0, 0.3, 301)\nA = np.array([2, 8]).reshape(-1, 1)\nf = np.array([150, 140]).reshape(-1, 1)\nxn = (A * np.exp(2j * np.pi * f * t)).sum(axis=0) + 5 * prng.randn(*t.shape)\n\nfig, (ax0, ax1) = plt.subplots(ncols=2, layout='constrained')\n\nyticks = np.arange(-50, 30, 10)\nyrange = (yticks[0], yticks[-1])\nxticks = np.arange(-500, 550, 200)\n\nax0.psd(xn, NFFT=301, Fs=fs, window=mlab.window_none, pad_to=1024,\n        scale_by_freq=True)\nax0.set_title('Periodogram')\nax0.set_yticks(yticks)\nax0.set_xticks(xticks)\nax0.grid(True)\nax0.set_ylim(yrange)\n\nax1.psd(xn, NFFT=150, Fs=fs, window=mlab.window_none, pad_to=512, noverlap=75,\n        scale_by_freq=True)\nax1.set_title('Welch')\nax1.set_xticks(xticks)\nax1.set_yticks(yticks)\nax1.set_ylabel('')  # overwrite the y-label added by `psd`\nax1.grid(True)\nax1.set_ylim(yrange)\n\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. tags::\n\n   domain: signal-processing\n   plot-type: line\n   level: intermediate\n\n"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "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.13.2"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}