fepydas

fepydas ("felix' python data analysis suite") is a Python 3 library for analysing and fitting optical-spectroscopy and photophysics data: photoluminescence (PL), PL excitation (PLE), Raman, time-resolved PL (TCSPC decays), external quantum efficiency (EQE) and spatially-resolved spectrum maps.

It gives you a small, unit-aware data model plus ready-made readers, fitters and plotters, so a typical analysis is just: load a file → manipulate the dataset → fit → plot/export.

from fepydas.constructors.datareaders.Custom import ASCII_Spectrum
from fepydas.workers.Fit import GaussianFit

spectrum = ASCII_Spectrum("measurement.txt")   # -> a Spectrum
spectrum.cutAxis(500, 600)                      # keep 500-600 nm
spectrum.normalize()                            # scale peak to 1

fit = GaussianFit()
fit.initializeAutoFromSpectrum(spectrum)        # guess starting parameters
fit.fitSpectrum(spectrum)                        # run the fit

spectrum.toPlot().save("spectrum.png")          # save a plot
spectrum.export("spectrum.txt")                 # save tab-separated data

• Install: pip install fepydas (from PyPI)
• Repository: https://git.tu-berlin.de/nippert/fepydas
• Runnable examples: examples/

Table of contents

1. Installation
2. Core concepts
3. Loading data
4. Working with spectra and series
5. Units, axes and calibration
6. Fitting
7. Time-resolved decays (IRF convolution)
8. Spectrum maps
9. PLE measurements
10. Decomposition and clustering
11. Plotting
12. Exporting and saving
13. Batch processing whole directories
14. Batch TRPL lifetime fitting
15. Worked examples
16. API map
17. Requirements

Installation

pip install fepydas

Python 3.9+ is required (since numpy >= 2.0 is a dependency). All other dependencies (scipy, matplotlib, lmfit, scikit-learn, pyamg, sympy, pathos, tqdm, ptufile) are installed automatically. See Requirements.

To work from a checkout instead:

git clone https://git.tu-berlin.de/nippert/fepydas
cd fepydas
pip install -e .

Core concepts

fepydas is organised in three layers. You usually only touch the top two.

The data model

• Data / Data1D / Data2D / Data3D — a numpy array (.values) plus a DataType (.datatype) and optional .errors. The dimensionality matches the structure: a single spectrum is Data1D, a stack of spectra is Data2D, a spatial cube is Data3D.
• DataType — a quantity name + a Unit (e.g. "Wavelength" in nm). Predefined ones live in fepydas.datatypes (WL_nm, E_eV, COUNTS, TIME_ns, TEMP_K, ...).
• Unit — wraps a sympy unit so values can be compared and converted (including the special wavelength ↔ energy case, with the Jacobian applied to intensities).

Datasets

A Dataset pairs an x-axis with data. The hierarchy specialises this:

Dataset
├── Spectrum                       one spectrum (axis + Data1D)
└── SpectrumSeries                 many spectra sharing an axis, indexed by a "key" axis
    ├── BinnedSpectrumSeries       a series rebinned along its key axis
    ├── JoinSpectrumSeries         several separate spectra stacked into a series
    ├── MergeSpectraSeries         several series concatenated
    └── SpectrumSeriesWithCalibration
        ├── PLE                    emission spectra keyed by excitation wavelength
        └── BaseSpectrumMap
            ├── SpectrumMap        dense rectangular map (spectrum at every pixel)
            └── SparseSpectrumMap  map where only some pixels were measured

Spectrum
└── SpectrumWithIRF                a decay histogram bundled with its instrument response
    └── SpectraWithCommonIRF       several decays sharing one IRF

Map                                a single scalar value per map pixel (e.g. a fit result)

Most dataset methods mutate in place and return None (e.g. cutAxis, normalize, convertToEnergy). Methods that derive a new object return it (e.g. integrate, average, collapse, toPlot).

The typical workflow

ds = SomeReader(path)            # constructors.datareaders.* -> a dataset
ds.cutAxis(...); ds.normalize()  # in-place dataset operations
fit = SomeFit(); fit.initializeAutoFromSpectrum(ds); fit.fitSpectrum(ds)
ds.toPlot().save("plot.png")     # constructors.Plots
ds.export("data.txt")            # workers.ASCIIWriter

Loading data

All readers are plain functions that return a dataset. Pick the one that matches your file.

Generic / in-memory (fepydas.constructors.datareaders.Generic)

from fepydas.constructors.datareaders.Generic import (
    GenericXY, GenericXMultiY, ASCII, Binary,
)
import numpy

# Build a Spectrum from arrays you already have. Units are abbreviation strings.
x = numpy.linspace(400, 700, 300)
y = numpy.exp(-((x - 550) ** 2) / 50)
spectrum = GenericXY(x, "Wavelength", "nm", y, "Counts", "")

# Build a SpectrumSeries from a shared axis and a 2D stack of y-data.
series = GenericXMultiY(x, "Wavelength", "nm", y2d, "Counts", "",
                        keys, "Temperature", "K")

# Read a raw numeric table (thin numpy.genfromtxt wrapper).
table = ASCII("data.txt", delimiter="\t", skip_header=0)

# Reload anything previously saved with .saveBinary()
obj = Binary("dataset.bin")

TCSPC histograms (Generic + Proprietary)

from fepydas.constructors.datareaders.Generic import Histogram, HistogramWithIRF

decay = Histogram("decay.phu")                      # .phu/.thi/.sdt/.ptu -> Spectrum
pair  = HistogramWithIRF("decay.sdt", "irf.sdt")    # -> SpectrumWithIRF

Histogram auto-detects the format from the extension: .phu (PicoHarp), .thi (TimeHarp), .sdt (Becker & Hickl), .ptu (PicoQuant TTTR). For large .ptu files the reader is chunked/low-RAM and takes options (forwarded through **kwargs): channel, bin_width_s, all_channels, number_of_bins, use_memmap, show_progress, progress_callback, ... See the Histogram/PTUHistogram docstrings.

Lab-specific ASCII formats (fepydas.constructors.datareaders.Custom)

LabVIEW exports (fepydas.constructors.datareaders.FromLabview)

from fepydas.constructors.datareaders.FromLabview import Automatic

ds = Automatic("measurement.plm")   # inspects the header and returns the right type

Automatic detects the measurement type from the file header and returns a Spectrum, SpectrumSeries, PLE, SparseSpectrumMap, Map, or a list[Spectrum] for PLE-response files. Dark-frame correction is applied automatically when present; tracked power is attached as dataset.power/infos["Power Track"] where supported so you can apply a power correction explicitly.

Working with spectra and series

These operations are available on Spectrum and (where it makes sense) SpectrumSeries. They mutate the dataset in place unless noted.

# --- cropping ---
spectrum.cutAxis(500, 600)          # keep an axis value range
spectrum.cutAxisByIdx(50, 3000)     # keep an axis index range
series.cutKeys(10, 300)             # keep a key range (series only)

# --- intensity operations ---
spectrum.normalize()                # scale max to 1
series.normalize(individual=True)   # scale each spectrum to its own max
series.subtractBaseline()           # remove a flat baseline (no dark spectrum available)
spectrum.antiBloom(45000, width=2)  # mask CCD blooming above a threshold
spectrum.divideBySpectrum(reference)

# --- reductions (these RETURN new datasets) ---
integral = series.integrate()       # integrated intensity vs key  -> Spectrum
avg      = series.average()         # mean spectrum                -> Spectrum
total    = series.collapse()        # summed spectrum              -> Spectrum
peakpos  = series.maximum()         # peak axis position vs key    -> Spectrum

# --- selecting / combining ---
series.filterByIntegral(0.5)        # drop weak spectra (< 50% of the strongest integral)
hdr = series.highDynamicRange()     # merge differently-exposed spectra -> Spectrum

identifyPeaks(threshold=10, width=10) returns index ranges around detected peaks (used by the calibration helpers).

Units, axes and calibration

Predefined data types

Import from fepydas.datatypes: WL_nm, RS_cm (Raman shift), E_eV, COUNTS, NUMBER, TIME_s, TIME_ns, TIME_ps, TEMP_K, ANGLE_au, I_mA, P_mW, POS_um, POS_mm. Build your own with generateDataType("Power", "mW") (the unit string is parsed by unitFromAbbrev).

Wavelength ↔ energy

spectrum.convertToEnergy("eV")   # converts the axis AND applies the Jacobian to intensities

This works on Spectrum, SpectrumSeries and the maps. At the Data1D level you can also call axis.transformTo(E_eV) for any compatible unit.

Vacuum correction

spectrum.applyVacuumCorrection(temperature=20, pressure=101325, humidity=0, co2=610)

Converts air wavelengths to vacuum wavelengths using the refractive index of moist air.

Spectral calibration

If a calibration spectrum (e.g. a lamp with known lines) is attached (PLE, the map types, or anything that is a SpectrumSeriesWithCalibration):

ds.calibrate(references=[None, 435.8335, 404.6565])  # known positions for detected peaks

Standalone, you can fit a calibration curve from measured vs reference line positions:

from fepydas.workers.Fit import Calibration

cal = Calibration(measured_positions, nist_positions)   # fits on construction; saves a plot
calibrated_axis = cal.apply(other_axis.values)

Fitting

Every fit model is a Fit subclass. The pattern is always: create → initialise parameters → fit → use the result.

from fepydas.workers.Fit import GaussianFit

fit = GaussianFit()
fit.initializeAutoFromSpectrum(spectrum)   # guess bg, I, x0, fwhm from the data
# (optional) constrain parameters before fitting:
fit.parameters["x0"].min = 419
fit.parameters["x0"].max = 421
fit.fitSpectrum(spectrum)                   # or fit.fit(x, y)

print(fit.result.params)                    # lmfit result
model_curve = fit.evaluate(spectrum.axis.values)

Available models

You can set parameters manually instead of initializeAuto*:

from fepydas.workers.Fit import Fit
from lmfit import Parameters

def model(x, a=1, b=0):
    return a * x + b

fit = Fit(model)
p = Parameters()
p.add("a", value=1, min=0)
p.add("b", value=0)
fit.initializeParameters(p)
fit.fitSpectrum(spectrum)

Fitting a whole series in parallel

fit = GaussianFit()
fit.initializeAutoFromSpectrum(series.average())
series.fit(fit)                               # fits every spectrum (uses a process pool)

Windows note: parallel fitting forks a process pool. On Windows, guard your script with if __name__ == "__main__": so child processes don't re-run the whole module.

Saving fit reports

fit.saveReport("fit.txt")          # parameter table + statistics
fit.saveLMfitReport("lmfit.txt")   # the raw lmfit fit_report

FitSimple — fitting with native lmfit models

For workflows built directly on lmfit.Model, FitSimple offers parameter guessing, optional parallel fitting and dynamic bound propagation:

from fepydas.workers.Fit import FitSimple
from lmfit.models import GaussianModel

fs = FitSimple(GaussianModel(), spectrum_or_series)
fs.initguess()                         # or fs.setinitparams(params)
fs.fit(multiprocessing=True)           # y_err=..., bounds_tolerance=..., cores=... all supported

Time-resolved decays (IRF convolution)

For TCSPC lifetime fitting, load the decay together with its instrument response function and fit a model convolved with the IRF:

from fepydas.constructors.datareaders.Generic import HistogramWithIRF
from fepydas.constructors.Plots import SpectrumWithIRFPlot
from fepydas.workers.Fit import ExponentialFit

s = HistogramWithIRF("Signal.sdt", "InstrumentResponseFunction.sdt")
s.normalize()

fit = ExponentialFit()
fit.initializeAutoFromSpectrum(s)
fit.parameters["rise"].value = 0.01
fit.parameters["rise"].vary = False
fit.convolutedFit(s.axis.values, s.data.values, s.IRF.values)

plot = SpectrumWithIRFPlot(s)
plot.addLine(fit.evaluateConvolution(s.axis.values, s.IRF.values), label="Convolution Fit")
plot.setYLog(); plot.legend()
plot.save("ConvolutionFit.png")

Spectrum maps

Maps store a spectrum at every spatial pixel. SpectrumMap is a dense rectangular grid; SparseSpectrumMap (e.g. from a masked LabVIEW scan) stores only measured pixels.

from fepydas.constructors.datareaders.Custom import Mapscan
from fepydas.workers.Fit import GaussianFit
from fepydas.datatypes.Data import DataType

if __name__ == "__main__":           # required for parallel fitting on Windows
    m = Mapscan("Example.map")
    m.cutAxis(290, 360)              # restrict the spectral window

    m.integrate().toPlot().save("integrated.png")   # integrated intensity image
    m.maximum().toPlot().save("peakpos.png")        # peak-position image
    m.average().toPlot().save("avg_spectrum.png")   # mean spectrum

    # fit every pixel and turn a fit parameter into a map
    fit = GaussianFit()
    fit.initializeAutoFromSpectrum(m.average())
    m.fit(fit)
    posMap = m.fitParameterToMap(fit, "x0", DataType("Position", "nm"))
    posMap.toPlot(zRange=[330, 350]).save("fitted_position.png")

Maps also support clustering-based segmentation:

clusters = m.fitMixture(n=2, nmf=3)               # spectral clustering (optionally NMF-reduced)
m.getClusterMapIdx(clusters).toPlot().save("clusters.png")
m.getClusterIntegratedSpectra(clusters).toPlot().save("cluster_spectra.png")

integrate(), maximum() and fitParameterToMap(...) return Map objects (one value per pixel), which plot as false-color images via toPlot(log=..., zRange=...).

PLE measurements

PLE is a series of emission spectra keyed by excitation wavelength, with extra correction steps. A full PLE workflow (from the PLE example):

from fepydas.constructors.datareaders.Custom import ASCII_PLE

ple = ASCII_PLE("PLE.txt", "Calibration.cal", "IntensityRef.dat")
ple.toPlot().save("01_raw.png")

ple.calibrate([None, 435.8335, 404.6565])     # calibrate the detection axis
ple.applyVacuumCorrection()
ple.subtractBaseline()
ple.applyLampResponseCorrection()              # divide by the lamp response
ple.calibrateLamp()                            # calibrate the excitation axis from the lamp line
ple.export("CalibratedPLE.txt")

ple.integrate().toPlot().save("integrated.png")
ple.normalize(individual=True)
ple.toPlot().save("normalized.png")            # false-color contour (log by default)

Decomposition and clustering

SpectrumSeries (and maps) can be decomposed into spectral components using scikit-learn (NMF / PCA / LDA) — useful for separating overlapping signals or removing a background.

# find / inspect components
decomp, components = series.getDominatingComponents(algo="LDA", n=3)
components.toPlot().save("components.png")
strongest = series.getDominatingComponent(algo="PCA")     # -> Spectrum

# remove a background learned from a reference, then apply it to the signal
nmf = reference.removeDominatingComponents(n=10)   # fit + subtract, returns the decomposition
signal.removeComponents(nmf)                        # subtract the same components

# project each series onto a fitted decomposition
weights = series.getProjection(decomp)              # -> SpectrumSeries of component weights

Plotting

Datasets expose toPlot() (and sometimes toScatterPlot()), which return a plot object. Call .save(path) to write an image, .show() to open a window, or styling methods first.

plot = series.toPlot()           # SpectrumSeriesPlot (overlaid, colored by key)
plot.setYLog()                   # log y-axis
plot.setXRange(400, 700)         # zoom x
plot.addLine(extra.data.values, "extra")   # overlay another curve
plot.legend()
plot.save("series.png", dpi=300, bbox_inches="tight")   # kwargs go to matplotlib savefig

You can also construct plot classes directly from fepydas.constructors.Plots, e.g. MultiSpectrumScatterPlot([a, b], ["A", "B"]) to overlay spectra with different axes, or SpectrumScatterPlotWithFit(spectrum, fit) to show data and fit together. Common styling methods on every plot: setXRange, setYRange, setXLog, setYLog, legend, show, save.

Exporting and saving

spectrum.export("spectrum.txt")     # tab-separated text (X + Y columns, with names/units)
series.export("series.txt")         # shared X plus one column per key
map.export("map.txt")               # one column per y position

spectrum.saveBinary("spectrum.bin") # pickle the whole object
# reload with Binary("spectrum.bin")

For fit output alongside the data, use the writers in fepydas.workers.ASCIIWriter:

from fepydas.workers.ASCIIWriter import SpectrumWithFitWriter
SpectrumWithFitWriter("plot.txt", spectrum, fit, withIRF=True)   # writes data + fit (+ IRF)

Batch processing whole directories

DirectoryLoader applies one reader to every file in a folder and lets you run the same operation on all of them:

from fepydas.workers.DirectoryLoader import DirectoryLoader
from fepydas.constructors.datareaders.Custom import ASCII_Gelbes
from fepydas.datatypes.Dataset import JoinSpectrumSeries

loader = DirectoryLoader(ASCII_Gelbes)
loader.loadFromDirectory("data/", lambda name: name.endswith(".dat"))
loader.callFunction("cutAxisByIdx", 50, 3000)      # call .cutAxisByIdx(50, 3000) on each

series = JoinSpectrumSeries(loader.getList(), interpol=True, names=loader.getKeys())

Rebinning a series

from fepydas.constructors.Binner import LinearBinner   # also LogBinner, InverseBinner
from fepydas.datatypes.Dataset import BinnedSpectrumSeries

binned = BinnedSpectrumSeries(series, LinearBinner(), num_bins=200)

Batch TRPL lifetime fitting

fepydas.workers.TRPL turns a folder of time-resolved measurements into a single tidy results table, choosing the number of exponential components per decay automatically by BIC (Bayesian Information Criterion). Files are loaded by extension (.plh → Automatic, .ptu → pulsed-decay extraction, .sdt/.phu/.thi → Histogram), so heterogeneous formats mix freely.

from fepydas.workers.TRPL import fit_trpl_batch

table = fit_trpl_batch(
    "trpl_data/",                 # a folder or an explicit list of files
    buckets=(("fast", 0, 2), ("mid", 2, 20), ("slow", 20, float("inf"))),
    fraction_floor=0.02,          # discard components below 2 % amplitude
    out_csv="trpl_results.csv",
    qc_pdf="trpl_qc.pdf",         # one data+fit+residual panel per measurement
    qc_redchi_max=2.0,            # flag poor fits (red frame) above this reduced chi-squared
)

The table has one row per (file × component): file, sample, T_K, feature_nm, component, bucket, tau, amp, fraction, tau_avg, N, BIC, redchi. Lifetimes are in the decay's own time unit (ns for histograms, seconds for gated .ptu), so choose buckets to match.

Fit a single decay directly with auto_exp_fit — it tries N = 1–3, keeps the lowest-BIC model, and guards against spurious components with per-component lifetime bounds and a minimum amplitude fraction:

from fepydas.constructors.datareaders.Generic import HistogramWithIRF
from fepydas.workers.TRPL import auto_exp_fit

decay = HistogramWithIRF("signal.sdt", "irf.sdt")   # the IRF is convolved automatically
result = auto_exp_fit(decay, fraction_floor=0.02)
print(result.N, [c.tau for c in result.components], result.redchi)

The tidy table feeds straight into your own trend plots (e.g. lifetime vs temperature, faceted by sample/bucket) with pandas/matplotlib or seaborn.

Gated .ptu measurements (pulsed excitation) are reduced to a single accumulated decay by Generic.PulseDecay / Proprietary.PTUDecay, which detect the pulse edges (threshold at the midpoint of the on/off levels) and sum the windows; the pure-NumPy steps live in fepydas.libs.pulse_windows. Pass diagnostic_plot=True (or a file path) to inspect the edge detection and the summed decay:

from fepydas.constructors.datareaders.Generic import PulseDecay

decay = PulseDecay("run.ptu", bin_width_s=0.01, diagnostic_plot="edges.png")

Worked examples

The examples/ directory contains complete, runnable scripts (each with its own data files). Install the package first with pip install -e . (or set PYTHONPATH to the repository root), then run a single example from its folder, e.g. cd examples/PLE && python3 PLE.py, or run them all with cd examples && python3 RunExamples.py.

API map

Every class and function carries a detailed docstring (visible as a tooltip in most editors) describing its parameters, options and return value.

Requirements

• Python 3 (3.9+ recommended)
• numpy ≥ 2.0, scipy, matplotlib, lmfit, scikit-learn, pyamg, sympy ≥ 1.13.3, pathos, tqdm, ptufile, pandas

These are installed automatically by pip install fepydas.

fepydas is developed at TU Berlin. Bug reports and contributions: https://git.tu-berlin.de/nippert/fepydas.