.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/contrib/sinetransform.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_examples_contrib_sinetransform.py: Sine Transform ============== Contributed by `@ShazAlvi `_. This is a simple test program to illustrate how the sine (or cosine as it works basically the same way ) Fourier transform works using `FFTLog`. The test provides as input as sine function and performs the sine Fourier transform. The input function is then recovered by performing an inverse Fourier transform. The inverse is performed using the following integral, .. math:: :label: sinetest F(t) = \sqrt{\frac{\pi}{2}}\int^\infty_0 A(f)\ \sin(ft) \ \text{d}f \ . .. GENERATED FROM PYTHON SOURCE LINES 19-24 .. code-block:: Python import pyfftlog import numpy as np import scipy as sp import matplotlib.pyplot as plt .. GENERATED FROM PYTHON SOURCE LINES 25-28 Define the parameters you wish to use ------------------------------------- The presets are the *Reasonable choices of parameters* from `fftlogtest.f`. .. GENERATED FROM PYTHON SOURCE LINES 28-64 .. code-block:: Python # Range of periodic interval logtmin = -3 # logtmax = 0.798 #2pi # 5(2pi) #Longer range in r gives you a better reconstruction. 10\pi will give # you a better reconstruction than 2\pi. logtmax = 1.497 # Number of points (Max 4096) # 1000 points give you a fairly smooth distribution of af in frequency, f. # However you can get a good, working fit for 300 points as well. n = 1000 # Order mu of Bessel function mu = 0.5 # Choose -0.5 for cosine fourier transform # Bias exponent: q = 0 is unbiased # The unbiased transforms give better results as far as I checked. q = 0 # Sensible approximate choice of f_c t_c # The output and the reconstruction is sensitive to the choice of this value # This value is found by trial and error. In this example, the input function # is a simple sine function which is not smooth in frequency space (as it # only has one frequency) because of this reason a better value of this # quantity is not found by the function fhti. For functions smooth # in both time and frequency domain, the fhti should return the best # value of the f_c t_c. ft = 0.016 # Tell fhti to change ft to low-ringing value # WARNING: kropt = 3 will fail, as interaction is not supported ftopt = 1 # Forward transform (changed from dir to tdir, as dir is a python fct) tdir = 1 .. GENERATED FROM PYTHON SOURCE LINES 65-67 Computation related to the logarithmic spacing ---------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 67-82 .. code-block:: Python # Central point log10(t_c) of periodic interval logtc = (logtmin + logtmax)/2 print(f"Central point of periodic interval at log10(t_c) = {logtc}") # Central index (1/2 integral if n is even) nc = (n + 1)/2.0 # Log-spacing of points dlogt = (logtmax - logtmin)/n dlnr = dlogt*np.log(10.0) .. rst-class:: sphx-glr-script-out .. code-block:: none Central point of periodic interval at log10(t_c) = -0.7515 .. GENERATED FROM PYTHON SOURCE LINES 83-85 Compute input function: :math:`\sin(t)` --------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 85-89 .. code-block:: Python t = 10**(logtc + (np.arange(1, n+1) - nc)*dlogt) a_t = np.sin(t) .. GENERATED FROM PYTHON SOURCE LINES 90-92 Initialize FFTLog transform - note `fhti` resets `ft` ----------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 92-95 .. code-block:: Python ft, xsave = pyfftlog.fhti(n, mu, dlnr, q, ft, ftopt) .. GENERATED FROM PYTHON SOURCE LINES 96-98 Call `pyfftlog.fhtl` -------------------- .. GENERATED FROM PYTHON SOURCE LINES 98-133 .. code-block:: Python logfc = np.log10(ft) - logtc # Fourier sine Transform a_f = pyfftlog.fftl(a_t.copy(), xsave, np.sqrt(2/np.pi), tdir) # Notice that np.sqrt(2/np.pi) is the normalization factor for the transform # Reconstruct the input function by taking the inverse fourier transform as # given in the description f = 10**(logfc + (np.arange(1, n+1) - nc)*dlogt) # Array to store the reconstructed function for each value of t Recon_Fun = np.zeros((len(t))) for i in range(len(t)): Recon_Fun[i] = (np.sqrt(2/np.pi)**-1) * \ sp.integrate.trapezoid(f, a_f*np.sin(t[i]*f)) # Plotting the input function and the reconstructed input function and also # the distribution of the a(f) vs f. plt.figure() ax1 = plt.subplot(121) plt.title(r'Frequency domain') plt.xlabel('f') plt.ylabel(r'$a_f(f)$') plt.semilogx(f, a_f, 'k') ax2 = plt.subplot(122) plt.title('Time domain') plt.xlabel("t") plt.ylabel("sin(t)") plt.semilogx(t, a_t, lw=2, label=r'$\sin(t)$') plt.semilogx(t, -Recon_Fun, '--', label='Reconstructed') plt.legend() plt.tight_layout() plt.show() .. image-sg:: /examples/contrib/images/sphx_glr_sinetransform_001.png :alt: Frequency domain, Time domain :srcset: /examples/contrib/images/sphx_glr_sinetransform_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.518 seconds) .. _sphx_glr_download_examples_contrib_sinetransform.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: sinetransform.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: sinetransform.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: sinetransform.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_