.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "_build/auto_examples/Numerical_tricks/FFT_sym.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__build_auto_examples_Numerical_tricks_FFT_sym.py: FFT and symmetries ================== A few numerical investigation of how symmetries of a signal determine the shape of its discrete Fourier transform, and how to take advantage of this fact to optimize performance. .. GENERATED FROM PYTHON SOURCE LINES 10-18 The spectral solvers in choreo make extensive use of the discrete Fourier transform to go back end forth between time representation (positions :math:`x_j = x(j/n)`) and frequency representation (Fourier coefficients :math:`c_k`) representation of orbits. Explicit formulae read: .. math:: c_k = \sum_{j=0}^{n-1} x_j \exp\left(-2i\pi k\tfrac{j}{n} \right) \\ x_k = \frac{1}{n}\sum_{j=0}^{n-1} c_j \exp\left(2i\pi k\tfrac{j}{n} \right) \\ The magic of FFT algorithms enables computation of **all** the positions from **all** the coefficients (and vice-versa) in :math:`\mathcal{O}(n\log(n))` time instead of the naive :math:`\mathcal{O}(n^2)`. .. GENERATED FROM PYTHON SOURCE LINES 21-22 Note that the total number of nodes needs to be divisible by 4 to properly account for symmetries .. GENERATED FROM PYTHON SOURCE LINES 22-32 .. code-block:: Python import numpy as np import scipy import matplotlib.pyplot as plt import sys if ("--no-show" in sys.argv): plt.show = (lambda : None) n_base = 8 # Base number of nodes. This could be any **even** integer > 0 .. GENERATED FROM PYTHON SOURCE LINES 33-37 Real-valued signals and RFFT **************************** Signals we use in choreo are coordinates of the orbits, which are real-valued. .. GENERATED FROM PYTHON SOURCE LINES 37-48 .. code-block:: Python n = n_base # Total number of nodes. t = np.array(range(n))/n x = np.array(range(n))+1 err = np.linalg.norm(x.imag) print(err) plt.plot(t,x) plt.show() .. image-sg:: /_build/auto_examples/Numerical_tricks/images/sphx_glr_FFT_sym_001.png :alt: FFT sym :srcset: /_build/auto_examples/Numerical_tricks/images/sphx_glr_FFT_sym_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none 0.0 .. GENERATED FROM PYTHON SOURCE LINES 49-50 Therefore the FFT has Hermitian symmetry, meaning that the Fourier coefficients of negative index are the complex conjugates of their positive counterparts .. GENERATED FROM PYTHON SOURCE LINES 50-57 .. code-block:: Python fft_c = scipy.fft.fft(x) err = ( abs(fft_c[0] - fft_c[0].conjugate()) + np.linalg.norm(fft_c[1:n//2] - fft_c[n//2+1:][::-1].conjugate()) ) print(err) .. rst-class:: sphx-glr-script-out .. code-block:: none 0.0 .. GENERATED FROM PYTHON SOURCE LINES 58-59 This symmetry can be leveraged to only compute coefficients of positive index using the RFFT .. GENERATED FROM PYTHON SOURCE LINES 59-64 .. code-block:: Python rfft_c = scipy.fft.rfft(x) err = np.linalg.norm(fft_c[:n//2+1] - rfft_c) print(err) .. rst-class:: sphx-glr-script-out .. code-block:: none 0.0 .. GENERATED FROM PYTHON SOURCE LINES 65-66 The positions are retrieved using the IRFFT. .. GENERATED FROM PYTHON SOURCE LINES 66-71 .. code-block:: Python irfft_c = scipy.fft.irfft(rfft_c) err = np.linalg.norm(irfft_c - x) print(err) .. rst-class:: sphx-glr-script-out .. code-block:: none 0.0 .. GENERATED FROM PYTHON SOURCE LINES 72-76 Real-valued even signals and DCT I ********************************** Suppose a coordinate of the orbit is constrained to be even. .. GENERATED FROM PYTHON SOURCE LINES 76-97 .. code-block:: Python n = 2*n_base # Total number of nodes. Necessarily even here. t = np.array(range(n))/n y = np.array(range(n//2+1))+1 x = np.zeros(n) for i in range(n//2+1): x[i] = y[i] for i in range(n//2+1,n): x[i] = y[n - i] z = np.zeros(n) for i in range(n): z[i] = x[i] - x[(n-i) % n] err = np.linalg.norm(z) print(err) plt.plot(t,x) plt.show() .. image-sg:: /_build/auto_examples/Numerical_tricks/images/sphx_glr_FFT_sym_002.png :alt: FFT sym :srcset: /_build/auto_examples/Numerical_tricks/images/sphx_glr_FFT_sym_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none 0.0 .. GENERATED FROM PYTHON SOURCE LINES 98-100 In this case, the Fourier transform is purely real, i.e. its imaginary part is zero. .. GENERATED FROM PYTHON SOURCE LINES 100-104 .. code-block:: Python rfft_c = scipy.fft.rfft(x) err = np.linalg.norm(rfft_c.imag) print(err) .. rst-class:: sphx-glr-script-out .. code-block:: none 8.881784197001252e-16 .. GENERATED FROM PYTHON SOURCE LINES 105-106 This symmetry can be leveraged to only compute the real part of the coefficients using the DCT I .. GENERATED FROM PYTHON SOURCE LINES 106-111 .. code-block:: Python dct_I_c = scipy.fft.dct(x[0:n//2+1],1) err = np.linalg.norm(rfft_c.real - dct_I_c) print(err) .. rst-class:: sphx-glr-script-out .. code-block:: none 0.0 .. GENERATED FROM PYTHON SOURCE LINES 112-113 Half the positions are retrieved using the IDCT I. .. GENERATED FROM PYTHON SOURCE LINES 113-118 .. code-block:: Python idct_I_x = scipy.fft.idct(dct_I_c,1) err = np.linalg.norm(x[0:n//2+1] - idct_I_x) print(err) .. rst-class:: sphx-glr-script-out .. code-block:: none 4.440892098500626e-16 .. GENERATED FROM PYTHON SOURCE LINES 119-123 Real-valued odd signals and DST I ********************************* Suppose a coordinate of the orbit is constrained to be odd. .. GENERATED FROM PYTHON SOURCE LINES 123-145 .. code-block:: Python n = 2*n_base # Total number of nodes. Necessarily even here. t = np.array(range(n))/n y = np.array(range(n//2-1))+1 x = np.zeros(n) for i in range(1,n//2): x[i] = y[i-1] for i in range(n//2+1,n): x[i] = - y[n - i - 1] z = np.zeros(n) for i in range(n): z[i] = x[i] + x[(n-i) % n] err = np.linalg.norm(z) print(err) plt.plot(t,x) plt.show() .. image-sg:: /_build/auto_examples/Numerical_tricks/images/sphx_glr_FFT_sym_003.png :alt: FFT sym :srcset: /_build/auto_examples/Numerical_tricks/images/sphx_glr_FFT_sym_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none 0.0 .. GENERATED FROM PYTHON SOURCE LINES 146-147 In this case, the Fourier transform is purely imaginary, i.e. its real part is zero. .. GENERATED FROM PYTHON SOURCE LINES 147-152 .. code-block:: Python rfft_c = scipy.fft.rfft(x) err = np.linalg.norm(rfft_c.real) print(err) .. rst-class:: sphx-glr-script-out .. code-block:: none 1.7763568394002505e-15 .. GENERATED FROM PYTHON SOURCE LINES 153-154 This symmetry can be leveraged to only compute the imaginary part of the coefficients using the DST I .. GENERATED FROM PYTHON SOURCE LINES 154-159 .. code-block:: Python dst_I_c = scipy.fft.dst(x[1:n//2],1) err = np.linalg.norm(rfft_c.imag[1:n//2] - (- dst_I_c)) print(err) .. rst-class:: sphx-glr-script-out .. code-block:: none 0.0 .. GENERATED FROM PYTHON SOURCE LINES 160-161 Half the positions are retrieved using the IDST I. .. GENERATED FROM PYTHON SOURCE LINES 161-166 .. code-block:: Python idst_I_x = scipy.fft.idst(dst_I_c,1) err = np.linalg.norm(x[1:n//2] - idst_I_x) print(err) .. rst-class:: sphx-glr-script-out .. code-block:: none 4.440892098500626e-16 .. GENERATED FROM PYTHON SOURCE LINES 167-171 Real-valued even-odd signals and DCT III **************************************** Suppose a coordinate of the orbit is constrained to be even arround the origin :math:`t=0` and odd arround a quarter period :math:`t=T/4`. .. GENERATED FROM PYTHON SOURCE LINES 171-205 .. code-block:: Python n = 4*n_base # Total number of nodes. Necessarily divisible by 4 here. t = np.array(range(n))/n y = np.array(range(n//4))+1 x = np.zeros(n) for i in range(n//4): x[i] = y[i] for i in range(n//4+1,n//2): x[i] = - y[n//2 - i] for i in range(n//2,n//2 + n//4): x[i] = - y[i-n//2] for i in range(n//2 + n//4 + 1,n): x[i] = y[n - i] z = np.zeros(n) for i in range(n): z[i] = x[i] - x[(n-i) % n] err = np.linalg.norm(z) print(err) for i in range(n): z[i] = x[i] + x[(n+n//2 - i) % n] err = np.linalg.norm(z) print(err) plt.plot(t,x) plt.show() .. image-sg:: /_build/auto_examples/Numerical_tricks/images/sphx_glr_FFT_sym_004.png :alt: FFT sym :srcset: /_build/auto_examples/Numerical_tricks/images/sphx_glr_FFT_sym_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none 0.0 0.0 .. GENERATED FROM PYTHON SOURCE LINES 206-207 In this case, the Fourier transform is purely real and its odd-numbered coefficients are zero. .. GENERATED FROM PYTHON SOURCE LINES 207-215 .. code-block:: Python rfft_c = scipy.fft.rfft(x) err = np.linalg.norm(rfft_c.imag) print(err) err = np.linalg.norm(rfft_c[::2]) print(err) .. rst-class:: sphx-glr-script-out .. code-block:: none 8.51910812704787e-15 0.0 .. GENERATED FROM PYTHON SOURCE LINES 216-217 This symmetry can be leveraged to only compute the non-zero coefficients using the DCT III .. GENERATED FROM PYTHON SOURCE LINES 217-222 .. code-block:: Python dct_III_c = scipy.fft.dct(x[0:n//4],3) err = np.linalg.norm(rfft_c.real[1::2] - 2*dct_III_c) print(err) .. rst-class:: sphx-glr-script-out .. code-block:: none 6.646518689688359e-15 .. GENERATED FROM PYTHON SOURCE LINES 223-224 A quarter of the positions are retrieved using the IDCT III. .. GENERATED FROM PYTHON SOURCE LINES 224-229 .. code-block:: Python idct_III_c = scipy.fft.idct(dct_III_c,3) err = np.linalg.norm(x[0:n//4] - idct_III_c) print(err) .. rst-class:: sphx-glr-script-out .. code-block:: none 1.831026719408895e-15 .. GENERATED FROM PYTHON SOURCE LINES 230-234 Real-valued odd-even signals and DST III **************************************** Suppose a coordinate of the orbit is constrained to be odd arround the origin :math:`t=0` and even arround a quarter period :math:`t=T/4`. .. GENERATED FROM PYTHON SOURCE LINES 234-269 .. code-block:: Python n = 4*n_base # Total number of nodes. Necessarily divisible by 4 here. t = np.array(range(n))/n y = np.array(range(n//4))+1 x = np.zeros(n) for i in range(1,n//4 + 1): x[i] = y[i-1] for i in range(n//4+1,n//2): x[i] = y[n//2-1 - i] for i in range(n//2+1,n//2 + n//4+1): x[i] = - y[i-(n//2+1)] for i in range(n//2 + n//4+1,n): x[i] = -y[n-1 - i] z = np.zeros(n) for i in range(n): z[i] = x[i] + x[(n-i) % n] err = np.linalg.norm(z) print(err) for i in range(n): z[i] = x[i] - x[(n+n//2 - i) % n] err = np.linalg.norm(z) print(err) plt.plot(t,x) plt.show() .. image-sg:: /_build/auto_examples/Numerical_tricks/images/sphx_glr_FFT_sym_005.png :alt: FFT sym :srcset: /_build/auto_examples/Numerical_tricks/images/sphx_glr_FFT_sym_005.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none 0.0 0.0 .. GENERATED FROM PYTHON SOURCE LINES 270-271 In this case, the Fourier transform is purely imaginary and its even-numbered coefficients are zero. .. GENERATED FROM PYTHON SOURCE LINES 271-278 .. code-block:: Python rfft_c = scipy.fft.rfft(x) err = np.linalg.norm(rfft_c.real) print(err) err = np.linalg.norm(rfft_c[::2]) print(err) .. rst-class:: sphx-glr-script-out .. code-block:: none 5.511875769711332e-15 0.0 .. GENERATED FROM PYTHON SOURCE LINES 279-280 This symmetry can be leveraged to only compute the non-zero coefficients using the DST III .. GENERATED FROM PYTHON SOURCE LINES 280-285 .. code-block:: Python dst_III_c = scipy.fft.dst(x[1:n//4+1],3) err = np.linalg.norm(rfft_c.imag[1::2] - (-2*dst_III_c)) print(err) .. rst-class:: sphx-glr-script-out .. code-block:: none 9.532413939823706e-15 .. GENERATED FROM PYTHON SOURCE LINES 286-287 A quarter of the positions are retrieved using the IDCT III. .. GENERATED FROM PYTHON SOURCE LINES 287-292 .. code-block:: Python idst_III_c = scipy.fft.idst(dst_III_c,3) err = np.linalg.norm(x[1:n//4+1] - idst_III_c) print(err) .. rst-class:: sphx-glr-script-out .. code-block:: none 2.3914935841127266e-15 .. GENERATED FROM PYTHON SOURCE LINES 293-297 Real-valued signals with a subperiod ************************************ Suppose a coordinate of the orbit is constrained to have a period that is a multiple of the base period .. GENERATED FROM PYTHON SOURCE LINES 297-319 .. code-block:: Python sub_period = 3 n = n_base * sub_period # Total number of nodes. Necessarily a multiple of the sub period. t = np.array(range(n))/n y = np.array(range(n_base)) x = np.zeros(n) for i in range(sub_period): for j in range(n_base): x[j+i*n_base] = y[j] z = np.zeros(n) for i_sub in range(sub_period-1): for i in range(n): z[i] = x[i] - x[(i+(i_sub+1)*n_base) % n] err = np.linalg.norm(z) print(err) plt.plot(t,x) plt.show() .. image-sg:: /_build/auto_examples/Numerical_tricks/images/sphx_glr_FFT_sym_006.png :alt: FFT sym :srcset: /_build/auto_examples/Numerical_tricks/images/sphx_glr_FFT_sym_006.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none 0.0 0.0 .. GENERATED FROM PYTHON SOURCE LINES 320-321 In this case, the coefficients of the Fourier transform that are not indexed by a multiple of the subperiod are zero. .. GENERATED FROM PYTHON SOURCE LINES 321-329 .. code-block:: Python rfft_c = scipy.fft.rfft(x) err = 0 for i in range(1,sub_period): err += np.linalg.norm(rfft_c[i::sub_period]) print(err) .. rst-class:: sphx-glr-script-out .. code-block:: none 0.0 .. GENERATED FROM PYTHON SOURCE LINES 330-331 This symmetry can be leveraged to only compute the non-zero coefficients using the RFFT on a sub-period only .. GENERATED FROM PYTHON SOURCE LINES 331-336 .. code-block:: Python rfftsub_c = scipy.fft.rfft(x[:n_base]) err = np.linalg.norm(rfft_c[::sub_period] - sub_period*rfftsub_c) print(err) .. rst-class:: sphx-glr-script-out .. code-block:: none 6.4047456679787536e-15 .. GENERATED FROM PYTHON SOURCE LINES 337-338 The positions on a sub-period are retrieved using the RFFT. .. GENERATED FROM PYTHON SOURCE LINES 338-342 .. code-block:: Python irfftsub_c = scipy.fft.irfft(rfftsub_c) err = np.linalg.norm(x[:n_base] - irfftsub_c) print(err) .. rst-class:: sphx-glr-script-out .. code-block:: none 4.440892098500626e-16 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.407 seconds) .. _sphx_glr_download__build_auto_examples_Numerical_tricks_FFT_sym.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: FFT_sym.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: FFT_sym.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_