.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "api/gallery/signal_processing/fourier/plot_02_cooley_tukey_fft.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_signal_processing_fourier_plot_02_cooley_tukey_fft.py: Cooley and Tukey's fast Fourier transform (1965) ================================================ A direct DFT of :math:`N` samples needs :math:`N^2` complex multiplications. James Cooley and John Tukey (1965) split it into the DFTs of the even- and odd-numbered samples, combined by "butterflies" :math:`E_k \pm W^k O_k`, and repeated the split down to single samples (:func:`~engineeringkit.signal_processing.radix2_fft`). The cost falls to :math:`(N/2)\log_2 N`. Gauss had used the same idea in 1805 to interpolate asteroid orbits. In 1965 it made spectral analysis on a computer routine, and it is what :func:`~engineeringkit.signal_processing.dft` (NumPy's FFT) does. .. GENERATED FROM PYTHON SOURCE LINES 17-22 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np from engineeringkit.signal_processing import convolution, dft, fft_frequencies, radix2_fft .. GENERATED FROM PYTHON SOURCE LINES 23-28 The butterflies give the DFT ---------------------------- A direct DFT is a matrix product with :math:`N^2` entries :math:`e^{-j2\pi kn/N}`. .. GENERATED FROM PYTHON SOURCE LINES 28-37 .. code-block:: Python rng = np.random.default_rng(1965) x = rng.normal(size=1024) + 1j * rng.normal(size=1024) n = np.arange(x.size) direct = np.exp(-2j * np.pi * np.outer(n, n) / x.size) @ x fast, count = radix2_fft(x, return_multiplications=True) assert np.allclose(fast, direct) and np.allclose(fast, dft(x)) print(f"N = {x.size}: direct {x.size**2} multiplications, radix-2 {count}, speed-up {x.size**2 / count:.0f}x") .. rst-class:: sphx-glr-script-out .. code-block:: none N = 1024: direct 1048576 multiplications, radix-2 5120, speed-up 205x .. GENERATED FROM PYTHON SOURCE LINES 38-40 Operation counts ---------------- .. GENERATED FROM PYTHON SOURCE LINES 40-51 .. code-block:: Python sizes = 2 ** np.arange(3, 17) counts = [radix2_fft(np.zeros(size), return_multiplications=True)[1] for size in sizes[:10]] assert counts == [size // 2 * int(np.log2(size)) for size in sizes[:10]] fig, axes = plt.subplots(1, 2, figsize=(10, 3.8)) axes[0].loglog(sizes, sizes**2, label="direct DFT, $N^2$") axes[0].loglog(sizes, sizes / 2 * np.log2(sizes), label=r"FFT, $(N/2)\log_2 N$") axes[0].loglog(sizes[:10], counts, "o", label="counted in radix2_fft") axes[0].set(xlabel="N", ylabel="complex multiplications") axes[0].legend() .. image-sg:: /api/gallery/signal_processing/fourier/images/sphx_glr_plot_02_cooley_tukey_fft_001.png :alt: plot 02 cooley tukey fft :srcset: /api/gallery/signal_processing/fourier/images/sphx_glr_plot_02_cooley_tukey_fft_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 52-58 Fast convolution ---------------- The FFT turns convolution into multiplication: zero-pad both sequences to the output length, multiply their transforms, and invert. For a long filter this beats the direct sum. .. GENERATED FROM PYTHON SOURCE LINES 58-76 .. code-block:: Python signal_in = rng.normal(size=4000) h = np.exp(-np.arange(300) / 50.0) length = int(2 ** np.ceil(np.log2(signal_in.size + h.size - 1))) fast_conv = np.fft.ifft(dft(np.pad(signal_in, (0, length - signal_in.size))) * dft(np.pad(h, (0, length - h.size)))).real fast_conv = fast_conv[: signal_in.size + h.size - 1] assert np.allclose(fast_conv, convolution(signal_in, h)) print(f"FFT convolution of {signal_in.size} by {h.size} samples matches the direct sum") fs = 1000.0 tone = np.sin(2 * np.pi * 125 * np.arange(64) / fs) spectrum = radix2_fft(tone) axes[1].stem(np.fft.fftshift(fft_frequencies(64, fs)), np.fft.fftshift(abs(spectrum)), basefmt=" ") axes[1].set(xlabel="frequency [Hz]", ylabel="|X|", title="radix2_fft of a 125 Hz tone") assert abs(fft_frequencies(64, fs)[np.argmax(abs(spectrum[:32]))] - 125) < 1e-9 fig.tight_layout() plt.show() .. rst-class:: sphx-glr-script-out .. code-block:: none FFT convolution of 4000 by 300 samples matches the direct sum .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.111 seconds) .. _sphx_glr_download_api_gallery_signal_processing_fourier_plot_02_cooley_tukey_fft.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/signal_processing/fourier/plot_02_cooley_tukey_fft.ipynb :alt: Launch JupyterLite :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_02_cooley_tukey_fft.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_02_cooley_tukey_fft.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_02_cooley_tukey_fft.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_