diff --git a/CHANGES.rst b/CHANGES.rst
index 50ecb8f7..f7cf3484 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -1,6 +1,9 @@
1.2.0 (unreleased)
==================
+- New ``filter_parameterization`` subpackage to handle filter parameterization,
+ adapted from ``tynt`` package written by Brett Morris. [#257]
+
- Dropped support for Python 3.6 and 3.7. Minimum supported Python
version is now 3.8. [#330]
diff --git a/docs/index.rst b/docs/index.rst
index 0801233d..74e45263 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -287,6 +287,7 @@ Using **synphot**
synphot/observation
synphot/formulae
synphot/units
+ synphot/filter_par
synphot/tutorials
.. _synphot_history:
@@ -325,6 +326,8 @@ Also imports this C-extension to local namespace:
.. automodapi:: synphot.exceptions
+.. automodapi:: synphot.filter_parameterization
+
.. automodapi:: synphot.models
.. automodapi:: synphot.observation
diff --git a/docs/synphot/filter_par.rst b/docs/synphot/filter_par.rst
new file mode 100644
index 00000000..699d100b
--- /dev/null
+++ b/docs/synphot/filter_par.rst
@@ -0,0 +1,124 @@
+.. _synphot_par_filters:
+
+Parameterized Filters
+=====================
+
+.. note::
+
+ The algorithm for parameterized filters here was originally developed by
+ Brett Morris for the `tynt `_ package.
+
+Filter responses can 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 reduces the
+storage size and increases 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
+
+ >>> delta_lambda # Median wavelength separation # doctest: +FLOAT_CMP
+
+ >>> tr_max # Peak value of throughput # doctest: +FLOAT_CMP
+
+ >>> 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_LT_1_17. 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: +SKIP
+
+ filter n_lambda ... fft_9
+ ...
+ str17 int... ... complex128
+ ----------------- -------- ... ---------------------------------------
+ HST/ACS/HRC/F555W 10000 ... (38.635087381362396-13.02803811279449j)
+ >>> filter_pars_table.write('my_filter_pars.fits') # doctest: +SKIP
+
+.. TODO: Only skipping the filter_pars_table comparison above because output
+ is slightly different for NUMPY_LT_1_17. Unskip it and replace with
+ +FLOAT_CMP +ELLIPSIS when Numpy minversion is 1.17.
+
+.. _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)
diff --git a/licenses/TYNT_LICENSE.txt b/licenses/TYNT_LICENSE.txt
new file mode 100644
index 00000000..6a3fc282
--- /dev/null
+++ b/licenses/TYNT_LICENSE.txt
@@ -0,0 +1,29 @@
+BSD 3-Clause License
+
+Copyright (c) 2019, Brett M. Morris
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+1. Redistributions of source code must retain the above copyright notice, this
+ list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright notice,
+ this list of conditions and the following disclaimer in the documentation
+ and/or other materials provided with the distribution.
+
+3. Neither the name of the copyright holder nor the names of its
+ contributors may be used to endorse or promote products derived from
+ this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
diff --git a/setup.cfg b/setup.cfg
index 7addd4ec..9cf8b6fb 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -61,6 +61,7 @@ docs =
[options.package_data]
synphot.tests = data/*
+synphot.filter_parameterization.tests = data/*
[flake8]
# Ignoring these for now:
diff --git a/synphot/compat.py b/synphot/compat.py
index 90e189f3..c72c9d7b 100644
--- a/synphot/compat.py
+++ b/synphot/compat.py
@@ -2,6 +2,7 @@
"""Module to handle backward-compatibility."""
import astropy
+import numpy
from astropy.utils.introspection import minversion
try:
@@ -18,11 +19,12 @@
else:
HAS_DUST_EXTINCTION = True
-
__all__ = ['ASTROPY_LT_5_0', 'ASTROPY_LT_4_3', 'ASTROPY_LT_4_1',
- 'ASTROPY_LT_4_0', 'HAS_SPECUTILS', 'HAS_DUST_EXTINCTION']
+ 'ASTROPY_LT_4_0', 'NUMPY_LT_1_17', 'HAS_SPECUTILS',
+ 'HAS_DUST_EXTINCTION']
ASTROPY_LT_5_0 = not minversion(astropy, '4.99') # astropy<5 but includes 5.0.dev # noqa
ASTROPY_LT_4_3 = not minversion(astropy, '4.3')
ASTROPY_LT_4_1 = not minversion(astropy, '4.1')
ASTROPY_LT_4_0 = not minversion(astropy, '4.0')
+NUMPY_LT_1_17 = not minversion(numpy, '1.17')
diff --git a/synphot/filter_parameterization/__init__.py b/synphot/filter_parameterization/__init__.py
new file mode 100644
index 00000000..7d05f9c3
--- /dev/null
+++ b/synphot/filter_parameterization/__init__.py
@@ -0,0 +1,9 @@
+"""This subpackage handles filter parameterization.
+
+The algorithms in this subpackage were originally developed by
+Brett Morris as part of the `tynt `_
+package.
+
+"""
+
+from .filter_fft import * # noqa
diff --git a/synphot/filter_parameterization/filter_fft.py b/synphot/filter_parameterization/filter_fft.py
new file mode 100644
index 00000000..909ec889
--- /dev/null
+++ b/synphot/filter_parameterization/filter_fft.py
@@ -0,0 +1,233 @@
+"""Handle Fast Fourier Transform (FFT) for filter parameterization."""
+
+import numpy as np
+from astropy import units as u
+from astropy.modeling.models import custom_model, Sine1D
+from astropy.table import Table
+
+from synphot.compat import NUMPY_LT_1_17
+from synphot.models import Empirical1D
+from synphot.spectrum import SpectralElement
+from synphot.units import validate_quantity
+
+__all__ = ['filter_to_fft', 'filter_from_fft', 'analytical_model_from_fft',
+ 'filters_to_fft_table']
+
+
+def _simplified_wavelength(n_lambda, lambda_0, delta_lambda):
+ # tynt assumed everything was in Angstrom, which coincides with
+ # synphot internal wavelength unit.
+ wave_unit = SpectralElement._internal_wave_unit
+
+ lambda_0 = validate_quantity(
+ lambda_0, wave_unit, equivalencies=u.spectral())
+ delta_lambda = validate_quantity(
+ delta_lambda, wave_unit, equivalencies=u.spectral())
+ lambda_max = (n_lambda + 1) * delta_lambda + lambda_0
+
+ return np.arange(lambda_0.value, lambda_max.value,
+ delta_lambda.value) * wave_unit
+
+
+def filter_to_fft(bp, wavelengths=None, n_terms=10):
+ """Calculate filter parameters using FFT.
+
+ Parameters
+ ----------
+ bp : `~synphot.spectrum.SpectralElement`
+ Filter to parameterize.
+
+ wavelengths : array-like or `~astropy.units.quantity.Quantity`
+ Wavelength values for sampling.
+ If not a Quantity, assumed to be in Angstrom.
+ If `None`, ``waveset`` is used.
+
+ n_terms : int
+ Number of FFT parameters to keep.
+
+ Returns
+ -------
+ n_lambda : int
+ Number of elements in ``wl``.
+
+ lambda_0 : `~astropy.units.quantity.Quantity`
+ Minimum value of ``wl``.
+
+ delta_lambda : `~astropy.units.quantity.Quantity`
+ Median delta wavelength.
+
+ tr_max : `~astropy.units.quantity.Quantity`
+ Maximum value of ``tr``.
+
+ fft_parameters : list of complex
+ List of complex values that are FFT parameters to keep.
+
+ """
+ wl = bp._validate_wavelengths(wavelengths)
+ tr = bp(wl)
+
+ diff_wl = np.diff(wl)
+
+ delta_lambda = np.nanmedian(diff_wl[diff_wl != 0])
+ lambda_0 = wl.min()
+ n_lambda = len(wl)
+
+ # Create a simplified wavelength grid
+ simplified_wavelength = _simplified_wavelength(
+ n_lambda, lambda_0, delta_lambda)
+
+ tr_max = tr.max()
+
+ # Interpolate transmittance onto simplified wavelength grid
+ tr_interp = np.interp(simplified_wavelength, wl, tr)
+
+ # Take the DFT of the interpolated transmittance curve
+ fft = np.fft.fft(tr_interp)[:n_terms]
+
+ 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):
+ """Reconstruct a filter from given FFT parameters.
+ The inputs for this function can be obtained from :func:`filter_to_fft`.
+
+ Parameters
+ ----------
+ n_lambda : int
+ Number of elements in original wavelength array.
+
+ lambda_0 : float or `~astropy.units.quantity.Quantity`
+ Minimum value of original wavelength array.
+ If not a Quantity, assumed to be in Angstrom.
+
+ delta_lambda : float or `~astropy.units.quantity.Quantity`
+ Median delta wavelength of original wavelength array.
+ If not a Quantity, assumed to be in Angstrom.
+
+ tr_max : float or `~astropy.units.quantity.Quantity`
+ Maximum value of transmittance curve.
+ If a Quantity, must be unitless.
+
+ fft_parameters : list of complex
+ List of complex values that are FFT parameters representing the
+ filter transmittance curve.
+
+ Returns
+ -------
+ bp : `~synphot.spectrum.SpectralElement`
+ Reconstructed filter.
+
+ """
+ wavelength = _simplified_wavelength(n_lambda, lambda_0, delta_lambda)
+ n_wave = len(wavelength)
+ ifft = np.fft.ifft(fft_parameters, n=n_wave)
+ transmittance = ((ifft.real - ifft.real.min()) * tr_max / ifft.real.ptp()) # noqa
+ return SpectralElement(
+ Empirical1D, points=wavelength, lookup_table=transmittance)
+
+
+def analytical_model_from_fft(n_lambda, lambda_0, delta_lambda, tr_max,
+ fft_parameters):
+ """Similar to :func:`filter_from_fft` except that this returns
+ an analytical model.
+
+ .. note::
+
+ This model needs to be sampled using the full range of
+ wavelength. See https://github.com/bmorris3/tynt/issues/9 .
+
+ Returns
+ -------
+ astropy_model : `~astropy.modeling.CompoundModel`
+ A compound model that consists of
+ `~astropy.modeling.functional_models.Sine1D` models.
+
+ """
+ wavelength = _simplified_wavelength(n_lambda, lambda_0, delta_lambda)
+ n_wave = len(wavelength)
+ n_fft_pars = len(fft_parameters)
+
+ m = (np.sum([Sine1D(amplitude=fft_parameters[i].real / n_wave,
+ frequency=i / n_wave, phase=0.25)
+ for i in range(n_fft_pars)]) -
+ np.sum([Sine1D(amplitude=fft_parameters[i].imag / n_wave,
+ frequency=i / n_wave)
+ for i in range(n_fft_pars)]))
+
+ @custom_model
+ def fft_model(x):
+ """Approximate Fourier reconstruction of an astronomical filter.
+
+ Parameters
+ ----------
+ x : `~astropy.units.quantity.Quantity`
+ Full wavelength range that samples the filter.
+
+ Returns
+ -------
+ transmittance : array-like or `~astropy.units.quantity.Quantity`
+ Transmittance curve. If ``tr_max`` is a Quantity, this will
+ be a Quantity as well.
+
+ """
+ wave_unit = SpectralElement._internal_wave_unit
+ x = validate_quantity(x, wave_unit, equivalencies=u.spectral())
+
+ mo = m((x - wavelength.min()) / (wavelength[1] - wavelength[0]))
+ return (mo - mo.min()) * tr_max / mo.ptp()
+
+ return fft_model()
+
+
+def filters_to_fft_table(filters_mapping, n_terms=10):
+ """Run :func:`filter_to_fft` on a list of filters
+ and store results in a table.
+
+ Parameters
+ ----------
+ 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': (, None),
+ 'Flat': (, )}
+
+ n_terms : int
+ Number of FFT parameters to keep.
+
+ Returns
+ -------
+ fft_table : `~astropy.table.Table`
+ Table storing FFT parameterization for the given filters.
+ Use its ``write`` method to save it to file.
+
+ """ # noqa
+ wave_unit = SpectralElement._internal_wave_unit
+ colnames = ['filter', 'n_lambda', 'lambda_0', 'delta_lambda',
+ 'tr_max'] + [f'fft_{i}' for i in range(n_terms)]
+ rows = []
+
+ for key, (bp, wavelengths) in filters_mapping.items():
+ n_lambda, lambda_0, delta_lambda, tr_max, fft_pars = filter_to_fft(
+ bp, wavelengths=wavelengths, n_terms=n_terms)
+ if not NUMPY_LT_1_17:
+ rows.append(tuple(
+ [key, n_lambda, lambda_0, delta_lambda, tr_max] + fft_pars))
+ else: # Numpy 1.16 cannot handle unit here
+ rows.append(tuple(
+ [key, n_lambda, lambda_0.value, delta_lambda.value,
+ tr_max.value] + fft_pars))
+
+ fft_table = Table(rows=rows, names=colnames)
+ fft_table['lambda_0'].unit = wave_unit
+ fft_table['delta_lambda'].unit = wave_unit
+
+ return fft_table
diff --git a/synphot/filter_parameterization/tests/__init__.py b/synphot/filter_parameterization/tests/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/synphot/filter_parameterization/tests/data/fft_test_data.fits b/synphot/filter_parameterization/tests/data/fft_test_data.fits
new file mode 100644
index 00000000..10860302
--- /dev/null
+++ b/synphot/filter_parameterization/tests/data/fft_test_data.fits
@@ -0,0 +1,7 @@
+SIMPLE = T / conforms to FITS standard BITPIX = 8 / array data type NAXIS = 0 / number of array dimensions EXTEND = T END XTENSION= 'BINTABLE' / binary table extension BITPIX = 8 / array data type NAXIS = 2 / number of array dimensions NAXIS1 = 193 / length of dimension 1 NAXIS2 = 13 / length of dimension 2 PCOUNT = 0 / number of group parameters GCOUNT = 1 / number of groups TFIELDS = 15 / number of table fields TTYPE1 = 'filter ' TFORM1 = '17A ' TTYPE2 = 'n_lambda' TFORM2 = 'J ' TTYPE3 = 'lambda_0' TFORM3 = 'E ' TUNIT3 = 'Angstrom' TTYPE4 = 'delta_lambda' TFORM4 = 'E ' TUNIT4 = 'Angstrom' TTYPE5 = 'tr_max ' TFORM5 = 'E ' TTYPE6 = 'fft_0 ' TFORM6 = 'M ' TTYPE7 = 'fft_1 ' TFORM7 = 'M ' TTYPE8 = 'fft_2 ' TFORM8 = 'M ' TTYPE9 = 'fft_3 ' TFORM9 = 'M ' TTYPE10 = 'fft_4 ' TFORM10 = 'M ' TTYPE11 = 'fft_5 ' TFORM11 = 'M ' TTYPE12 = 'fft_6 ' TFORM12 = 'M ' TTYPE13 = 'fft_7 ' TFORM13 = 'M ' TTYPE14 = 'fft_8 ' TFORM14 = 'M ' TTYPE15 = 'fft_9 ' TFORM15 = 'M ' END SLOAN/SDSS.u /E:@ A =5?@cS@ rP" ?>
+οmc?Ne2)?$?[6 ? ,V}sT~ٿʾטqj3# `?
+8ah?oKΣ0SLOAN/SDSS.g YEb A >P@0;x !WUǁ|BP_lI<"M6@5؆Y "C[J5S 7@?se$?oؒX+[(?{nI:մ_A Rt?ڼ? @?_0܀?+R,X_@?@@#U?ރ+.SLOAN/SDSS.i YE A >غ@1?# "-mܿbE3_Q? 0nz:@H;E?x=R?Q -j?Ȣ6傭{8GxcMkD?؈r$?Zhkn'?àHԳ?`piL['IHWSLOAN/SDSS.z E A =-w@RT_ =N i> ȋ~?
+L$? Q?vI>~?b.u?H~?rܿ;x@G~ .{?/^TEBEQ*&Ez@C"'2?z?62MASS/2MASS.J kF% A ? @J'6 3֕fS@!=<^I@Jk@"ז$|K`Q^@wãT?2W?<{ȑq-G?_?]r)?겎?'=8?0 0]ª2MASS/2MASS.H :FIh B ? @;֕7 >$@/V.;;KnB4?ܙp$f ?5*'?p??@\?2Aw}h?\0? <