Source code for synphot.filter_parameterization.filter_fft

"""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.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',

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

[docs] 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
[docs] 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 / np.ptp(ifft.real)) # noqa return SpectralElement( Empirical1D, points=wavelength, lookup_table=transmittance)
[docs] 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 . 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 / np.ptp(mo) return fft_model()
[docs] 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': (<SpectralElement ...>, None), 'Flat': (<SpectralElement ...>, <Quantity [1000., ..., 9999.] Angstrom>)} 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) rows.append(tuple( [key, n_lambda, lambda_0, delta_lambda, tr_max] + 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