Tutorials

Gaussian Peak Fitting

When a mass spectrometer, fluorescence detector, or particle counter produces a peak, you want more than "there is a peak somewhere around 450 nm" — you need the exact centre (wavelength of maximum intensity), the full width at half maximum (FWHM, which relates to resolution), and the area (proportional to concentration). Simply reading the maximum of noisy data gives you an imprecise estimate; fitting a Gaussian model gives you a statistically optimal estimate *and* uncertainty bounds for each parameter. `scipy.optimize.curve_fit` solves this in a few lines: define the model function, provide initial guesses, and get back the fitted parameters plus their covariance matrix.

### Simulating a noisy Gaussian peak

A single Gaussian peak with known parameters provides a ground truth to verify that the fitting recovers the correct values.


import numpy as np

rng = np.random.default_rng(42)
x = np.linspace(-5, 5, 200)

# True parameters
true_amp, true_center, true_sigma = 4.0, 0.5, 0.8

signal_clean = true_amp * np.exp(-(x - true_center)**2 / (2 * true_sigma**2))
signal_noisy = signal_clean + rng.normal(0, 0.25, len(x))

fwhm = 2 * np.sqrt(2 * np.log(2)) * true_sigma
print(f"True amplitude:  {true_amp}")
print(f"True centre:     {true_center}")
print(f"True sigma:      {true_sigma:.2f}  → FWHM = {fwhm:.3f}")
print(f"Noise std:       0.25  (SNR ≈ {true_amp / 0.25:.1f}x)")
True amplitude:  4.0
True centre:     0.5
True sigma:      0.80  → FWHM = 1.884
Noise std:       0.25  (SNR ≈ 16.0x)

- The Gaussian model `A * exp(-(x-μ)²/(2σ²))` has three parameters: amplitude `A`, centre `μ`, and width `σ`.
- FWHM (full width at half maximum) = `2√(2 ln 2) × σ ≈ 2.355 × σ` is the standard width measure in spectroscopy and chromatography.
- An SNR of 16:1 represents a typical detector with moderate noise — low enough that reading the peak position from the raw data directly would introduce noticeable error.

### Fitting the Gaussian model with curve_fit

`curve_fit` minimises the sum of squared residuals between the model and data, returning the best-fit parameters and their covariance matrix.


import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

rng = np.random.default_rng(42)
x = np.linspace(-5, 5, 200)
true_amp, true_center, true_sigma = 4.0, 0.5, 0.8
signal_clean = true_amp * np.exp(-(x - true_center)**2 / (2 * true_sigma**2))
signal_noisy = signal_clean + rng.normal(0, 0.25, len(x))

def gaussian(x, amplitude, center, sigma):
    return amplitude * np.exp(-(x - center)**2 / (2 * sigma**2))

p0 = [1.0, 0.0, 1.0]  # initial guesses: amp, center, sigma
popt, pcov = curve_fit(gaussian, x, signal_noisy, p0=p0)
perr = np.sqrt(np.diag(pcov))  # 1-sigma parameter uncertainties

print(f"{'Parameter':>12}  {'True':>8}  {'Fitted':>8}  {'±1σ':>8}")
print("-" * 44)
for name, true, fit, err in zip(
        ["Amplitude", "Centre", "Sigma"],
        [true_amp, true_center, true_sigma],
        popt, perr):
    print(f"{name:>12}  {true:>8.3f}  {fit:>8.3f}  {err:>8.4f}")

fitted_fwhm = 2 * np.sqrt(2 * np.log(2)) * popt[2]
print(f"\nFitted FWHM: {fitted_fwhm:.3f}  (true: {2*np.sqrt(2*np.log(2))*true_sigma:.3f})")

fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(x, signal_noisy,  color="lightgray",  linewidth=0.8, label="Noisy data")
ax.plot(x, signal_clean,  color="black",      linewidth=1.0, linestyle="--", label="True signal")
ax.plot(x, gaussian(x, *popt), color="tomato", linewidth=2, label="Gaussian fit")
ax.set_xlabel("x"); ax.set_ylabel("Intensity")
ax.set_title("Gaussian peak fitting")
ax.legend()
plt.tight_layout()
plt.show()
   Parameter      True    Fitted       ±1σ
--------------------------------------------
   Amplitude     4.000     3.954    0.0508
      Centre     0.500     0.514    0.0118
       Sigma     0.800     0.795    0.0118

Fitted FWHM: 1.872  (true: 1.884)

- `p0` provides initial guesses for [amplitude, centre, sigma]. Good guesses — approximately the right magnitude — help the optimiser converge. Bad guesses can lead to wrong local minima.
- `pcov` is the covariance matrix of the fitted parameters; `np.sqrt(np.diag(pcov))` gives the 1-sigma standard errors. Multiply by 1.96 for 95% confidence intervals.
- The fitted parameters should be close to the true values, with uncertainties that reflect how well the noisy data constrains each parameter.

### Fitting overlapping peaks

When two peaks overlap, fitting each independently gives wrong results. Fitting the sum of two Gaussians simultaneously recovers both peaks correctly.


import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

rng = np.random.default_rng(42)
x = np.linspace(-3, 7, 300)

def two_gaussians(x, a1, c1, s1, a2, c2, s2):
    g1 = a1 * np.exp(-(x - c1)**2 / (2 * s1**2))
    g2 = a2 * np.exp(-(x - c2)**2 / (2 * s2**2))
    return g1 + g2

# True: peak 1 at x=0, peak 2 at x=2 (overlapping)
y_true = two_gaussians(x, 3.0, 0.0, 0.6, 2.0, 2.0, 0.7)
y_noisy = y_true + rng.normal(0, 0.2, len(x))

p0 = [2.5, -0.5, 0.5, 1.5, 2.5, 0.5]
popt, _ = curve_fit(two_gaussians, x, y_noisy, p0=p0)

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

fig, ax = plt.subplots(figsize=(9, 4))
ax.plot(x, y_noisy, color="lightgray", linewidth=0.8, label="Data")
ax.plot(x, two_gaussians(x, *popt), color="tomato", linewidth=2, label="Combined fit")
ax.plot(x, gaussian(x, *popt[:3]),  color="steelblue", linewidth=1.5, linestyle="--",
        label=f"Peak 1: c={popt[1]:.2f}, σ={popt[2]:.2f}")
ax.plot(x, gaussian(x, *popt[3:]),  color="green",     linewidth=1.5, linestyle="--",
        label=f"Peak 2: c={popt[4]:.2f}, σ={popt[5]:.2f}")
ax.set_xlabel("x"); ax.set_ylabel("Intensity")
ax.set_title("Two overlapping Gaussian peaks")
ax.legend(fontsize=9)
plt.tight_layout()
plt.show()
- `two_gaussians` takes 6 parameters (3 per peak); `curve_fit` fits all simultaneously.
- For heavily overlapping peaks, initial guesses are critical — start with your best visual estimate of each peak's location and height.
- The dashed curves show the deconvolved individual peaks, which could not be resolved by eye from the overlapping raw data.

### Conclusion

Gaussian fitting is the standard approach for peak analysis in spectroscopy, chromatography, and any measurement where features have a bell-shaped profile. Always check the residuals (data minus fit) — systematic patterns in the residuals mean the Gaussian assumption doesn't hold and a different model may be needed. For peaks with heavier tails than a Gaussian predicts, see [Lorentzian peak fitting](/tutorials/lorentzian-peak-fitting).

For detecting peak positions without fitting a model, see [finding peaks with scipy](/tutorials/find-peaks-with-scipy). For fitting a Lorentzian profile instead, see [Lorentzian peak fitting](/tutorials/lorentzian-peak-fitting).