.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "intro/scipy/auto_examples/plot_fftpack.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. .. rst-class:: sphx-glr-example-title .. _sphx_glr_intro_scipy_auto_examples_plot_fftpack.py: ============================================= Plotting and manipulating FFTs for filtering ============================================= Plot the power of the FFT of a signal and inverse FFT back to reconstruct a signal. This example demonstrate :func:`scipy.fft.fft`, :func:`scipy.fft.fftfreq` and :func:`scipy.fft.ifft`. It implements a basic filter that is very suboptimal, and should not be used. .. GENERATED FROM PYTHON SOURCE LINES 15-20 .. code-block:: Python import numpy as np import scipy as sp import matplotlib.pyplot as plt .. GENERATED FROM PYTHON SOURCE LINES 21-23 Generate the signal ########################################################### .. GENERATED FROM PYTHON SOURCE LINES 23-36 .. code-block:: Python # Seed the random number generator rng = np.random.default_rng(27446968) time_step = 0.02 period = 5.0 time_vec = np.arange(0, 20, time_step) sig = np.sin(2 * np.pi / period * time_vec) + 0.5 * rng.normal(size=time_vec.size) plt.figure(figsize=(6, 5)) plt.plot(time_vec, sig, label="Original signal") .. image-sg:: /intro/scipy/auto_examples/images/sphx_glr_plot_fftpack_001.png :alt: plot fftpack :srcset: /intro/scipy/auto_examples/images/sphx_glr_plot_fftpack_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [] .. GENERATED FROM PYTHON SOURCE LINES 37-39 Compute and plot the power ########################################################### .. GENERATED FROM PYTHON SOURCE LINES 39-73 .. code-block:: Python # The FFT of the signal sig_fft = sp.fft.fft(sig) # And the power (sig_fft is of complex dtype) power = np.abs(sig_fft) ** 2 # The corresponding frequencies sample_freq = sp.fft.fftfreq(sig.size, d=time_step) # Plot the FFT power plt.figure(figsize=(6, 5)) plt.plot(sample_freq, power) plt.xlabel("Frequency [Hz]") plt.ylabel("plower") # Find the peak frequency: we can focus on only the positive frequencies pos_mask = np.where(sample_freq > 0) freqs = sample_freq[pos_mask] peak_freq = freqs[power[pos_mask].argmax()] # Check that it does indeed correspond to the frequency that we generate # the signal with np.allclose(peak_freq, 1.0 / period) # An inner plot to show the peak frequency axes = plt.axes((0.55, 0.3, 0.3, 0.5)) plt.title("Peak frequency") plt.plot(freqs[:8], power[pos_mask][:8]) plt.setp(axes, yticks=[]) # scipy.signal.find_peaks_cwt can also be used for more advanced # peak detection .. image-sg:: /intro/scipy/auto_examples/images/sphx_glr_plot_fftpack_002.png :alt: Peak frequency :srcset: /intro/scipy/auto_examples/images/sphx_glr_plot_fftpack_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [] .. GENERATED FROM PYTHON SOURCE LINES 74-79 Remove all the high frequencies ########################################################### We now remove all the high frequencies and transform back from frequencies to signal. .. GENERATED FROM PYTHON SOURCE LINES 79-92 .. code-block:: Python high_freq_fft = sig_fft.copy() high_freq_fft[np.abs(sample_freq) > peak_freq] = 0 filtered_sig = sp.fft.ifft(high_freq_fft) plt.figure(figsize=(6, 5)) plt.plot(time_vec, sig, label="Original signal") plt.plot(time_vec, filtered_sig, linewidth=3, label="Filtered signal") plt.xlabel("Time [s]") plt.ylabel("Amplitude") plt.legend(loc="best") .. image-sg:: /intro/scipy/auto_examples/images/sphx_glr_plot_fftpack_003.png :alt: plot fftpack :srcset: /intro/scipy/auto_examples/images/sphx_glr_plot_fftpack_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none /home/docs/checkouts/readthedocs.org/user_builds/scientific-python-lectures-ja/envs/latest/lib/python3.12/site-packages/matplotlib/cbook.py:1719: ComplexWarning: Casting complex values to real discards the imaginary part return math.isfinite(val) /home/docs/checkouts/readthedocs.org/user_builds/scientific-python-lectures-ja/envs/latest/lib/python3.12/site-packages/matplotlib/cbook.py:1355: ComplexWarning: Casting complex values to real discards the imaginary part return np.asarray(x, float) .. GENERATED FROM PYTHON SOURCE LINES 93-97 **Note** This is actually a bad way of creating a filter: such brutal cut-off in frequency space does not control distortion on the signal. Filters should be created using the SciPy filter design code .. GENERATED FROM PYTHON SOURCE LINES 98-99 .. code-block:: Python plt.show() .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.148 seconds) .. _sphx_glr_download_intro_scipy_auto_examples_plot_fftpack.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_fftpack.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_fftpack.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_fftpack.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_