Derivatives and Peak Analysis¶

This notebook demonstrates:

1. Computing spectral derivatives (Savitzky-Golay and gap-segment)
2. Finding peaks in original and derivative spectra
3. Integrating peak areas

In [ ]:

import matplotlib.pyplot as plt
import numpy as np

from spectrakit import (
    derivative_gap_segment,
    derivative_savgol,
    peaks_find,
    peaks_integrate,
    smooth_savgol,
)

import matplotlib.pyplot as plt import numpy as np from spectrakit import ( derivative_gap_segment, derivative_savgol, peaks_find, peaks_integrate, smooth_savgol, )

Create Test Spectrum with Overlapping Peaks¶

In [ ]:

wavenumbers = np.linspace(400, 4000, 1000)


def gaussian(x, c, a, s):
    return a * np.exp(-((x - c) ** 2) / (2 * s**2))


# Two overlapping peaks + one isolated peak
spectrum = (
    gaussian(wavenumbers, 1600, 2.0, 30)
    + gaussian(wavenumbers, 1680, 1.5, 25)  # overlapping with first
    + gaussian(wavenumbers, 2900, 3.0, 50)  # isolated
)

# Add some noise
rng = np.random.default_rng(42)
noisy = spectrum + rng.normal(0, 0.03, len(spectrum))

plt.figure(figsize=(10, 4))
plt.plot(wavenumbers, noisy, alpha=0.7, label="Noisy")
plt.plot(wavenumbers, spectrum, "--", label="True")
plt.legend()
plt.xlabel("Wavenumber")
plt.gca().invert_xaxis()
plt.title("Test Spectrum with Overlapping Peaks")
plt.show()

wavenumbers = np.linspace(400, 4000, 1000) def gaussian(x, c, a, s): return a * np.exp(-((x - c) ** 2) / (2 * s**2)) # Two overlapping peaks + one isolated peak spectrum = ( gaussian(wavenumbers, 1600, 2.0, 30) + gaussian(wavenumbers, 1680, 1.5, 25) # overlapping with first + gaussian(wavenumbers, 2900, 3.0, 50) # isolated ) # Add some noise rng = np.random.default_rng(42) noisy = spectrum + rng.normal(0, 0.03, len(spectrum)) plt.figure(figsize=(10, 4)) plt.plot(wavenumbers, noisy, alpha=0.7, label="Noisy") plt.plot(wavenumbers, spectrum, "--", label="True") plt.legend() plt.xlabel("Wavenumber") plt.gca().invert_xaxis() plt.title("Test Spectrum with Overlapping Peaks") plt.show()

Savitzky-Golay Derivatives¶

In [ ]:

# Smooth first to reduce noise amplification
smoothed = smooth_savgol(noisy, window_length=15, polyorder=3)

# First and second derivatives
d1 = derivative_savgol(smoothed, window_length=15, polyorder=3, deriv=1)
d2 = derivative_savgol(smoothed, window_length=15, polyorder=3, deriv=2)

fig, axes = plt.subplots(3, 1, figsize=(10, 8), sharex=True)

axes[0].plot(wavenumbers, smoothed)
axes[0].set_ylabel("Intensity")
axes[0].set_title("Smoothed Spectrum")

axes[1].plot(wavenumbers, d1)
axes[1].axhline(0, color="gray", linestyle="--", linewidth=0.5)
axes[1].set_ylabel("1st Derivative")
axes[1].set_title("First Derivative (zero crossings = peak positions)")

axes[2].plot(wavenumbers, d2)
axes[2].axhline(0, color="gray", linestyle="--", linewidth=0.5)
axes[2].set_ylabel("2nd Derivative")
axes[2].set_title("Second Derivative (minima = peak positions)")
axes[2].set_xlabel("Wavenumber")

for ax in axes:
    ax.invert_xaxis()

plt.tight_layout()
plt.show()

# Smooth first to reduce noise amplification smoothed = smooth_savgol(noisy, window_length=15, polyorder=3) # First and second derivatives d1 = derivative_savgol(smoothed, window_length=15, polyorder=3, deriv=1) d2 = derivative_savgol(smoothed, window_length=15, polyorder=3, deriv=2) fig, axes = plt.subplots(3, 1, figsize=(10, 8), sharex=True) axes[0].plot(wavenumbers, smoothed) axes[0].set_ylabel("Intensity") axes[0].set_title("Smoothed Spectrum") axes[1].plot(wavenumbers, d1) axes[1].axhline(0, color="gray", linestyle="--", linewidth=0.5) axes[1].set_ylabel("1st Derivative") axes[1].set_title("First Derivative (zero crossings = peak positions)") axes[2].plot(wavenumbers, d2) axes[2].axhline(0, color="gray", linestyle="--", linewidth=0.5) axes[2].set_ylabel("2nd Derivative") axes[2].set_title("Second Derivative (minima = peak positions)") axes[2].set_xlabel("Wavenumber") for ax in axes: ax.invert_xaxis() plt.tight_layout() plt.show()

Gap-Segment Derivative Comparison¶

In [ ]:

d1_sg = derivative_savgol(smoothed, window_length=15, polyorder=3, deriv=1)
d1_gap = derivative_gap_segment(smoothed, gap=7, segment=7, deriv=1)

fig, axes = plt.subplots(1, 2, figsize=(14, 4), sharex=True)

axes[0].plot(wavenumbers, d1_sg)
axes[0].set_title("Savitzky-Golay 1st Derivative")
axes[0].axhline(0, color="gray", linestyle="--", linewidth=0.5)

axes[1].plot(wavenumbers, d1_gap)
axes[1].set_title("Gap-Segment 1st Derivative")
axes[1].axhline(0, color="gray", linestyle="--", linewidth=0.5)

for ax in axes:
    ax.set_xlabel("Wavenumber")
    ax.invert_xaxis()

plt.tight_layout()
plt.show()

d1_sg = derivative_savgol(smoothed, window_length=15, polyorder=3, deriv=1) d1_gap = derivative_gap_segment(smoothed, gap=7, segment=7, deriv=1) fig, axes = plt.subplots(1, 2, figsize=(14, 4), sharex=True) axes[0].plot(wavenumbers, d1_sg) axes[0].set_title("Savitzky-Golay 1st Derivative") axes[0].axhline(0, color="gray", linestyle="--", linewidth=0.5) axes[1].plot(wavenumbers, d1_gap) axes[1].set_title("Gap-Segment 1st Derivative") axes[1].axhline(0, color="gray", linestyle="--", linewidth=0.5) for ax in axes: ax.set_xlabel("Wavenumber") ax.invert_xaxis() plt.tight_layout() plt.show()

Peak Finding¶

In [ ]:

result = peaks_find(smoothed, wavenumbers, prominence=0.3, distance=20)

plt.figure(figsize=(10, 4))
plt.plot(wavenumbers, smoothed)
plt.plot(result.wavenumbers, result.heights, "rv", markersize=10, label="Peaks")

for wn, h in zip(result.wavenumbers, result.heights, strict=True):
    plt.annotate(
        f"{wn:.0f}", (wn, h), textcoords="offset points", xytext=(0, 10), ha="center", fontsize=9
    )

plt.xlabel("Wavenumber")
plt.ylabel("Intensity")
plt.title(f"Found {len(result.indices)} Peaks")
plt.legend()
plt.gca().invert_xaxis()
plt.show()

print(f"Peak positions: {result.wavenumbers}")
print(f"Peak heights: {result.heights}")

result = peaks_find(smoothed, wavenumbers, prominence=0.3, distance=20) plt.figure(figsize=(10, 4)) plt.plot(wavenumbers, smoothed) plt.plot(result.wavenumbers, result.heights, "rv", markersize=10, label="Peaks") for wn, h in zip(result.wavenumbers, result.heights, strict=True): plt.annotate( f"{wn:.0f}", (wn, h), textcoords="offset points", xytext=(0, 10), ha="center", fontsize=9 ) plt.xlabel("Wavenumber") plt.ylabel("Intensity") plt.title(f"Found {len(result.indices)} Peaks") plt.legend() plt.gca().invert_xaxis() plt.show() print(f"Peak positions: {result.wavenumbers}") print(f"Peak heights: {result.heights}")

Peak Integration¶

In [ ]:

# Define integration ranges around each peak
ranges = [(1500, 1750), (2800, 3050)]
areas = peaks_integrate(smoothed, wavenumbers, ranges=ranges)

plt.figure(figsize=(10, 4))
plt.plot(wavenumbers, smoothed)

colors = ["skyblue", "salmon"]
for (lo, hi), area, color in zip(ranges, areas, colors, strict=True):
    mask = (wavenumbers >= lo) & (wavenumbers <= hi)
    plt.fill_between(
        wavenumbers[mask],
        smoothed[mask],
        alpha=0.3,
        color=color,
        label=f"{lo}-{hi}: area={area:.2f}",
    )

plt.xlabel("Wavenumber")
plt.ylabel("Intensity")
plt.title("Peak Integration")
plt.legend()
plt.gca().invert_xaxis()
plt.show()

# Define integration ranges around each peak ranges = [(1500, 1750), (2800, 3050)] areas = peaks_integrate(smoothed, wavenumbers, ranges=ranges) plt.figure(figsize=(10, 4)) plt.plot(wavenumbers, smoothed) colors = ["skyblue", "salmon"] for (lo, hi), area, color in zip(ranges, areas, colors, strict=True): mask = (wavenumbers >= lo) & (wavenumbers <= hi) plt.fill_between( wavenumbers[mask], smoothed[mask], alpha=0.3, color=color, label=f"{lo}-{hi}: area={area:.2f}", ) plt.xlabel("Wavenumber") plt.ylabel("Intensity") plt.title("Peak Integration") plt.legend() plt.gca().invert_xaxis() plt.show()