.. note:: :class: sphx-glr-download-link-note Click :ref:`here ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_gallery_lines_bars_and_markers_psd_demo.py: ======== Psd Demo ======== Plotting Power Spectral Density (PSD) in Matplotlib. The PSD is a common plot in the field of signal processing. NumPy has many useful libraries for computing a PSD. Below we demo a few examples of how this can be accomplished and visualized with Matplotlib. .. code-block:: python import matplotlib.pyplot as plt import numpy as np import matplotlib.mlab as mlab # Fixing random state for reproducibility np.random.seed(19680801) dt = 0.01 t = np.arange(0, 10, dt) nse = np.random.randn(len(t)) r = np.exp(-t / 0.05) cnse = np.convolve(nse, r) * dt cnse = cnse[:len(t)] s = 0.1 * np.sin(2 * np.pi * t) + cnse plt.subplot(211) plt.plot(t, s) plt.subplot(212) plt.psd(s, 512, 1 / dt) plt.show() .. image:: /gallery/lines_bars_and_markers/images/sphx_glr_psd_demo_001.png :class: sphx-glr-single-img Compare this with the equivalent Matlab code to accomplish the same thing:: dt = 0.01; t = [0:dt:10]; nse = randn(size(t)); r = exp(-t/0.05); cnse = conv(nse, r)*dt; cnse = cnse(1:length(t)); s = 0.1*sin(2*pi*t) + cnse; subplot(211) plot(t,s) subplot(212) psd(s, 512, 1/dt) Below we'll show a slightly more complex example that demonstrates how padding affects the resulting PSD. .. code-block:: python dt = np.pi / 100. fs = 1. / dt t = np.arange(0, 8, dt) y = 10. * np.sin(2 * np.pi * 4 * t) + 5. * np.sin(2 * np.pi * 4.25 * t) y = y + np.random.randn(*t.shape) # Plot the raw time series fig = plt.figure() fig.subplots_adjust(hspace=0.45, wspace=0.3) ax = fig.add_subplot(2, 1, 1) ax.plot(t, y) # Plot the PSD with different amounts of zero padding. This uses the entire # time series at once ax2 = fig.add_subplot(2, 3, 4) ax2.psd(y, NFFT=len(t), pad_to=len(t), Fs=fs) ax2.psd(y, NFFT=len(t), pad_to=len(t) * 2, Fs=fs) ax2.psd(y, NFFT=len(t), pad_to=len(t) * 4, Fs=fs) plt.title('zero padding') # Plot the PSD with different block sizes, Zero pad to the length of the # original data sequence. ax3 = fig.add_subplot(2, 3, 5, sharex=ax2, sharey=ax2) ax3.psd(y, NFFT=len(t), pad_to=len(t), Fs=fs) ax3.psd(y, NFFT=len(t) // 2, pad_to=len(t), Fs=fs) ax3.psd(y, NFFT=len(t) // 4, pad_to=len(t), Fs=fs) ax3.set_ylabel('') plt.title('block size') # Plot the PSD with different amounts of overlap between blocks ax4 = fig.add_subplot(2, 3, 6, sharex=ax2, sharey=ax2) ax4.psd(y, NFFT=len(t) // 2, pad_to=len(t), noverlap=0, Fs=fs) ax4.psd(y, NFFT=len(t) // 2, pad_to=len(t), noverlap=int(0.05 * len(t) / 2.), Fs=fs) ax4.psd(y, NFFT=len(t) // 2, pad_to=len(t), noverlap=int(0.2 * len(t) / 2.), Fs=fs) ax4.set_ylabel('') plt.title('overlap') plt.show() .. image:: /gallery/lines_bars_and_markers/images/sphx_glr_psd_demo_002.png :class: sphx-glr-single-img This is a ported version of a MATLAB example from the signal processing toolbox that showed some difference at one time between Matplotlib's and MATLAB's scaling of the PSD. .. code-block:: python fs = 1000 t = np.linspace(0, 0.3, 301) A = np.array([2, 8]).reshape(-1, 1) f = np.array([150, 140]).reshape(-1, 1) xn = (A * np.sin(2 * np.pi * f * t)).sum(axis=0) xn += 5 * np.random.randn(*t.shape) fig, (ax0, ax1) = plt.subplots(ncols=2) fig.subplots_adjust(hspace=0.45, wspace=0.3) yticks = np.arange(-50, 30, 10) yrange = (yticks[0], yticks[-1]) xticks = np.arange(0, 550, 100) ax0.psd(xn, NFFT=301, Fs=fs, window=mlab.window_none, pad_to=1024, scale_by_freq=True) ax0.set_title('Periodogram') ax0.set_yticks(yticks) ax0.set_xticks(xticks) ax0.grid(True) ax0.set_ylim(yrange) ax1.psd(xn, NFFT=150, Fs=fs, window=mlab.window_none, pad_to=512, noverlap=75, scale_by_freq=True) ax1.set_title('Welch') ax1.set_xticks(xticks) ax1.set_yticks(yticks) ax1.set_ylabel('') # overwrite the y-label added by `psd` ax1.grid(True) ax1.set_ylim(yrange) plt.show() .. image:: /gallery/lines_bars_and_markers/images/sphx_glr_psd_demo_003.png :class: sphx-glr-single-img This is a ported version of a MATLAB example from the signal processing toolbox that showed some difference at one time between Matplotlib's and MATLAB's scaling of the PSD. It uses a complex signal so we can see that complex PSD's work properly. .. code-block:: python prng = np.random.RandomState(19680801) # to ensure reproducibility fs = 1000 t = np.linspace(0, 0.3, 301) A = np.array([2, 8]).reshape(-1, 1) f = np.array([150, 140]).reshape(-1, 1) xn = (A * np.exp(2j * np.pi * f * t)).sum(axis=0) + 5 * prng.randn(*t.shape) fig, (ax0, ax1) = plt.subplots(ncols=2) fig.subplots_adjust(hspace=0.45, wspace=0.3) yticks = np.arange(-50, 30, 10) yrange = (yticks[0], yticks[-1]) xticks = np.arange(-500, 550, 200) ax0.psd(xn, NFFT=301, Fs=fs, window=mlab.window_none, pad_to=1024, scale_by_freq=True) ax0.set_title('Periodogram') ax0.set_yticks(yticks) ax0.set_xticks(xticks) ax0.grid(True) ax0.set_ylim(yrange) ax1.psd(xn, NFFT=150, Fs=fs, window=mlab.window_none, pad_to=512, noverlap=75, scale_by_freq=True) ax1.set_title('Welch') ax1.set_xticks(xticks) ax1.set_yticks(yticks) ax1.set_ylabel('') # overwrite the y-label added by `psd` ax1.grid(True) ax1.set_ylim(yrange) plt.show() .. image:: /gallery/lines_bars_and_markers/images/sphx_glr_psd_demo_004.png :class: sphx-glr-single-img .. _sphx_glr_download_gallery_lines_bars_and_markers_psd_demo.py: .. only :: html .. container:: sphx-glr-footer :class: sphx-glr-footer-example .. container:: sphx-glr-download :download:`Download Python source code: psd_demo.py ` .. container:: sphx-glr-download :download:`Download Jupyter notebook: psd_demo.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature Keywords: matplotlib code example, codex, python plot, pyplot `Gallery generated by Sphinx-Gallery `_