.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "api/gallery/special_functions/convolution/plot_01_fast_convolution.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code or to run this example in your browser via JupyterLite. .. rst-class:: sphx-glr-example-title .. _sphx_glr_api_gallery_special_functions_convolution_plot_01_fast_convolution.py: Fast convolution: the convolution theorem and the FFT ===================================================== Stockham's 1966 fast convolution: multiplying zero-padded FFT spectra gives the same linear convolution as the direct O(nm) sum, at O((n+m) log(n+m)) cost. Without the zero padding the product gives the *circular* convolution instead. .. GENERATED FROM PYTHON SOURCE LINES 12-17 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np from mathematicskit.special_functions import circular_convolve, compare_convolution_methods, convolve_direct, convolve_fft .. GENERATED FROM PYTHON SOURCE LINES 18-21 Same result, two algorithms --------------------------- Smoothing a noisy step with a Gaussian kernel. .. GENERATED FROM PYTHON SOURCE LINES 21-37 .. code-block:: Python rng = np.random.default_rng(0) x = np.concatenate([np.zeros(200), np.ones(200)]) + 0.2 * rng.normal(size=400) kernel = np.exp(-0.5 * (np.arange(-30, 31) / 8.0) ** 2) kernel /= kernel.sum() direct = convolve_direct(x, kernel, mode="same") fast = convolve_fft(x, kernel, mode="same") print(f"max |direct - fft| = {np.max(np.abs(direct - fast)):.2e}") fig, ax = plt.subplots() ax.plot(x, lw=0.6, alpha=0.6, label="noisy step") ax.plot(fast, lw=2, label="Gaussian-smoothed (FFT convolution)") ax.legend() ax.set_title("Convolution with a Gaussian kernel") .. image-sg:: /api/gallery/special_functions/convolution/images/sphx_glr_plot_01_fast_convolution_001.png :alt: Convolution with a Gaussian kernel :srcset: /api/gallery/special_functions/convolution/images/sphx_glr_plot_01_fast_convolution_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none max |direct - fft| = 5.55e-16 Text(0.5, 1.0, 'Convolution with a Gaussian kernel') .. GENERATED FROM PYTHON SOURCE LINES 38-42 Linear vs. circular convolution ------------------------------- The DFT product without padding wraps the tail of the linear convolution back onto its start. .. GENERATED FROM PYTHON SOURCE LINES 42-48 .. code-block:: Python a = np.array([1.0, 2.0, 3.0, 4.0]) h = np.array([1.0, 1.0, 0.0, 0.0]) print("linear :", convolve_direct(a, h)) print("circular:", np.round(circular_convolve(a, h), 12)) .. rst-class:: sphx-glr-script-out .. code-block:: none linear : [1. 3. 5. 7. 4. 0. 0.] circular: [5. 3. 5. 7.] .. GENERATED FROM PYTHON SOURCE LINES 49-51 O(nm) vs. O((n+m) log(n+m)) --------------------------- .. GENERATED FROM PYTHON SOURCE LINES 51-67 .. code-block:: Python kernel_sizes = [8, 32, 128, 512, 2048, 8192] direct_times, fft_times = [], [] for m in kernel_sizes: result = compare_convolution_methods(20000, m, seed=0) direct_times.append(result.direct_time) fft_times.append(result.fft_time) print(f"m={m:>5}: direct={result.direct_time * 1000:8.3f} ms, fft={result.fft_time * 1000:7.3f} ms, max error={result.max_error:.1e}") fig, ax = plt.subplots() ax.loglog(kernel_sizes, direct_times, "o-", label="direct: O(nm)") ax.loglog(kernel_sizes, fft_times, "o-", label="FFT: O((n+m) log(n+m))") ax.set_xlabel("kernel length m (signal length n = 20000)") ax.set_ylabel("time (s)") ax.set_title("Direct vs. FFT convolution") ax.legend() .. image-sg:: /api/gallery/special_functions/convolution/images/sphx_glr_plot_01_fast_convolution_002.png :alt: Direct vs. FFT convolution :srcset: /api/gallery/special_functions/convolution/images/sphx_glr_plot_01_fast_convolution_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none m= 8: direct= 0.043 ms, fft= 0.774 ms, max error=6.2e-15 m= 32: direct= 0.218 ms, fft= 0.616 ms, max error=1.4e-14 m= 128: direct= 0.376 ms, fft= 0.584 ms, max error=3.2e-14 m= 512: direct= 1.582 ms, fft= 0.607 ms, max error=6.4e-14 m= 2048: direct= 6.314 ms, fft= 0.744 ms, max error=1.4e-13 m= 8192: direct= 29.173 ms, fft= 0.977 ms, max error=3.8e-13 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.289 seconds) .. _sphx_glr_download_api_gallery_special_functions_convolution_plot_01_fast_convolution.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: lite-badge .. image:: images/jupyterlite_badge_logo.svg :target: ../../../../lite/lab/index.html?path=api/gallery/special_functions/convolution/plot_01_fast_convolution.ipynb :alt: Launch JupyterLite :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_01_fast_convolution.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_01_fast_convolution.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_01_fast_convolution.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_