Add documentation. Fix bugs.

docs/index.rst

@@ -279,6 +279,7 @@ Using **synphot**
    synphot/observation
    synphot/formulae
    synphot/units
+   synphot/filter_par
    synphot/tutorials
 
 .. _synphot_history:

docs/synphot/filter_par.rst

@@ -0,0 +1,120 @@
+.. _synphot_par_filters:
+
+Parameterized Filters
+=====================
+
+.. note::
+
+    The algorithm for parameterized filters here was originally developed by
+    Brett Morris for the `tynt <https://github.com/bmorris3/tynt/>`_ package.
+
+Some filters could be approximated using Fast Fourier Transform (FFT).
+If a filter is approximated this way, one only needs to store its FFT
+parameters instead of all the sampled data points. This would reduce the
+storage size and increase performance, at the cost of reduced accuracy.
+If you decide to use the parameterization functions provided here,
+it is up to you to decide whether the results are good enough for your
+use cases or not.
+
+.. _filter_fft_generation:
+
+Generating FFT
+--------------
+
+.. testsetup::
+
+    >>> import os
+    >>> from astropy.utils.data import get_pkg_data_filename
+    >>> filename = get_pkg_data_filename(
+    ...     os.path.join('data', 'hst_acs_hrc_f555w.fits'),
+    ...     package='synphot.tests')
+
+You could parameterize a given filter using
+:func:`~synphot.filter_parameterization.filter_to_fft` as follows.
+By default, 10 FFT parameters are returned as complex numbers::
+
+    >>> from synphot import SpectralElement
+    >>> from synphot.filter_parameterization import filter_to_fft
+    >>> filename = 'hst_acs_hrc_f555w.fits'  # doctest: +SKIP
+    >>> bp = SpectralElement.from_file(filename)
+    >>> n_lambda, lambda_0, delta_lambda, tr_max, fft_pars = filter_to_fft(bp)
+    >>> n_lambda  # Number of elements in wavelengths
+    10000
+    >>> lambda_0  # Starting value of wavelengths  # doctest: +FLOAT_CMP
+    <Quantity 3479.999 Angstrom>
+    >>> delta_lambda  # Median wavelength separation  # doctest: +FLOAT_CMP
+    <Quantity 0.66748047 Angstrom>
+    >>> tr_max  # Peak value of throughput  # doctest: +FLOAT_CMP
+    <Quantity 0.241445>
+    >>> fft_pars  # FFT parameters  # doctest: +SKIP
+    [(407.5180314841658+7.494005416219807e-16j),
+     (-78.52240189503877-376.53990235136575j),
+     (-294.86589196496584+127.25464850352665j),
+     (130.20273803287864+190.84263652863257j),
+     (96.62299079012317-91.70087676328245j),
+     (-32.572468348727654-34.227696019221035j),
+     (-8.051741476066471-21.354793540998294j),
+     (-51.708676896903725+6.883836090870033j),
+     (13.08719675518801+54.48177212720124j),
+     (38.635087381362396-13.02803811279449j)]
+
+.. TODO: Only skipping the fft_pars comparison above because output is very
+   different for Numpy 1.16. Unskip it and replace with +FLOAT_CMP when
+   Numpy minversion is 1.17.
+
+It is up to you to decide how to store this data, though storing it in a
+table format is recommended. In fact, if you have many filters to parameterize,
+:func:`~synphot.filter_parameterization.filters_to_fft_table`
+will store the results in a table for you::
+
+    >>> from synphot.filter_parameterization import filters_to_fft_table
+    >>> mapping = {'HST/ACS/HRC/F555W': (bp, None)}
+    >>> filter_pars_table = filters_to_fft_table(mapping)
+    >>> filter_pars_table  # doctest: +FLOAT_CMP +ELLIPSIS
+    <Table length=1>
+          filter      n_lambda ...                  fft_9
+                               ...
+          str17        int32   ...                complex128
+    ----------------- -------- ... ---------------------------------------
+    HST/ACS/HRC/F555W    10000 ... (38.635087381362396-13.02803811279449j)
+    >>> filter_pars_table.write('my_filter_pars.fits')  # doctest: +SKIP
+
+.. _filter_fft_construction:
+
+Reconstructing Filter from FFT
+------------------------------
+
+Once you have a parameterized filter (see :ref:`filter_fft_generation`),
+you can reconstruct it for use using
+:func:`~synphot.filter_parameterization.filter_from_fft`.
+Following from the example above::
+
+    >>> from synphot.filter_parameterization import filter_from_fft
+    >>> reconstructed_bp = filter_from_fft(
+    ...     n_lambda, lambda_0, delta_lambda, tr_max, fft_pars)
+
+For this particular example using HST ACS/HRC F555W filter, perhaps 10
+parameters are not quite sufficient. Therefore, caution needs to be exercised
+if you opt to parameterize your filters using this method.
+
+.. plot::
+
+    import os
+    import matplotlib.pyplot as plt
+    from astropy.utils.data import get_pkg_data_filename
+    from synphot import SpectralElement
+    from synphot.filter_parameterization import filter_to_fft, filter_from_fft
+    filename = get_pkg_data_filename(
+        os.path.join('data', 'hst_acs_hrc_f555w.fits'),
+        package='synphot.tests')
+    bp = SpectralElement.from_file(filename)
+    fit_result = filter_to_fft(bp)
+    reconstructed_bp = filter_from_fft(*fit_result)
+    w = bp.waveset
+    plt.plot(w, bp(w), 'b-', label='Original')
+    plt.plot(w, reconstructed_bp(w), 'r--', label='Reconstructed')
+    plt.xlim(3500, 8000)
+    plt.xlabel('Wavelength (Angstrom)')
+    plt.ylabel('Throughput')
+    plt.title('HST ACS/HRC F555W')
+    plt.legend(loc='upper right', numpoints=1)

synphot/filter_parameterization/filter_fft.py

@@ -83,7 +83,12 @@ def filter_to_fft(bp, wavelengths=None, n_terms=10):
     # Take the DFT of the interpolated transmittance curve
     fft = np.fft.fft(tr_interp)[:n_terms]
 
-    return n_lambda, lambda_0, delta_lambda, tr_max, fft.tolist()
+    if isinstance(fft, u.Quantity):
+        fft_parameters = fft.value.tolist()
+    else:  # Older Numpy does not return Quantity
+        fft_parameters = fft.tolist()
+
+    return n_lambda, lambda_0, delta_lambda, tr_max, fft_parameters
 
 
 def filter_from_fft(n_lambda, lambda_0, delta_lambda, tr_max, fft_parameters):
@@ -187,6 +192,8 @@ def filters_to_fft_table(filters_mapping, n_terms=10):
     filters_mapping : dict
         Dictionary mapping human-readable filter name to its
         `~synphot.spectrum.SpectralElement` and wavelengths, if applicable.
+        If the filter object has a valid ``waveset``, just provide `None`
+        for wavelengths; otherwise provide a Quantity array for sampling.
         For example::
 
             {'JOHNSON/V': (<SpectralElement ...>, None),