From 440823b750f80b99dc6e4e1921ea773a490fef66 Mon Sep 17 00:00:00 2001
From: Pey Lian Lim <2090236+pllim@users.noreply.github.com>
Date: Tue, 5 May 2020 18:05:42 -0400
Subject: [PATCH] Add documentation. Fix bugs.
---
docs/index.rst | 1 +
docs/synphot/filter_par.rst | 117 ++++++++++++++++++
synphot/filter_parameterization/filter_fft.py | 4 +-
3 files changed, 121 insertions(+), 1 deletion(-)
create mode 100644 docs/synphot/filter_par.rst
diff --git a/docs/index.rst b/docs/index.rst
index 086c702d..3089d270 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -278,6 +278,7 @@ Using **synphot**
synphot/observation
synphot/formulae
synphot/units
+ synphot/filter_fft
synphot/tutorials
.. _synphot_history:
diff --git a/docs/synphot/filter_par.rst b/docs/synphot/filter_par.rst
new file mode 100644
index 00000000..6efe71d6
--- /dev/null
+++ b/docs/synphot/filter_par.rst
@@ -0,0 +1,117 @@
+.. _synphot_par_filters:
+
+Parameterized Filters
+=====================
+
+.. note::
+
+ The algorithm for parameterized filters here was originally developed by
+ Brett Morris for the `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_fft.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: +FLOAT_CMP
+ [(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)]
+
+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.filter_fft.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
+
+ 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_fft.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::
+ :include-source:
+
+ 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/synphot/filter_parameterization/filter_fft.py b/synphot/filter_parameterization/filter_fft.py
index 8f8731c7..4dd8ab03 100644
--- a/synphot/filter_parameterization/filter_fft.py
+++ b/synphot/filter_parameterization/filter_fft.py
@@ -83,7 +83,7 @@ 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()
+ return n_lambda, lambda_0, delta_lambda, tr_max, fft.value.tolist()
def filter_from_fft(n_lambda, lambda_0, delta_lambda, tr_max, fft_parameters):
@@ -187,6 +187,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': (, None),