Spectral Operations (spectral_ops)

Overview

The spectral_ops subpackage provides tools for performing spectral analysis and manipulation, specifically focusing on the convolution of spectral data with slit functions. This subpackage is designed to facilitate the preprocessing and analysis of spectral data by applying mathematical transformations that mimic the effects of instrumental slit functions on measured spectra.

Features

Convolution with Slit Functions: Apply a specified slit function to spectral data, simulating the instrumental broadening effects.
Support for Equidistant and Non-Equidistant Spectral Data: Functions within the subpackage can handle both equidistant and non-equidistant spectral wavelength grids, providing flexibility for various types of spectral data.
Spline Representation for Slit Functions: Utilize spline representations (specifically, PPoly objects from SciPy) to define slit functions, allowing for precise and customizable convolution operations.

Reference

To use the spectral_ops subpackage, import the desired function(s) and pass your spectral data along with the slit function spline. Ensure your spectral data and slit function are prepared according to the function requirements, particularly regarding the equidistance of wavelength grids and the spline representation of the slit function.

spec_ops: Provide several convolution of spectrum with a polynomial.

`conv_eqi_spec_with_slit(input_wavelengths, input_intensities, output_wavelengths, slit_function_spline)`

Convolve the input spectrum with the slit function described by a PPoly spline onto the specified output wavelengths.

Parameters:

Name	Type	Description	Default
`input_wavelengths`	`ndarray`	Wavelengths of the input spectrum.	required
`input_intensities`	`ndarray`	Intensity values of the input spectrum.	required
`output_wavelengths`	`ndarray`	Output wavelengths for the convolved spectrum.	required
`slit_function_spline`	`PPoly`	Spline representation of the slit function. Use `scipy.interpolate`	required

Returns:

Type	Description
`ndarray`	Intensity values of the convolved spectrum at output wavelengths in as numpy `ndarray`.

Source code in src/pyshic/spectral_ops/spec_ops.py

def conv_eqi_spec_with_slit(input_wavelengths: np.ndarray, 
                            input_intensities: np.ndarray, 
                            output_wavelengths: np.ndarray, 
                            slit_function_spline: PPoly) -> np.ndarray:
    """
    Convolve the input spectrum with the slit function described by a PPoly spline
    onto the specified output wavelengths.

    Args:
        input_wavelengths (np.ndarray): Wavelengths of the input spectrum.
        input_intensities (np.ndarray): Intensity values of the input spectrum.
        output_wavelengths (np.ndarray): Output wavelengths for the convolved spectrum.
        slit_function_spline (PPoly): Spline representation of the slit function. 
                                    Use `scipy.interpolate`

    Returns:
        Intensity values of the convolved spectrum at output wavelengths in as numpy `ndarray`.
    """
    # Extracting spline properties
    breaks = slit_function_spline.x
    spmin = np.min(breaks)
    spmax = np.max(breaks)

    # Check for equidistant wavelength vectors
    d1 = np.diff(input_wavelengths)
    d2 = np.diff(output_wavelengths)
    if not np.allclose(d1, d1[0], atol=1e-6) or not np.allclose(d2, d2[0], atol=1e-6):
        raise ValueError('Wavelength vectors are not equidistant.')

    # Check step size requirement
    min_d1 = np.min(d1)
    min_d2 = np.min(d2)
    period = min_d2 / min_d1
    if not np.isclose(period, round(period), atol=1e-6):
        raise ValueError('Step size of output_wavelengths must be a multiple of the step size of input_wavelengths.')
    period = int(round(period))

    offset = (output_wavelengths[0] - input_wavelengths[0]) / min_d1
    delta = (offset - np.floor(offset)) * min_d1
    starting_index = 1 + np.floor(offset)

    slit_start = np.floor((spmin + delta) / min_d1) - 1
    slit_end = np.ceil((spmax + delta) / min_d1) + 1
    slit_wavelengths = np.arange(slit_end, slit_start - 1, -1) * min_d1 - delta

    # Evaluate spline
    slit_values = slit_function_spline(slit_wavelengths) * (slit_wavelengths >= spmin) * (slit_wavelengths <= spmax)
    slit_values = slit_values / np.sum(slit_values)

    starting_index += slit_end
    convolved = convolve(input_intensities, slit_values, mode='full')
    indices = np.arange(0, len(output_wavelengths)) * period + starting_index
    convolved_intensities = convolved[indices.astype(int)-1]

    return convolved_intensities

`conv_spec_with_slit(spec_wavelengths, spec_intensities, output_wavelengths, slit_spline)`

Convolve spectrum with a slit function on specified output wavelengths.

Parameters:

Name	Type	Description	Default
`spec_wavelengths`	`ndarray`	Wavelengths of the input spectrum.	required
`spec_intensities`	`ndarray`	Intensity values of the input spectrum.	required
`output_wavelengths`	`ndarray`	Output wavelengths for the convolved spectrum.	required
`slit_spline`	`PPoly`	Spline representation of the slit function. Use `scipy.interpolate`	required

Returns: Intensity values of the convolved spectrum at output wavelengths.

Source code in src/pyshic/spectral_ops/spec_ops.py

def conv_spec_with_slit(spec_wavelengths: np.ndarray, 
                        spec_intensities: np.ndarray, 
                        output_wavelengths: np.ndarray, 
                        slit_spline: PPoly) -> np.ndarray:
    """
    Convolve spectrum with a slit function on specified output wavelengths.

    Args:
        spec_wavelengths (np.ndarray): Wavelengths of the input spectrum.
        spec_intensities (np.ndarray): Intensity values of the input spectrum.
        output_wavelengths (np.ndarray): Output wavelengths for the convolved spectrum.
        slit_spline (PPoly): Spline representation of the slit function.
                             Use `scipy.interpolate`

    Returns:
    Intensity values of the convolved spectrum at output wavelengths.
    """
    # Check if wavelength vectors are equidistant and meet step size requirements
    diff_spec_wl = np.diff(spec_wavelengths)
    diff_output_wl = np.diff(output_wavelengths)
    min_diff_spec_wl = np.min(diff_spec_wl)
    min_diff_output_wl = np.min(diff_output_wl)
    period = min_diff_output_wl / min_diff_spec_wl

    if (np.max(diff_spec_wl) - np.min(diff_spec_wl) > 1e-6 or
        np.max(diff_output_wl) - np.min(diff_output_wl) > 1e-6 or
        abs(period - round(period)) > 1e-6):
        eqi = False
    else:
        eqi = True

    valid_indices = ~np.isnan(spec_intensities)

    if eqi:
        convolved_intensities = conv_eqi_spec_with_slit(spec_wavelengths, spec_intensities, output_wavelengths, slit_spline)
    else:
        # For non-equidistant grid: convolve then interpolate
        temp_convolution = conv_eqi_spec_with_slit(spec_wavelengths[valid_indices], spec_intensities[valid_indices], spec_wavelengths, slit_spline)
        cs = CubicSpline(spec_wavelengths[valid_indices], temp_convolution)
        convolved_intensities = cs(output_wavelengths)

    return convolved_intensities

Examples

from spectral_ops import conv_spec_with_slit
import numpy as np
from scipy.interpolate import CubicSpline

# Example spectral data and slit function
spec_wavelengths = np.linspace(400, 700, 300)
spec_intensities = np.sin(2 * np.pi * spec_wavelengths / 100) + 0.5
output_wavelengths = np.linspace(400, 700, 150)
x = np.linspace(-5, 5, 100)
y = np.exp(-x**2)
slit_spline = CubicSpline(x, y)

# Perform convolution
convolved_intensities = conv_spec_with_slit(spec_wavelengths, spec_intensities, output_wavelengths, slit_spline)

For detailed examples and advanced usage, refer to the function documentation within the subpackage.

This README template provides a basic introduction to the spectral_ops subpackage, its mathematical foundations, and a simple usage example. You can expand this documentation with more specific examples, installation instructions, and any dependencies or prerequisites as needed for your package.

MATH

The main mathematical operation performed by this subpackage is the convolution of spectral data with a specified slit function. Convolution is a mathematical operation that produces a new function showing how the shape of one function is modified by the other. In the context of spectral analysis, convolution with a slit function simulates how an instrument's finite resolution broadens the spectral lines.

Convolution Operation

Given a spectrum \(I(\lambda)\) and a slit function \(S(\lambda)\), the convolution operation is defined as:

\[ I'(\lambda) = (I * S)(\lambda) = \int_{-\infty}^{\infty} I(\lambda') S(\lambda - \lambda') d\lambda' \]

where \(I'(\lambda)\) is the convolved spectrum, showing the instrumental broadening effect on the original spectral data.

Spline Representation of Slit Functions

Slit functions are represented using spline interpolation (specifically, piecewise polynomials or PPoly in SciPy), which allows for a flexible and accurate description of the slit function's shape. This representation is crucial for simulating various instrumental profiles accurately.

Handling of Equidistant and Non-Equidistant Grids

The subpackage accommodates both equidistant and non-equidistant spectral grids by adjusting the convolution process accordingly. For non-equidistant grids, additional steps such as interpolation are performed to ensure the convolution operation accurately reflects the instrumental effects on the spectral data.