"""
Spectrum convolution with lineshapes
Modern implementation with multi-dimensional support
"""
import matplotlib.pyplot as plt
import numpy as np
from scipy import signal
from ..logging_config import get_logger
logger = get_logger(__name__)
[docs]
def convspec(spectrum, step_size, width, derivative=0, alpha=1.0, phase=0.0):
"""
Convolve spectrum with lineshape functions.
Applies broadening to stick spectra or other discrete data by
convolution with Gaussian, Lorentzian, or pseudo-Voigt profiles.
Parameters:
-----------
spectrum : array
Input spectrum to convolve
step_size : float or array
Abscissa step size for each dimension
width : float or array
Full width at half maximum for lineshape
derivative : int or array, default=0
Derivative order (0=function, 1=first deriv, 2=second deriv)
alpha : float or array, default=1.0
Shape parameter (1=Gaussian, 0=Lorentzian, 0-1=pseudo-Voigt)
phase : float or array, default=0.0
Phase (0=absorption, π/2=dispersion)
Returns:
--------
array
Convolved spectrum with same shape as input
Examples:
---------
>>> # Simple 1D convolution
>>> x = np.linspace(0, 100, 1000)
>>> stick_spec = np.zeros_like(x)
>>> stick_spec[500] = 1.0 # Delta peak at center
>>> broadened = convspec(stick_spec, 0.1, 2.0) # Gaussian, FWHM=2
>>>
>>> # Lorentzian broadening
>>> lorentz = convspec(stick_spec, 0.1, 2.0, alpha=0.0)
>>>
>>> # First derivative
>>> deriv = convspec(stick_spec, 0.1, 2.0, derivative=1)
"""
spectrum = np.asarray(
spectrum, dtype=complex if np.iscomplexobj(spectrum) else float
)
# Handle multi-dimensional parameters
ndim = spectrum.ndim
step_size = _expand_parameter(step_size, ndim)
width = _expand_parameter(width, ndim)
derivative = _expand_parameter(derivative, ndim)
alpha = _expand_parameter(alpha, ndim)
phase = _expand_parameter(phase, ndim)
# Validate inputs
_validate_convspec_inputs(spectrum, step_size, width, derivative, alpha, phase)
# Perform convolution
result = _convolve_spectrum(spectrum, step_size, width, derivative, alpha, phase)
# Preserve real/complex nature of input
if np.isrealobj(spectrum) and np.iscomplexobj(result):
result = np.real(result)
return result
def _expand_parameter(param, ndim):
"""Expand scalar parameters to match number of dimensions"""
param = np.asarray(param)
if param.ndim == 0:
return np.full(ndim, param.item())
elif len(param) == ndim:
return param
else:
raise ValueError(
f"Parameter length {len(param)} doesn't match spectrum dimensions {ndim}"
)
def _validate_convspec_inputs(spectrum, step_size, width, derivative, alpha, phase):
"""Validate convolution parameters"""
if np.any(step_size <= 0):
raise ValueError("step_size must be positive")
if np.any(width < 0):
raise ValueError("width must be non-negative")
if np.any((derivative < -1) | (derivative > 2)):
raise ValueError("derivative must be -1, 0, 1, or 2")
if np.any((alpha < 0) | (alpha > 1)):
raise ValueError("alpha must be between 0 and 1")
if np.any(~np.isfinite([step_size, width, derivative, alpha, phase])):
raise ValueError("All parameters must be finite")
def _convolve_spectrum(spectrum, step_size, width, derivative, alpha, phase):
"""Core convolution implementation using FFT"""
# For simplicity, use basic convolution with scipy.signal
# Full implementation would use lshape function
from .gaussian import gaussian
from .lorentzian import lorentzian
# Create convolution kernel
n_kernel = min(len(spectrum) // 2, 500) # Reasonable kernel size
x_kernel = np.arange(-n_kernel, n_kernel + 1) * step_size[0]
if alpha[0] == 1.0:
# Pure Gaussian
kernel = gaussian(
x_kernel, 0, width[0], derivative=int(derivative[0]), phase=phase[0]
)
elif alpha[0] == 0.0:
# Pure Lorentzian
kernel = lorentzian(
x_kernel, 0, width[0], derivative=int(derivative[0]), phase=phase[0]
)
else:
# Mixed (pseudo-Voigt)
gauss_part = gaussian(
x_kernel, 0, width[0], derivative=int(derivative[0]), phase=phase[0]
)
lorentz_part = lorentzian(
x_kernel, 0, width[0], derivative=int(derivative[0]), phase=phase[0]
)
kernel = alpha[0] * gauss_part + (1 - alpha[0]) * lorentz_part
# Normalize kernel
if derivative[0] == 0 and np.sum(kernel) != 0:
kernel = kernel / np.sum(kernel)
# Convolve
result = signal.convolve(spectrum, kernel, mode="same")
return result
def demo():
"""Interactive demonstration of spectrum convolution"""
# Modern colors
colors = ["#e41a1c", "#377eb8", "#4daf4a", "#984ea3", "#ff7f00"]
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
# 1D convolution examples
ax = axes[0, 0]
# Create stick spectrum
x = np.linspace(0, 100, 1000)
stick = np.zeros_like(x)
stick[200] = 1.0 # Peak at x=20
stick[500] = 1.5 # Peak at x=50
stick[800] = 0.8 # Peak at x=80
# Different convolutions
try:
gauss = convspec(stick, 0.1, 4.0, alpha=1.0)
lorentz = convspec(stick, 0.1, 4.0, alpha=0.0)
voigt = convspec(stick, 0.1, 4.0, alpha=0.5)
ax.plot(x, stick, "k-", linewidth=3, alpha=0.7, label="Stick spectrum")
ax.plot(x, gauss, color=colors[0], linewidth=2, label="Gaussian")
ax.plot(x, lorentz, color=colors[1], linewidth=2, label="Lorentzian")
ax.plot(x, voigt, color=colors[2], linewidth=2, label="Pseudo-Voigt")
except ImportError:
# Fallback if lineshapes not available
from scipy import ndimage
gauss = ndimage.gaussian_filter1d(
stick, sigma=4.0 / 0.1 / (2 * np.sqrt(2 * np.log(2)))
)
ax.plot(x, stick, "k-", linewidth=3, alpha=0.7, label="Stick spectrum")
ax.plot(x, gauss, color=colors[0], linewidth=2, label="Gaussian (scipy)")
ax.set_title("Spectrum Convolution", fontweight="bold")
ax.legend()
ax.grid(alpha=0.3)
# Simple broadening demonstration
ax = axes[0, 1]
# Single peak
single_peak = np.zeros_like(x)
single_peak[500] = 1.0
# Different widths
widths = [2, 4, 8]
for i, width in enumerate(widths):
try:
convolved = convspec(single_peak, 0.1, width)
ax.plot(
x, convolved, color=colors[i], linewidth=2.5, label=f"FWHM = {width}"
)
except Exception:
# Fallback
from scipy import ndimage
sigma = width / 0.1 / (2 * np.sqrt(2 * np.log(2)))
convolved = ndimage.gaussian_filter1d(single_peak, sigma=sigma)
ax.plot(
x, convolved, color=colors[i], linewidth=2.5, label=f"FWHM ≈ {width}"
)
ax.set_title("Different Widths", fontweight="bold")
ax.legend()
ax.grid(alpha=0.3)
# Conceptual comparison of shapes
ax = axes[1, 0]
# Create theoretical shapes for comparison
x_theory = np.linspace(-10, 10, 200)
# Gaussian
sigma = 2 / (2 * np.sqrt(2 * np.log(2)))
gauss_theory = np.exp(-(x_theory**2) / (2 * sigma**2)) / (
sigma * np.sqrt(2 * np.pi)
)
# Lorentzian
gamma = 1 # half-width
lorentz_theory = (1 / np.pi) * gamma / (x_theory**2 + gamma**2)
# Pseudo-Voigt (50/50 mix)
voigt_theory = 0.5 * gauss_theory + 0.5 * lorentz_theory
ax.plot(x_theory, gauss_theory, color=colors[0], linewidth=2.5, label="Gaussian")
ax.plot(
x_theory, lorentz_theory, color=colors[1], linewidth=2.5, label="Lorentzian"
)
ax.plot(
x_theory, voigt_theory, color=colors[2], linewidth=2.5, label="Pseudo-Voigt"
)
ax.set_title("Lineshape Comparison", fontweight="bold")
ax.legend()
ax.grid(alpha=0.3)
# Usage example
ax = axes[1, 1]
# Create EPR-like spectrum with multiple peaks
positions = [25, 35, 45, 55, 65, 75] # Peak positions
intensities = [0.8, 1.2, 1.0, 0.6, 1.1, 0.9] # Relative intensities
epr_stick = np.zeros_like(x)
for pos, intensity in zip(positions, intensities):
idx = int(pos * 10) # Convert to array index
if 0 <= idx < len(epr_stick):
epr_stick[idx] = intensity
# Apply broadening
try:
epr_broadened = convspec(epr_stick, 0.1, 3.0, alpha=0.7) # Mostly Gaussian
ax.plot(x, epr_broadened, color=colors[1], linewidth=2.5, label="Broadened")
except Exception:
from scipy import ndimage
sigma = 3.0 / 0.1 / (2 * np.sqrt(2 * np.log(2)))
epr_broadened = ndimage.gaussian_filter1d(epr_stick, sigma=sigma)
ax.plot(x, epr_broadened, color=colors[1], linewidth=2.5, label="Broadened")
# Show stick spectrum
ax.stem(
x[epr_stick > 0],
epr_stick[epr_stick > 0],
linefmt="gray",
markerfmt="ko",
basefmt=" ",
label="Stick spectrum",
)
ax.set_title("EPR-like Spectrum", fontweight="bold")
ax.legend()
ax.grid(alpha=0.3)
# Style all subplots
for ax in axes.flat:
ax.set_xlabel("Position/Field (mT)")
ax.set_ylabel("Intensity")
ax.spines["top"].set_visible(False)
ax.spines["right"].set_visible(False)
plt.tight_layout()
plt.show()
logger.info("\nConvolution Demo:")
logger.info("- convspec() applies lineshape broadening to spectra")
logger.info("- Converts stick spectra to realistic lineshapes")
logger.info("- Supports Gaussian, Lorentzian, and pseudo-Voigt profiles")
logger.info("- Essential for EPR spectrum simulation")
if __name__ == "__main__":
demo()