Peak Analysis Workflow
This tutorial shows an end-to-end peak-analysis workflow in SpectroChemPy:
prepare a spectrum,
detect peaks with
find_peaks(),inspect the result through
PeakFindingResultandPeakTable,export the table,
write and validate a fitting script,
fit the spectrum and inspect the resulting
FitResult.
It complements the dedicated tutorials on peak finding and fitting by focusing on the bridge between detection and fitting.
[1]:
from pathlib import Path
from tempfile import TemporaryDirectory
import spectrochempy as scp
Load and prepare a spectrum
We build a small synthetic spectrum with two broad OH-like peaks on top of a gentle baseline and a little Gaussian noise. This keeps the tutorial self-contained while still showing the full workflow on a slightly more realistic profile.
[2]:
prefs = scp.preferences
prefs.figure.figsize = (7, 3)
x = scp.Coord.linspace(3700.0, 3300.0, 1200, title="wavenumber", units="cm^-1")
baseline = scp.polynomial(
x,
offset=0.015,
slope=0.00002,
ampl=1.0,
c_2=1.5e-6,
)
peak_1 = scp.gaussian(x, ampl=0.95, pos=3624.0, width=42.39, normalized=False)
peak_2 = scp.gaussian(x, ampl=0.32, pos=3542.0, width=51.81, normalized=False)
noise = scp.normal(loc=0.0, scale=0.007, size=x.size)
nd_oh = baseline + peak_1 + peak_2 + noise
nd_oh.title = "Synthetic OH region"
nd_oh_corr = scp.basc(
nd_oh,
[3700.0, 3670.0],
[3490.0, 3300.0],
model="polynomial",
order=2,
)
ax = nd_oh.plot(label="Synthetic spectrum")
_ = nd_oh_corr.plot(clear=False, label="Baseline-corrected spectrum")
ax.legend()
[2]:
<matplotlib.legend.Legend at 0x7f687ad37770>
Detect peaks and inspect the structured result
find_peaks(..., as_result=True) returns a PeakFindingResult instead of the historical (peaks, properties) tuple. The result keeps the detected peak dataset and exposes a stable tabular view through result.table.
[3]:
result = nd_oh_corr.find_peaks(
height=0.05,
distance="20 cm^-1",
prominence=0.1,
width="10 cm^-1",
as_result=True,
)
result
[3]:
PeakFindingResult(n_peaks=2)
[4]:
table = result.table
table
[4]:
PeakTable(n_peaks=2)
PeakTable gives us a dependency-light view of the detected peaks. The raw peak dataset and the raw SciPy-style property dictionary are still available on result, but the table is often a better starting point for inspection, export, and later workflow steps.
[5]:
rows = table.to_dict()
rows[:4]
[5]:
[{'index': 0,
'position': <Quantity(3624.693, '1 / centimeter')>,
'height': 0.9517825919901952,
'peak_height': 0.9512735027031923,
'prominence': 0.9725659664101931,
'left_base': <Quantity(3693.995, '1 / centimeter')>,
'right_base': <Quantity(3421.435, '1 / centimeter')>,
'width': <Quantity(42.9368656, '1 / centimeter')>,
'width_height': 0.4649905194980958,
'left_ip': <Quantity(3645.48428, '1 / centimeter')>,
'right_ip': <Quantity(3602.54564, '1 / centimeter')>},
{'index': 1,
'position': <Quantity(3541.506, '1 / centimeter')>,
'height': 0.32834620657376945,
'peak_height': 0.3282738996567547,
'prominence': 0.22395564656545003,
'left_base': <Quantity(3579.9, '1 / centimeter')>,
'right_base': <Quantity(3421.435, '1 / centimeter')>,
'width': <Quantity(36.9047718, '1 / centimeter')>,
'width_height': 0.21629607637402967,
'left_ip': <Quantity(3561.18478, '1 / centimeter')>,
'right_ip': <Quantity(3524.27849, '1 / centimeter')>}]
We can also look at the available table columns:
[6]:
table.columns
[6]:
('index',
'position',
'height',
'peak_height',
'prominence',
'left_base',
'right_base',
'width',
'width_height',
'left_ip',
'right_ip')
Export the peak table
PeakTable.to_csv() writes a simple CSV file without adding any optional dependency such as pandas.
[7]:
with TemporaryDirectory() as tmpdir:
csv_path = Path(tmpdir) / "nh4y-oh-peaks.csv"
_ = table.to_csv(csv_path)
preview = "\n".join(csv_path.read_text(encoding="utf-8").splitlines()[:4])
print(preview)
index,position,height,peak_height,prominence,left_base,right_base,width,width_height,left_ip,right_ip
0,3624.693 cm⁻¹,0.9517825919901952,0.9512735027031923,0.9725659664101931,3693.995 cm⁻¹,3421.435 cm⁻¹,42.936865629178655 cm⁻¹,0.4649905194980958,3645.4842806647957 cm⁻¹,3602.5456396140867 cm⁻¹
1,3541.506 cm⁻¹,0.32834620657376945,0.3282738996567547,0.22395564656545003,3579.9 cm⁻¹,3421.435 cm⁻¹,36.904771783932375 cm⁻¹,0.21629607637402967,3561.184783028414 cm⁻¹,3524.278486185313 cm⁻¹
Select starting candidates for fitting
Peak detection gives geometric candidates. Fitting still requires a modeling decision: which peaks do we want to fit, with which line shape, and with which bounds? Here we keep the two strongest detected peaks and use their positions as initial guesses in a manually written fitting script.
[8]:
selected_table = table.top(2, by="height").sort_by(
"position",
reverse=True,
unit="cm^-1",
)
positions = selected_table.column("position", unit="cm^-1", as_float=True)
heights = selected_table.column("height", as_float=True)
widths = selected_table.column("width", unit="cm^-1", as_float=True)
for position, height, width in zip(positions, heights, widths, strict=False):
print(
f"candidate peak at {position:.2f} cm^-1 "
f"with height {height:.3f} and width {width:.2f} cm^-1"
)
candidate peak at 3624.69 cm^-1 with height 0.952 and width 42.94 cm^-1
candidate peak at 3541.51 cm^-1 with height 0.328 and width 36.90 cm^-1
[9]:
script = f"""
COMMON:
$ gratio: 0.1, 0.0, 1.0
$ gasym: 0.1, 0.0, 1.0
MODEL: LINE_1
shape: asymmetricvoigtmodel
$ ampl: 1.0, 0.0, none
$ pos: {positions[0]:.2f}, {positions[0] - widths[0]:.2f}, {positions[0] + widths[0]:.2f}
> ratio: gratio
> asym: gasym
$ width: {widths[0]:.2f}, {0.5 * widths[0]:.2f}, {2.0 * widths[0]:.2f}
MODEL: LINE_2
shape: asymmetricvoigtmodel
$ ampl: 0.2, 0.0, none
$ pos: {positions[1]:.2f}, {positions[1] - widths[1]:.2f}, {positions[1] + widths[1]:.2f}
> ratio: gratio
> asym: gasym
$ width: {widths[1]:.2f}, {0.5 * widths[1]:.2f}, {2.0 * widths[1]:.2f}
"""
print(script)
COMMON:
$ gratio: 0.1, 0.0, 1.0
$ gasym: 0.1, 0.0, 1.0
MODEL: LINE_1
shape: asymmetricvoigtmodel
$ ampl: 1.0, 0.0, none
$ pos: 3624.69, 3581.76, 3667.63
> ratio: gratio
> asym: gasym
$ width: 42.94, 21.47, 85.87
MODEL: LINE_2
shape: asymmetricvoigtmodel
$ ampl: 0.2, 0.0, none
$ pos: 3541.51, 3504.60, 3578.41
> ratio: gratio
> asym: gasym
$ width: 36.90, 18.45, 73.81
Validate the script before fitting
Optimize.validate_script() lets us check the script before launching the optimization.
[10]:
opt = scp.Optimize(log_level="INFO")
errors = opt.validate_script(script)
errors
[10]:
[]
An empty list means that the script is structurally valid and all referenced models are recognized.
Fit the spectrum and inspect the result
[11]:
opt.script = script
opt.max_iter = 20000
_ = opt.fit(nd_oh_corr)
**************************************************
Result:
**************************************************
COMMON:
$ gratio: 1.0000, 0.0, 1.0
$ gasym: 0.0156, 0.0, 1.0
MODEL: line_1
shape: asymmetricvoigtmodel
$ ampl: 0.9428, 0.0, none
> asym:gasym
$ pos: 3623.9841, 3581.76, 3667.63
> ratio:gratio
$ width: 42.2079, 21.47, 85.87
MODEL: line_2
shape: asymmetricvoigtmodel
$ ampl: 0.3163, 0.0, none
> asym:gasym
$ pos: 3541.9529, 3504.6, 3578.41
> ratio:gratio
$ width: 50.8726, 18.45, 73.81
[12]:
fit_result = opt.result
fit_result
[12]:
FitResult — Optimize
Parameters (5)
Outputs (3)
NDDataset [polynomial_Optimize.inverse_transform] — float64, shape: (u:1, x:1200)
Data
Dimension `u`
Dimension `x`
NDDataset [polynomial_Optimize.components] — float64, shape: (k:2, x:1200)
Data
[2.018e-13 2.278e-13 ... 1.689e-26 1.447e-26]]
Dimension `k`
Dimension `x`
NDDataset [polynomial] — float64, shape: (u:1, x:1200)
2026-07-07 17:32:34+00:00> Binary operation add with `gaussian` has been performed
2026-07-07 17:32:34+00:00> Binary operation add with `gaussian` has been performed
2026-07-07 17:32:34+00:00> Binary operation add with `NDDataset_d6cc5656` has been performed
2026-07-07 17:32:35+00:00> Binary operation sub with `polynomial_Baseline.baseline` has been performed
2026-07-07 17:32:36+00:00> Binary operation sub with `polynomial_Optimize.inverse_transform` has been performed
Data
Dimension `u`
Dimension `x`
Diagnostics (15)
FitResult groups fitted outputs and diagnostics. The existing estimator surface (opt.predict(), opt.components, plotting helpers, and so on) remains available, but opt.result is the stable result object for inspection.
[13]:
components = fit_result.components
{
"r_squared": fit_result.diagnostics["r_squared"],
"adjusted_r_squared": fit_result.diagnostics["adjusted_r_squared"],
"rmse": fit_result.diagnostics["rmse"],
"degrees_of_freedom": fit_result.diagnostics["degrees_of_freedom"],
"reduced_chi_square": fit_result.diagnostics["reduced_chi_square"],
"success": fit_result.diagnostics["success"],
}
[13]:
{'r_squared': 0.9991175009565901,
'adjusted_r_squared': 0.9991123184957648,
'rmse': 0.007220378188779769,
'degrees_of_freedom': 1192,
'reduced_chi_square': 5.248375287483887e-05,
'success': True}
The fitted components remain regular datasets, so they can be plotted directly against the corrected spectrum.
[14]:
_ = nd_oh_corr.plot()
ax = components[:].plot(clear=False)
ax.autoscale(enable=True, axis="y")
[15]:
_ = opt.plotmerit(offset=0, kind="scatter")
Summary
This workflow now has a clear progression:
find_peaks()detects candidates,PeakFindingResultstores the structured detection output,PeakTableprovides a stable tabular view for inspection and export,Optimize.validate_script()checks the fitting DSL before optimization,FitResultgroups fitted outputs and diagnostics.
This is the current reference workflow when moving from peak detection to peak fitting in SpectroChemPy.