This video covers the same content as the text below. Choose your preferred learning format.
What you will learn in this chapter: This chapter integrates all the spectroscopic analysis techniques learned in Chapters 1-5 (IR, Raman, UV-Vis, XPS) and builds practical Python data analysis workflows. You will develop a universal spectral data loader, automated peak detection algorithms, machine learning-based classification, batch processing systems, and interactive visualization tools. All code is provided in reusable module format that can be readily applied to your own research data.
Spectrometers export data in a wide variety of formats (CSV, TXT, instrument-specific binary formats, and so on). By designing a general-purpose loader class, all of these formats can be handled through a single, uniform interface.
import numpy as np
import pandas as pd
from pathlib import Path
from typing import Union, Tuple, Optional
class SpectralData:
"""
Universal spectral data container.
Attributes:
-----------
x : array
Abscissa (wavelength, wavenumber, binding energy, etc.)
y : array
Ordinate (intensity, absorbance, etc.)
x_label : str
Label of the horizontal axis
y_label : str
Label of the vertical axis
metadata : dict
Metadata (acquisition date, instrument, conditions, etc.)
"""
def __init__(self, x: np.ndarray, y: np.ndarray,
x_label: str = "X", y_label: str = "Intensity",
metadata: Optional[dict] = None):
"""
Initialize a spectrum.
Parameters:
-----------
x : array
Abscissa data
y : array
Ordinate data
x_label : str
Label of the horizontal axis
y_label : str
Label of the vertical axis
metadata : dict, optional
Metadata
"""
self.x = np.asarray(x)
self.y = np.asarray(y)
self.x_label = x_label
self.y_label = y_label
self.metadata = metadata or {}
# Validate data size
if len(self.x) != len(self.y):
raise ValueError(f"x and y must have the same length: {len(x)} != {len(y)}")
def __repr__(self):
return (f"SpectralData(n_points={len(self.x)}, "
f"x_range=[{self.x.min():.2f}, {self.x.max():.2f}], "
f"y_range=[{self.y.min():.2f}, {self.y.max():.2f}])")
def copy(self):
"""Create a copy of the data."""
return SpectralData(self.x.copy(), self.y.copy(),
self.x_label, self.y_label,
self.metadata.copy())
def trim(self, x_min: float, x_max: float):
"""
Crop the data to the specified range.
Parameters:
-----------
x_min, x_max : float
Range to keep
"""
mask = (self.x >= x_min) & (self.x <= x_max)
self.x = self.x[mask]
self.y = self.y[mask]
def normalize(self, method: str = 'max'):
"""
Normalize the spectrum.
Parameters:
-----------
method : str
'max': normalize by the maximum value
'minmax': rescale to the 0-1 range
'area': normalize by the integrated area
"""
if method == 'max':
self.y = self.y / self.y.max()
elif method == 'minmax':
self.y = (self.y - self.y.min()) / (self.y.max() - self.y.min())
elif method == 'area':
area = np.trapz(self.y, self.x)
self.y = self.y / area
else:
raise ValueError(f"Unknown normalization method: {method}")
class SpectralDataLoader:
"""
Universal spectral data loader.
"""
@staticmethod
def load_csv(filepath: Union[str, Path], x_col: int = 0, y_col: int = 1,
delimiter: str = ',', skiprows: int = 0,
x_label: str = "X", y_label: str = "Intensity") -> SpectralData:
"""
Load spectral data from a CSV file.
Parameters:
-----------
filepath : str or Path
Path to the file
x_col : int
Column index of the abscissa data
y_col : int
Column index of the ordinate data
delimiter : str
Field delimiter
skiprows : int
Number of header rows to skip
x_label, y_label : str
Axis labels
Returns:
--------
data : SpectralData
Spectral data object
"""
df = pd.read_csv(filepath, delimiter=delimiter, skiprows=skiprows, header=None)
x = df.iloc[:, x_col].values
y = df.iloc[:, y_col].values
metadata = {
'filename': Path(filepath).name,
'source': 'CSV',
'columns': f"x={x_col}, y={y_col}"
}
return SpectralData(x, y, x_label, y_label, metadata)
@staticmethod
def load_txt(filepath: Union[str, Path], x_col: int = 0, y_col: int = 1,
skiprows: int = 0, x_label: str = "X", y_label: str = "Intensity") -> SpectralData:
"""
Load spectral data from a whitespace-delimited text file.
Parameters:
-----------
filepath : str or Path
Path to the file
x_col, y_col : int
Column indices of the abscissa and ordinate data
skiprows : int
Number of header rows to skip
x_label, y_label : str
Axis labels
Returns:
--------
data : SpectralData
Spectral data object
"""
data_array = np.loadtxt(filepath, skiprows=skiprows)
x = data_array[:, x_col]
y = data_array[:, y_col]
metadata = {
'filename': Path(filepath).name,
'source': 'TXT',
'columns': f"x={x_col}, y={y_col}"
}
return SpectralData(x, y, x_label, y_label, metadata)
@staticmethod
def auto_load(filepath: Union[str, Path], **kwargs) -> SpectralData:
"""
Select the appropriate loader automatically from the file extension.
Parameters:
-----------
filepath : str or Path
Path to the file
**kwargs
Extra arguments forwarded to the individual loaders
Returns:
--------
data : SpectralData
Spectral data object
"""
filepath = Path(filepath)
ext = filepath.suffix.lower()
if ext == '.csv':
return SpectralDataLoader.load_csv(filepath, **kwargs)
elif ext in ['.txt', '.dat']:
return SpectralDataLoader.load_txt(filepath, **kwargs)
else:
raise ValueError(f"Unsupported file format: {ext}")
# Usage example
if __name__ == "__main__":
# Generate and save simulated data
x_sim = np.linspace(400, 700, 300)
y_sim = 0.8 * np.exp(-((x_sim - 550)**2) / (2 * 40**2)) + np.random.normal(0, 0.02, len(x_sim))
# Save as a CSV file
import tempfile
with tempfile.NamedTemporaryFile(mode='w', suffix='.csv', delete=False) as f:
for xi, yi in zip(x_sim, y_sim):
f.write(f"{xi},{yi}\n")
temp_path = f.name
# Load the data
loader = SpectralDataLoader()
spectrum = loader.auto_load(temp_path, x_label="Wavelength (nm)", y_label="Absorbance")
print(spectrum)
print(f"Metadata: {spectrum.metadata}")
# Manipulate the data
spectrum_copy = spectrum.copy()
spectrum_copy.trim(500, 600)
spectrum_copy.normalize(method='max')
print(f"Trimmed and normalized: {spectrum_copy}")
# Clean up
import os
os.unlink(temp_path)
Detecting peaks automatically is essential when many samples must be analyzed efficiently. Here we extend scipy.signal.find_peaks into a peak detection class tailored to spectroscopic analysis.
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import find_peaks, peak_widths, peak_prominences
from scipy.optimize import curve_fit
from dataclasses import dataclass
from typing import List, Tuple
@dataclass
class Peak:
"""
Data class holding the properties of a single peak.
Attributes:
-----------
position : float
Peak position (abscissa coordinate)
height : float
Peak height
width : float
Full width at half maximum (FWHM)
area : float
Peak area
prominence : float
Peak prominence
"""
position: float
height: float
width: float
area: float
prominence: float
def __repr__(self):
return (f"Peak(pos={self.position:.2f}, height={self.height:.2f}, "
f"FWHM={self.width:.2f}, area={self.area:.1f})")
class PeakDetector:
"""
Detect peaks automatically in spectral data.
"""
def __init__(self, spectral_data: SpectralData):
"""
Parameters:
-----------
spectral_data : SpectralData
Spectral data object
"""
self.data = spectral_data
self.peaks: List[Peak] = []
def detect_peaks(self, height: float = None, prominence: float = None,
distance: int = None, width: Tuple[float, float] = None) -> List[Peak]:
"""
Detect peaks.
Parameters:
-----------
height : float, optional
Minimum peak height
prominence : float, optional
Minimum prominence
distance : int, optional
Minimum separation between peaks (in data points)
width : tuple, optional
Allowed peak width range (min_width, max_width)
Returns:
--------
peaks : list of Peak
List of detected peaks
"""
# Peak detection
peak_indices, properties = find_peaks(self.data.y,
height=height,
prominence=prominence,
distance=distance,
width=width)
# Peak widths (FWHM)
widths, width_heights, left_ips, right_ips = peak_widths(
self.data.y, peak_indices, rel_height=0.5
)
# Peak prominences
prominences, _, _ = peak_prominences(self.data.y, peak_indices)
# Build the list of Peak objects
self.peaks = []
for i, idx in enumerate(peak_indices):
# Peak position (abscissa coordinate)
position = self.data.x[idx]
# Peak height
height_val = self.data.y[idx]
# FWHM (convert from data points to abscissa units)
dx = np.mean(np.diff(self.data.x))
fwhm = widths[i] * dx
# Approximate peak area (Gaussian approximation:
# area ~= height * FWHM * sqrt(pi / (4 ln2)))
area = height_val * fwhm * np.sqrt(np.pi / (4 * np.log(2)))
peak = Peak(
position=position,
height=height_val,
width=fwhm,
area=area,
prominence=prominences[i]
)
self.peaks.append(peak)
return self.peaks
def fit_peaks_gaussian(self, peak_region_width: float = 5.0) -> List[dict]:
"""
Fit each detected peak with a Gaussian function.
Parameters:
-----------
peak_region_width : float
Half-width of the fitting window around each peak
Returns:
--------
fit_results : list of dict
Fitting result for each peak
"""
fit_results = []
for peak in self.peaks:
# Extract the fitting region
mask = (self.data.x >= peak.position - peak_region_width) & \
(self.data.x <= peak.position + peak_region_width)
x_fit = self.data.x[mask]
y_fit = self.data.y[mask]
if len(x_fit) < 5:
continue
# Gaussian fit
def gaussian(x, A, mu, sigma):
return A * np.exp(-((x - mu)**2) / (2 * sigma**2))
try:
# Initial parameters
p0 = [peak.height, peak.position, peak.width / (2 * np.sqrt(2 * np.log(2)))]
popt, pcov = curve_fit(gaussian, x_fit, y_fit, p0=p0, maxfev=5000)
A_fit, mu_fit, sigma_fit = popt
fwhm_fit = 2 * np.sqrt(2 * np.log(2)) * sigma_fit
area_fit = A_fit * sigma_fit * np.sqrt(2 * np.pi)
fit_results.append({
'peak_position': mu_fit,
'amplitude': A_fit,
'fwhm': fwhm_fit,
'area': area_fit,
'fit_quality': 'success'
})
except RuntimeError:
fit_results.append({
'peak_position': peak.position,
'fit_quality': 'failed'
})
return fit_results
def plot_detected_peaks(self, show_labels: bool = True):
"""
Visualize the detected peaks.
Parameters:
-----------
show_labels : bool
Whether to annotate the peak positions
"""
fig, ax = plt.subplots(figsize=(12, 6))
# Plot the spectrum
ax.plot(self.data.x, self.data.y, 'b-', linewidth=1.5, label='Spectrum', alpha=0.7)
# Mark the peak positions
if self.peaks:
peak_positions = [p.position for p in self.peaks]
peak_heights = [p.height for p in self.peaks]
ax.plot(peak_positions, peak_heights, 'ro', markersize=10,
label=f'Detected Peaks (n={len(self.peaks)})', zorder=5)
# Peak labels
if show_labels:
for peak in self.peaks:
ax.annotate(f'{peak.position:.1f}',
xy=(peak.position, peak.height),
xytext=(0, 10), textcoords='offset points',
ha='center', fontsize=9, fontweight='bold',
bbox=dict(boxstyle='round,pad=0.3', facecolor='yellow', alpha=0.7))
ax.set_xlabel(self.data.x_label, fontsize=12)
ax.set_ylabel(self.data.y_label, fontsize=12)
ax.set_title('Automated Peak Detection Result', fontsize=14, fontweight='bold')
ax.legend()
ax.grid(alpha=0.3)
plt.tight_layout()
plt.show()
# Usage example
if __name__ == "__main__":
# Simulated spectrum with three peaks
x = np.linspace(400, 700, 600)
y = (0.8 * np.exp(-((x - 450)**2) / (2 * 30**2)) +
0.6 * np.exp(-((x - 550)**2) / (2 * 40**2)) +
0.9 * np.exp(-((x - 620)**2) / (2 * 25**2)) +
np.random.normal(0, 0.02, len(x)))
spectrum = SpectralData(x, y, x_label="Wavelength (nm)", y_label="Absorbance")
# Peak detection
detector = PeakDetector(spectrum)
peaks = detector.detect_peaks(height=0.2, prominence=0.15, distance=50)
print(f"Number of detected peaks: {len(peaks)}")
for i, peak in enumerate(peaks, 1):
print(f" Peak {i}: {peak}")
# Gaussian fitting
fit_results = detector.fit_peaks_gaussian(peak_region_width=50)
print("\nGaussian fitting results:")
for i, result in enumerate(fit_results, 1):
if result['fit_quality'] == 'success':
print(f" Peak {i}: pos={result['peak_position']:.2f}, "
f"FWHM={result['fwhm']:.2f}, area={result['area']:.1f}")
# Plot
detector.plot_detected_peaks(show_labels=True)
Before a spectrum can be fed to a machine learning model, suitable features must be extracted from it.
import numpy as np
import matplotlib.pyplot as plt
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.metrics import classification_report, confusion_matrix
from sklearn.decomposition import PCA
from typing import List
import seaborn as sns
class SpectralClassifier:
"""
Classify spectral data with machine learning.
"""
def __init__(self, n_estimators: int = 100, use_pca: bool = False, n_components: int = 10):
"""
Parameters:
-----------
n_estimators : int
Number of trees in the random forest
use_pca : bool
Whether to apply PCA preprocessing
n_components : int
Number of PCA components
"""
self.classifier = RandomForestClassifier(n_estimators=n_estimators,
max_depth=10,
random_state=42)
self.use_pca = use_pca
self.pca = PCA(n_components=n_components) if use_pca else None
self.class_names = None
def extract_features(self, spectra_list: List[SpectralData]) -> np.ndarray:
"""
Extract features from a list of spectra.
Parameters:
-----------
spectra_list : list of SpectralData
List of spectral data objects
Returns:
--------
features : array, shape (n_samples, n_features)
Feature matrix
"""
features = []
for spectrum in spectra_list:
# Feature set 1: the whole spectrum (resampled)
# Resample every spectrum onto a fixed-length grid (e.g. 100 points)
x_uniform = np.linspace(spectrum.x.min(), spectrum.x.max(), 100)
y_resampled = np.interp(x_uniform, spectrum.x, spectrum.y)
# Feature set 2: statistical descriptors
mean_val = np.mean(spectrum.y)
std_val = np.std(spectrum.y)
max_val = np.max(spectrum.y)
min_val = np.min(spectrum.y)
# Concatenate
feature_vector = np.concatenate([
y_resampled,
[mean_val, std_val, max_val, min_val]
])
features.append(feature_vector)
return np.array(features)
def train(self, X: np.ndarray, y: np.ndarray, class_names: List[str]):
"""
Train the classifier.
Parameters:
-----------
X : array, shape (n_samples, n_features)
Feature matrix
y : array, shape (n_samples,)
Class labels
class_names : list of str
Names of the classes
"""
self.class_names = class_names
# PCA preprocessing
if self.use_pca:
X = self.pca.fit_transform(X)
print(f"PCA: {X.shape[1]} components explain "
f"{self.pca.explained_variance_ratio_.sum():.2%} of variance")
# Training
self.classifier.fit(X, y)
# Cross-validation
cv_scores = cross_val_score(self.classifier, X, y, cv=5)
print(f"Cross-validation accuracy: {cv_scores.mean():.3f} +/- {cv_scores.std():.3f}")
def predict(self, X: np.ndarray) -> np.ndarray:
"""
Predict class labels.
Parameters:
-----------
X : array, shape (n_samples, n_features)
Feature matrix
Returns:
--------
y_pred : array, shape (n_samples,)
Predicted class labels
"""
if self.use_pca:
X = self.pca.transform(X)
return self.classifier.predict(X)
def evaluate(self, X_test: np.ndarray, y_test: np.ndarray):
"""
Evaluate the classifier on test data.
Parameters:
-----------
X_test : array
Test features
y_test : array
Test labels
"""
y_pred = self.predict(X_test)
# Classification report
print("\nClassification report:")
print(classification_report(y_test, y_pred, target_names=self.class_names))
# Confusion matrix
cm = confusion_matrix(y_test, y_pred)
plt.figure(figsize=(8, 6))
sns.heatmap(cm, annot=True, fmt='d', cmap='Blues',
xticklabels=self.class_names, yticklabels=self.class_names)
plt.xlabel('Predicted label', fontsize=12)
plt.ylabel('True label', fontsize=12)
plt.title('Confusion Matrix', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()
# Usage example: classifying the spectra of three materials
if __name__ == "__main__":
np.random.seed(42)
# Generate simulated data (3 classes x 30 samples each)
def generate_spectrum_class(peak_positions, peak_widths, n_samples=30):
spectra = []
x = np.linspace(400, 700, 300)
for _ in range(n_samples):
y = np.zeros_like(x)
for pos, width in zip(peak_positions, peak_widths):
# Add random jitter to the peak position
pos_noise = pos + np.random.normal(0, 5)
amplitude = np.random.uniform(0.7, 1.0)
y += amplitude * np.exp(-((x - pos_noise)**2) / (2 * width**2))
y += np.random.normal(0, 0.03, len(x))
spectra.append(SpectralData(x, y, "Wavelength (nm)", "Intensity"))
return spectra
# Class 1: single peak (550 nm)
class1_spectra = generate_spectrum_class([550], [40], n_samples=30)
# Class 2: two peaks (450 nm, 620 nm)
class2_spectra = generate_spectrum_class([450, 620], [30, 35], n_samples=30)
# Class 3: three peaks (480 nm, 550 nm, 650 nm)
class3_spectra = generate_spectrum_class([480, 550, 650], [25, 30, 30], n_samples=30)
# Combine the data
all_spectra = class1_spectra + class2_spectra + class3_spectra
labels = np.array([0]*30 + [1]*30 + [2]*30)
class_names = ['Material A (1 peak)', 'Material B (2 peaks)', 'Material C (3 peaks)']
# Feature extraction
classifier = SpectralClassifier(n_estimators=100, use_pca=False)
X = classifier.extract_features(all_spectra)
# Train / test split
X_train, X_test, y_train, y_test = train_test_split(X, labels, test_size=0.3,
random_state=42, stratify=labels)
# Training
classifier.train(X_train, y_train, class_names)
# Evaluation
classifier.evaluate(X_test, y_test)
print(f"\nTotal samples: {len(all_spectra)}")
print(f"Feature dimension: {X.shape[1]}")
print(f"Training samples: {len(X_train)}, Test samples: {len(X_test)}")
In a real laboratory it is routine to process tens or hundreds of spectral files, which calls for an efficient batch processing system.
import numpy as np
import pandas as pd
from pathlib import Path
from typing import List, Dict, Callable, Union
from concurrent.futures import ProcessPoolExecutor, as_completed
from tqdm import tqdm
import json
class BatchProcessor:
"""
Batch-process multiple spectral files.
"""
def __init__(self, data_dir: Union[str, Path]):
"""
Parameters:
-----------
data_dir : str or Path
Data directory
"""
self.data_dir = Path(data_dir)
self.results = []
def find_files(self, pattern: str = "*.csv") -> List[Path]:
"""
Search for files matching a pattern.
Parameters:
-----------
pattern : str
File pattern (e.g. "*.csv", "sample_*.txt")
Returns:
--------
files : list of Path
List of matching files
"""
files = list(self.data_dir.glob(pattern))
print(f"Found {len(files)} files matching '{pattern}'")
return files
def process_file(self, filepath: Path,
processing_func: Callable[[SpectralData], Dict]) -> Dict:
"""
Process a single file.
Parameters:
-----------
filepath : Path
Path to the file
processing_func : callable
Processing function (takes a SpectralData, returns a result dict)
Returns:
--------
result : dict
Processing result
"""
try:
# Load the data
loader = SpectralDataLoader()
spectrum = loader.auto_load(filepath)
# Run the processing function
result = processing_func(spectrum)
# Attach file information
result['filename'] = filepath.name
result['status'] = 'success'
return result
except Exception as e:
return {
'filename': filepath.name,
'status': 'failed',
'error': str(e)
}
def batch_process(self, files: List[Path],
processing_func: Callable[[SpectralData], Dict],
n_workers: int = 4) -> pd.DataFrame:
"""
Process multiple files in parallel.
Parameters:
-----------
files : list of Path
Files to process
processing_func : callable
Processing function
n_workers : int
Number of parallel workers
Returns:
--------
results_df : DataFrame
DataFrame of processing results
"""
results = []
# Parallel processing
with ProcessPoolExecutor(max_workers=n_workers) as executor:
# Submit the tasks
futures = {executor.submit(self.process_file, f, processing_func): f
for f in files}
# Progress display
for future in tqdm(as_completed(futures), total=len(files),
desc="Processing files"):
result = future.result()
results.append(result)
# Convert to a DataFrame
results_df = pd.DataFrame(results)
self.results = results_df
# Success / failure summary
success_count = (results_df['status'] == 'success').sum()
print(f"\nProcessing complete: {success_count}/{len(files)} successful")
return results_df
def save_results(self, output_path: Union[str, Path], format: str = 'csv'):
"""
Save the results.
Parameters:
-----------
output_path : str or Path
Output file path
format : str
One of 'csv', 'json', 'excel'
"""
if self.results is None or len(self.results) == 0:
print("No results to save")
return
output_path = Path(output_path)
if format == 'csv':
self.results.to_csv(output_path, index=False)
elif format == 'json':
self.results.to_json(output_path, orient='records', indent=2)
elif format == 'excel':
self.results.to_excel(output_path, index=False)
else:
raise ValueError(f"Unknown format: {format}")
print(f"Results saved to {output_path}")
# Usage example: implementing a batch job
def example_processing_function(spectrum: SpectralData) -> Dict:
"""
Example: peak detection plus summary statistics.
Parameters:
-----------
spectrum : SpectralData
Spectral data
Returns:
--------
result : dict
Analysis result
"""
# Peak detection
detector = PeakDetector(spectrum)
peaks = detector.detect_peaks(height=0.1, prominence=0.05, distance=20)
# Summary statistics
mean_intensity = np.mean(spectrum.y)
max_intensity = np.max(spectrum.y)
std_intensity = np.std(spectrum.y)
# Peak information
peak_positions = [p.position for p in peaks]
peak_heights = [p.height for p in peaks]
return {
'n_peaks': len(peaks),
'peak_positions': peak_positions,
'peak_heights': peak_heights,
'mean_intensity': mean_intensity,
'max_intensity': max_intensity,
'std_intensity': std_intensity,
'x_range_min': spectrum.x.min(),
'x_range_max': spectrum.x.max()
}
if __name__ == "__main__":
# Simulation: generate and process multiple files
import tempfile
import shutil
# Create a temporary directory
temp_dir = Path(tempfile.mkdtemp())
print(f"Temporary directory: {temp_dir}")
# Generate simulated data (10 files)
for i in range(10):
x = np.linspace(400, 700, 300)
# Random peaks
n_peaks = np.random.randint(1, 4)
y = np.zeros_like(x)
for _ in range(n_peaks):
peak_pos = np.random.uniform(450, 650)
peak_width = np.random.uniform(20, 40)
y += np.random.uniform(0.5, 1.0) * np.exp(-((x - peak_pos)**2) / (2 * peak_width**2))
y += np.random.normal(0, 0.02, len(x))
# Save as CSV
filepath = temp_dir / f"sample_{i:03d}.csv"
with open(filepath, 'w') as f:
for xi, yi in zip(x, y):
f.write(f"{xi},{yi}\n")
# Batch processing
processor = BatchProcessor(temp_dir)
files = processor.find_files("*.csv")
results_df = processor.batch_process(files, example_processing_function, n_workers=2)
print("\nExcerpt of the results:")
print(results_df.head())
# Save the results
output_path = temp_dir / "batch_results.csv"
processor.save_results(output_path, format='csv')
# Clean up
shutil.rmtree(temp_dir)
print(f"\nTemporary directory removed")
In addition to static Matplotlib figures, the interactive visualization library Plotly lets you build dynamic plots with zooming, panning, and hover readouts.
import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots
from typing import List, Dict
class InteractiveSpectralViewer:
"""
Interactive spectral visualization (built on Plotly).
"""
def __init__(self):
self.fig = None
def plot_single_spectrum(self, spectrum: SpectralData, title: str = "Spectrum"):
"""
Plot a single spectrum.
Parameters:
-----------
spectrum : SpectralData
Spectral data
title : str
Plot title
"""
fig = go.Figure()
fig.add_trace(go.Scatter(
x=spectrum.x,
y=spectrum.y,
mode='lines',
name='Spectrum',
line=dict(color='blue', width=2),
hovertemplate='%{x:.2f}<br>%{y:.3f}<extra></extra>'
))
fig.update_layout(
title=title,
xaxis_title=spectrum.x_label,
yaxis_title=spectrum.y_label,
hovermode='x unified',
template='plotly_white',
width=1000,
height=600
)
self.fig = fig
return fig
def plot_multiple_spectra(self, spectra_list: List[SpectralData],
labels: List[str], title: str = "Multiple Spectra"):
"""
Overlay several spectra in one plot.
Parameters:
-----------
spectra_list : list of SpectralData
List of spectra
labels : list of str
Label for each spectrum
title : str
Plot title
"""
fig = go.Figure()
colors = ['blue', 'red', 'green', 'orange', 'purple', 'brown']
for i, (spectrum, label) in enumerate(zip(spectra_list, labels)):
color = colors[i % len(colors)]
fig.add_trace(go.Scatter(
x=spectrum.x,
y=spectrum.y,
mode='lines',
name=label,
line=dict(color=color, width=2),
hovertemplate=f'{label}<br>%{{x:.2f}}<br>%{{y:.3f}}<extra></extra>'
))
fig.update_layout(
title=title,
xaxis_title=spectra_list[0].x_label,
yaxis_title=spectra_list[0].y_label,
hovermode='x unified',
template='plotly_white',
width=1000,
height=600,
legend=dict(x=0.02, y=0.98, bgcolor='rgba(255,255,255,0.8)')
)
self.fig = fig
return fig
def plot_with_peaks(self, spectrum: SpectralData, peaks: List[Peak],
title: str = "Spectrum with Detected Peaks"):
"""
Plot a spectrum together with its detected peaks.
Parameters:
-----------
spectrum : SpectralData
Spectral data
peaks : list of Peak
Detected peaks
title : str
Plot title
"""
fig = go.Figure()
# Spectrum
fig.add_trace(go.Scatter(
x=spectrum.x,
y=spectrum.y,
mode='lines',
name='Spectrum',
line=dict(color='blue', width=2),
hovertemplate='%{x:.2f}<br>%{y:.3f}<extra></extra>'
))
# Peaks
if peaks:
peak_x = [p.position for p in peaks]
peak_y = [p.height for p in peaks]
peak_info = [f"Position: {p.position:.2f}<br>Height: {p.height:.3f}<br>FWHM: {p.width:.2f}"
for p in peaks]
fig.add_trace(go.Scatter(
x=peak_x,
y=peak_y,
mode='markers',
name=f'Peaks (n={len(peaks)})',
marker=dict(color='red', size=12, symbol='diamond'),
text=peak_info,
hovertemplate='%{text}<extra></extra>'
))
fig.update_layout(
title=title,
xaxis_title=spectrum.x_label,
yaxis_title=spectrum.y_label,
hovermode='closest',
template='plotly_white',
width=1000,
height=600
)
self.fig = fig
return fig
def plot_comparison_grid(self, spectra_dict: Dict[str, SpectralData],
title: str = "Spectral Comparison"):
"""
Display several spectra on a grid of subplots.
Parameters:
-----------
spectra_dict : dict
Dictionary of {label: SpectralData}
title : str
Plot title
"""
n_spectra = len(spectra_dict)
rows = (n_spectra + 1) // 2
cols = 2
fig = make_subplots(
rows=rows, cols=cols,
subplot_titles=list(spectra_dict.keys()),
vertical_spacing=0.12,
horizontal_spacing=0.1
)
for i, (label, spectrum) in enumerate(spectra_dict.items(), 1):
row = (i - 1) // cols + 1
col = (i - 1) % cols + 1
fig.add_trace(
go.Scatter(
x=spectrum.x,
y=spectrum.y,
mode='lines',
name=label,
line=dict(width=2),
showlegend=False
),
row=row, col=col
)
fig.update_xaxes(title_text=spectrum.x_label, row=row, col=col)
fig.update_yaxes(title_text=spectrum.y_label, row=row, col=col)
fig.update_layout(
title_text=title,
template='plotly_white',
width=1200,
height=300 * rows
)
self.fig = fig
return fig
def show(self):
"""Display the current plot."""
if self.fig:
self.fig.show()
# Usage example
if __name__ == "__main__":
# Generate simulated data
x = np.linspace(400, 700, 500)
# Spectrum 1
y1 = 0.8 * np.exp(-((x - 500)**2) / (2 * 40**2)) + np.random.normal(0, 0.02, len(x))
spectrum1 = SpectralData(x, y1, "Wavelength (nm)", "Absorbance")
# Spectrum 2
y2 = 0.6 * np.exp(-((x - 550)**2) / (2 * 35**2)) + np.random.normal(0, 0.02, len(x))
spectrum2 = SpectralData(x, y2, "Wavelength (nm)", "Absorbance")
# Spectrum 3
y3 = 0.9 * np.exp(-((x - 600)**2) / (2 * 30**2)) + np.random.normal(0, 0.02, len(x))
spectrum3 = SpectralData(x, y3, "Wavelength (nm)", "Absorbance")
# Interactive viewer
viewer = InteractiveSpectralViewer()
# Single spectrum
print("1. Plotting a single spectrum")
fig1 = viewer.plot_single_spectrum(spectrum1, title="Sample Spectrum")
# fig1.show() # Opens in the browser
# Overlay of several spectra
print("2. Overlaying several spectra")
fig2 = viewer.plot_multiple_spectra(
[spectrum1, spectrum2, spectrum3],
['Sample A (500 nm)', 'Sample B (550 nm)', 'Sample C (600 nm)'],
title="Multiple Spectra Overlay"
)
# fig2.show()
# With detected peaks
print("3. Showing detected peaks")
detector = PeakDetector(spectrum1)
peaks = detector.detect_peaks(height=0.2, prominence=0.1)
fig3 = viewer.plot_with_peaks(spectrum1, peaks, title="Spectrum with Peaks")
# fig3.show()
# Grid display
print("4. Grid display")
spectra_dict = {
'Sample A': spectrum1,
'Sample B': spectrum2,
'Sample C': spectrum3
}
fig4 = viewer.plot_comparison_grid(spectra_dict, title="Spectral Comparison Grid")
# fig4.show()
print("\nHow to use the interactive plots:")
print(" - Hover with the mouse to read out data point values")
print(" - Drag to zoom, double-click to reset")
print(" - Click a legend entry to show/hide that series")
We now combine the components built so far into a complete workflow that runs from data loading through to result export.
import numpy as np
import pandas as pd
from pathlib import Path
from typing import Dict, List
import json
class SpectralAnalysisPipeline:
"""
Integrated spectral analysis pipeline.
"""
def __init__(self, config: Dict):
"""
Parameters:
-----------
config : dict
Pipeline configuration
{
'data_dir': Path,
'output_dir': Path,
'file_pattern': str,
'preprocessing': {...},
'peak_detection': {...},
'analysis_type': 'classification' or 'quantification'
}
"""
self.config = config
self.data_dir = Path(config['data_dir'])
self.output_dir = Path(config['output_dir'])
self.output_dir.mkdir(parents=True, exist_ok=True)
self.spectra = []
self.results = None
def run_pipeline(self):
"""
Run the whole pipeline.
"""
print("="*60)
print("Spectral Analysis Pipeline")
print("="*60)
# Step 1: load the data
print("\n[Step 1/5] Loading data...")
self._load_data()
# Step 2: preprocessing
print("\n[Step 2/5] Preprocessing...")
self._preprocess()
# Step 3: peak detection
print("\n[Step 3/5] Peak detection...")
self._detect_peaks()
# Step 4: analysis
print("\n[Step 4/5] Analysis...")
self._analyze()
# Step 5: export the results
print("\n[Step 5/5] Exporting results...")
self._export_results()
print("\n" + "="*60)
print("Pipeline completed successfully!")
print("="*60)
def _load_data(self):
"""Load the data."""
processor = BatchProcessor(self.data_dir)
files = processor.find_files(self.config.get('file_pattern', '*.csv'))
loader = SpectralDataLoader()
for file in files:
try:
spectrum = loader.auto_load(file)
spectrum.metadata['source_file'] = file.name
self.spectra.append(spectrum)
except Exception as e:
print(f" Warning: Failed to load {file.name}: {e}")
print(f" Loaded {len(self.spectra)} spectra")
def _preprocess(self):
"""Preprocess the spectra."""
preproc_config = self.config.get('preprocessing', {})
for spectrum in self.spectra:
# Denoising
if preproc_config.get('smooth', False):
from scipy.signal import savgol_filter
spectrum.y = savgol_filter(spectrum.y, window_length=11, polyorder=3)
# Trimming
if 'trim_range' in preproc_config:
x_min, x_max = preproc_config['trim_range']
spectrum.trim(x_min, x_max)
# Normalization
if 'normalize' in preproc_config:
spectrum.normalize(method=preproc_config['normalize'])
print(f" Preprocessed {len(self.spectra)} spectra")
def _detect_peaks(self):
"""Detect peaks."""
peak_config = self.config.get('peak_detection', {})
for spectrum in self.spectra:
detector = PeakDetector(spectrum)
peaks = detector.detect_peaks(
height=peak_config.get('height'),
prominence=peak_config.get('prominence'),
distance=peak_config.get('distance')
)
spectrum.metadata['peaks'] = peaks
total_peaks = sum(len(s.metadata['peaks']) for s in self.spectra)
print(f" Detected {total_peaks} peaks across {len(self.spectra)} spectra")
def _analyze(self):
"""Analyze the spectra."""
analysis_type = self.config.get('analysis_type', 'quantification')
if analysis_type == 'quantification':
self._analyze_quantification()
elif analysis_type == 'classification':
self._analyze_classification()
else:
raise ValueError(f"Unknown analysis type: {analysis_type}")
def _analyze_quantification(self):
"""Quantitative analysis."""
results = []
for spectrum in self.spectra:
peaks = spectrum.metadata.get('peaks', [])
result = {
'filename': spectrum.metadata.get('source_file', 'unknown'),
'n_peaks': len(peaks),
'peak_positions': [p.position for p in peaks],
'peak_areas': [p.area for p in peaks],
'total_area': sum(p.area for p in peaks)
}
results.append(result)
self.results = pd.DataFrame(results)
print(f" Quantification completed for {len(results)} samples")
def _analyze_classification(self):
"""Classification analysis."""
# Feature extraction
classifier = SpectralClassifier(n_estimators=100, use_pca=False)
X = classifier.extract_features(self.spectra)
# For simplicity we use clustering here (unsupervised classification)
from sklearn.cluster import KMeans
n_clusters = self.config.get('n_clusters', 3)
kmeans = KMeans(n_clusters=n_clusters, random_state=42)
labels = kmeans.fit_predict(X)
results = []
for i, spectrum in enumerate(self.spectra):
result = {
'filename': spectrum.metadata.get('source_file', 'unknown'),
'cluster': int(labels[i]),
'n_peaks': len(spectrum.metadata.get('peaks', []))
}
results.append(result)
self.results = pd.DataFrame(results)
print(f" Classification completed: {n_clusters} clusters identified")
def _export_results(self):
"""Export the results."""
if self.results is None:
print(" No results to export")
return
# CSV export
csv_path = self.output_dir / 'results.csv'
self.results.to_csv(csv_path, index=False)
print(f" Saved CSV: {csv_path}")
# JSON export
json_path = self.output_dir / 'results.json'
self.results.to_json(json_path, orient='records', indent=2)
print(f" Saved JSON: {json_path}")
# Summary report
summary = {
'total_samples': len(self.spectra),
'analysis_type': self.config.get('analysis_type'),
'output_files': [str(csv_path), str(json_path)]
}
summary_path = self.output_dir / 'summary.json'
with open(summary_path, 'w') as f:
json.dump(summary, f, indent=2)
print(f" Saved summary: {summary_path}")
# Usage example
if __name__ == "__main__":
import tempfile
import shutil
# Create a temporary directory
temp_dir = Path(tempfile.mkdtemp())
data_dir = temp_dir / 'data'
output_dir = temp_dir / 'output'
data_dir.mkdir()
# Generate simulated data
for i in range(15):
x = np.linspace(400, 700, 300)
y = np.zeros_like(x)
n_peaks = np.random.randint(1, 4)
for _ in range(n_peaks):
pos = np.random.uniform(450, 650)
width = np.random.uniform(20, 40)
y += np.random.uniform(0.5, 1.0) * np.exp(-((x - pos)**2) / (2 * width**2))
y += np.random.normal(0, 0.02, len(x))
filepath = data_dir / f"sample_{i:03d}.csv"
with open(filepath, 'w') as f:
for xi, yi in zip(x, y):
f.write(f"{xi},{yi}\n")
# Pipeline configuration
config = {
'data_dir': data_dir,
'output_dir': output_dir,
'file_pattern': '*.csv',
'preprocessing': {
'smooth': True,
'trim_range': (420, 680),
'normalize': 'max'
},
'peak_detection': {
'height': 0.2,
'prominence': 0.1,
'distance': 20
},
'analysis_type': 'quantification', # or 'classification'
'n_clusters': 3
}
# Run the pipeline
pipeline = SpectralAnalysisPipeline(config)
pipeline.run_pipeline()
# Inspect the results
print("\nExcerpt of the results:")
print(pipeline.results.head())
# Clean up
shutil.rmtree(temp_dir)
The following template lays out a project structure for analyzing real research data with the tools built in this chapter.
import os
from pathlib import Path
from typing import Optional
class SpectralProjectTemplate:
"""
Generate a template for a spectral analysis project.
"""
@staticmethod
def create_project(project_name: str, base_dir: Optional[Path] = None):
"""
Create the project directory structure.
Parameters:
-----------
project_name : str
Project name
base_dir : Path, optional
Base directory (defaults to the current directory)
"""
if base_dir is None:
base_dir = Path.cwd()
project_dir = base_dir / project_name
project_dir.mkdir(parents=True, exist_ok=True)
# Directory structure
dirs = [
'data/raw',
'data/processed',
'results',
'figures',
'notebooks',
'scripts',
'config'
]
for dir_path in dirs:
(project_dir / dir_path).mkdir(parents=True, exist_ok=True)
# README file
readme_content = f"""# {project_name}
## Project Structure
- `data/raw/`: raw data files (originals, never edit)
- `data/processed/`: preprocessed data
- `results/`: analysis results (CSV, JSON, ...)
- `figures/`: plots and figures
- `notebooks/`: Jupyter notebooks
- `scripts/`: Python analysis scripts
- `config/`: configuration files
## Workflow
1. Place the raw data in `data/raw/`
2. Run the preprocessing with `scripts/preprocess.py`
3. Run the analysis with `scripts/analyze.py`
4. Results are written to `results/` and figures to `figures/`
## Requirements
```bash
pip install numpy scipy matplotlib pandas scikit-learn plotly
```
## Quick Start
```python
python scripts/analyze.py --config config/config.json
```
"""
with open(project_dir / 'README.md', 'w') as f:
f.write(readme_content)
# Configuration file template
config_content = """{
"data_dir": "data/raw",
"output_dir": "results",
"file_pattern": "*.csv",
"preprocessing": {
"smooth": true,
"trim_range": [400, 700],
"normalize": "max"
},
"peak_detection": {
"height": 0.2,
"prominence": 0.1,
"distance": 20
},
"analysis_type": "quantification"
}
"""
with open(project_dir / 'config/config.json', 'w') as f:
f.write(config_content)
# Analysis script template
analyze_script = """#!/usr/bin/env python
# -*- coding: utf-8 -*-
\"\"\"
Main spectral analysis script
\"\"\"
import argparse
import json
from pathlib import Path
# Import the SpectralAnalysisPipeline class from this chapter below
# from spectral_tools import SpectralAnalysisPipeline
def main():
parser = argparse.ArgumentParser(description='Spectral Analysis Pipeline')
parser.add_argument('--config', type=str, required=True,
help='Path to config JSON file')
args = parser.parse_args()
# Load the configuration
with open(args.config, 'r') as f:
config = json.load(f)
# Run the pipeline
# pipeline = SpectralAnalysisPipeline(config)
# pipeline.run_pipeline()
print("Analysis completed!")
if __name__ == "__main__":
main()
"""
with open(project_dir / 'scripts/analyze.py', 'w') as f:
f.write(analyze_script)
# .gitignore file
gitignore_content = """# Data files
data/raw/*
!data/raw/.gitkeep
*.csv
*.txt
*.dat
# Python
__pycache__/
*.py[cod]
*.so
.ipynb_checkpoints/
# Results
results/*
!results/.gitkeep
figures/*
!figures/.gitkeep
"""
with open(project_dir / '.gitignore', 'w') as f:
f.write(gitignore_content)
# .gitkeep files (so that empty directories stay under version control)
for dir_path in ['data/raw', 'results', 'figures']:
(project_dir / dir_path / '.gitkeep').touch()
print(f"Project '{project_name}' created successfully!")
print(f"Location: {project_dir}")
print(f"\nNext steps:")
print(f" 1. cd {project_dir}")
print(f" 2. Put your data files in data/raw/")
print(f" 3. Edit config/config.json")
print(f" 4. python scripts/analyze.py --config config/config.json")
# Usage example
if __name__ == "__main__":
# Create the project
SpectralProjectTemplate.create_project("my_spectroscopy_project")
Create a SpectralData object holding 100 points over the 400-700 nm range, trim it to the 500-600 nm range, and then normalize it by its maximum value.
Solution:
import numpy as np
# Generate the data
x = np.linspace(400, 700, 100)
y = np.exp(-((x - 550)**2) / (2 * 50**2))
# Create the SpectralData object
spectrum = SpectralData(x, y, "Wavelength (nm)", "Intensity")
# Trim
spectrum.trim(500, 600)
# Normalize
spectrum.normalize(method='max')
print(f"Number of points after trimming: {len(spectrum.x)}")
print(f"Maximum after normalization: {spectrum.y.max()}") # 1.0
From simulated data, detect the peaks with a height of at least 0.3 and a prominence of at least 0.2, and print the position and height of each peak.
Solution:
# Simulated data
x = np.linspace(400, 700, 500)
y = (0.8 * np.exp(-((x - 500)**2) / (2 * 30**2)) +
0.5 * np.exp(-((x - 600)**2) / (2 * 40**2)))
spectrum = SpectralData(x, y, "Wavelength (nm)", "Intensity")
# Peak detection
detector = PeakDetector(spectrum)
peaks = detector.detect_peaks(height=0.3, prominence=0.2)
print(f"Number of detected peaks: {len(peaks)}")
for i, peak in enumerate(peaks, 1):
print(f" Peak {i}: position={peak.position:.1f} nm, height={peak.height:.3f}")
Load spectral data from a two-column CSV file (wavelength, intensity) and print the number of data points and the wavelength range.
Solution:
# Data loader
loader = SpectralDataLoader()
spectrum = loader.load_csv('sample.csv', x_col=0, y_col=1,
x_label="Wavelength (nm)", y_label="Intensity")
print(f"Number of points: {len(spectrum.x)}")
print(f"Wavelength range: {spectrum.x.min():.1f} - {spectrum.x.max():.1f} nm")
print(f"Intensity range: {spectrum.y.min():.3f} - {spectrum.y.max():.3f}")
Write a function that applies the following preprocessing steps in sequence: (1) Savitzky-Golay smoothing, (2) linear baseline correction, (3) normalization by the maximum value.
Solution:
from scipy.signal import savgol_filter
from scipy.stats import linregress
def custom_preprocessing_pipeline(spectrum: SpectralData) -> SpectralData:
"""
Custom preprocessing pipeline.
Parameters:
-----------
spectrum : SpectralData
Input spectrum
Returns:
--------
processed_spectrum : SpectralData
Preprocessed spectrum
"""
# Work on a copy
processed = spectrum.copy()
# (1) Savitzky-Golay smoothing
processed.y = savgol_filter(processed.y, window_length=11, polyorder=3)
# (2) Baseline correction (linear fit)
# Estimate the baseline from the outer 10% of the data on each side
n = len(processed.x)
baseline_indices = list(range(int(n * 0.1))) + list(range(int(n * 0.9), n))
x_base = processed.x[baseline_indices]
y_base = processed.y[baseline_indices]
slope, intercept, _, _, _ = linregress(x_base, y_base)
baseline = slope * processed.x + intercept
processed.y = processed.y - baseline
# (3) Normalization by the maximum value
processed.normalize(method='max')
return processed
# Usage example
x = np.linspace(400, 700, 300)
y = 0.8 * np.exp(-((x - 550)**2) / (2 * 40**2)) + 0.02 * x + np.random.normal(0, 0.02, len(x))
spectrum = SpectralData(x, y, "Wavelength (nm)", "Intensity")
processed_spectrum = custom_preprocessing_pipeline(spectrum)
print("Preprocessing complete:")
print(f" Original y range: [{spectrum.y.min():.3f}, {spectrum.y.max():.3f}]")
print(f" Processed y range: [{processed_spectrum.y.min():.3f}, {processed_spectrum.y.max():.3f}]")
Note that the linear baseline used here is only the simplest option. Curved baselines are usually handled with asymmetric least squares (ALS), SNIP, or rubberband algorithms. None of these live in scipy.signal - they must be implemented by hand or taken from a dedicated package such as pybaselines.
Load several spectral files, compute the number of peaks, the maximum peak height, and the mean intensity for each file, and return the results as a DataFrame.
Solution:
from typing import Dict
def batch_statistical_analysis(spectrum: SpectralData) -> Dict:
"""
Statistical analysis of a spectrum.
Parameters:
-----------
spectrum : SpectralData
Spectral data
Returns:
--------
stats : dict
Summary statistics
"""
# Peak detection
detector = PeakDetector(spectrum)
peaks = detector.detect_peaks(height=0.1, prominence=0.05)
# Statistics
max_peak_height = max([p.height for p in peaks]) if peaks else 0.0
mean_intensity = np.mean(spectrum.y)
return {
'n_peaks': len(peaks),
'max_peak_height': max_peak_height,
'mean_intensity': mean_intensity,
'intensity_std': np.std(spectrum.y)
}
# Batch processing
processor = BatchProcessor('data_directory')
files = processor.find_files("*.csv")
results_df = processor.batch_process(files, batch_statistical_analysis, n_workers=2)
print("Statistical analysis results:")
print(results_df.describe())
From a series of spectra measured consecutively in time, plot how the height of a particular peak (e.g. at 550 nm) evolves.
Solution:
import matplotlib.pyplot as plt
def analyze_time_series_spectra(spectra_list: List[SpectralData],
target_wavelength: float = 550.0) -> np.ndarray:
"""
Extract the intensity at a given wavelength from a time series of spectra.
Parameters:
-----------
spectra_list : list of SpectralData
Time-ordered list of spectra
target_wavelength : float
Wavelength to track (nm)
Returns:
--------
intensities : array
Intensity at each time point
"""
intensities = []
for spectrum in spectra_list:
# Index of the closest wavelength
idx = np.argmin(np.abs(spectrum.x - target_wavelength))
intensities.append(spectrum.y[idx])
return np.array(intensities)
# Usage example (simulation)
time_points = np.arange(0, 60, 2) # 0-60 min, 2 min intervals
spectra_series = []
for t in time_points:
x = np.linspace(400, 700, 300)
# The peak height decays with time
decay = np.exp(-t / 30)
y = 0.8 * decay * np.exp(-((x - 550)**2) / (2 * 40**2))
spectra_series.append(SpectralData(x, y, "Wavelength (nm)", "Intensity"))
# Track the 550 nm peak
intensities = analyze_time_series_spectra(spectra_series, target_wavelength=550)
# Plot
plt.figure(figsize=(10, 6))
plt.plot(time_points, intensities, 'o-', linewidth=2, markersize=8)
plt.xlabel('Time (min)', fontsize=12)
plt.ylabel('Intensity @ 550 nm', fontsize=12)
plt.title('Time-Series Spectral Analysis: Peak Intensity vs Time',
fontsize=14, fontweight='bold')
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()
print(f"Initial intensity: {intensities[0]:.3f}")
print(f"Final intensity: {intensities[-1]:.3f}")
print(f"Decay: {(1 - intensities[-1]/intensities[0]) * 100:.1f}%")
Build a one-dimensional convolutional neural network (CNN) and train it to classify spectra into three classes. Use TensorFlow/Keras.
Solution:
# Requires TensorFlow/Keras: pip install tensorflow
import numpy as np
import matplotlib.pyplot as plt
import tensorflow as tf
from tensorflow import keras
from tensorflow.keras import layers
def build_1d_cnn(input_length: int, n_classes: int):
"""
Build a 1D CNN model.
Parameters:
-----------
input_length : int
Length of the input spectrum
n_classes : int
Number of classes
Returns:
--------
model : keras.Model
CNN model
"""
model = keras.Sequential([
# Input layer
layers.Input(shape=(input_length, 1)),
# Convolution block 1
layers.Conv1D(filters=32, kernel_size=7, activation='relu'),
layers.MaxPooling1D(pool_size=2),
layers.Dropout(0.2),
# Convolution block 2
layers.Conv1D(filters=64, kernel_size=5, activation='relu'),
layers.MaxPooling1D(pool_size=2),
layers.Dropout(0.2),
# Convolution block 3
layers.Conv1D(filters=128, kernel_size=3, activation='relu'),
layers.GlobalAveragePooling1D(),
# Fully connected layer
layers.Dense(64, activation='relu'),
layers.Dropout(0.3),
# Output layer
layers.Dense(n_classes, activation='softmax')
])
model.compile(
optimizer='adam',
loss='sparse_categorical_crossentropy',
metrics=['accuracy']
)
return model
# Data preparation (simulation)
X_train = [] # shape: (n_samples, seq_length, 1)
y_train = []
for class_label in range(3):
for _ in range(100):
x = np.linspace(0, 1, 100)
# Each class has a different spectral shape
if class_label == 0:
y = np.exp(-((x - 0.5)**2) / 0.1)
elif class_label == 1:
y = np.exp(-((x - 0.3)**2) / 0.05) + np.exp(-((x - 0.7)**2) / 0.05)
else:
y = np.sin(x * 10) * 0.5 + 0.5
y += np.random.normal(0, 0.05, len(x))
X_train.append(y)
y_train.append(class_label)
X_train = np.array(X_train).reshape(-1, 100, 1)
y_train = np.array(y_train)
# Train the model
model = build_1d_cnn(input_length=100, n_classes=3)
history = model.fit(X_train, y_train, epochs=50, batch_size=32,
validation_split=0.2, verbose=0)
print(f"Final training accuracy: {history.history['accuracy'][-1]:.3f}")
print(f"Final validation accuracy: {history.history['val_accuracy'][-1]:.3f}")
# Plot the learning curves
plt.figure(figsize=(10, 5))
plt.plot(history.history['accuracy'], label='Training accuracy')
plt.plot(history.history['val_accuracy'], label='Validation accuracy')
plt.xlabel('Epoch')
plt.ylabel('Accuracy')
plt.title('CNN Learning Curves')
plt.legend()
plt.grid(alpha=0.3)
plt.show()
Answer: the CNN reaches high accuracy (above 90%) on this classification task.
Implement a system that analyzes every new spectral file as soon as it appears in a watched directory and flags anomalies (deviations from the normal pattern).
Solution:
import time
from watchdog.observers import Observer
from watchdog.events import FileSystemEventHandler
class SpectralMonitor(FileSystemEventHandler):
"""
Spectral file watcher with anomaly detection.
"""
def __init__(self, reference_spectra: List[SpectralData], threshold: float = 0.3):
"""
Parameters:
-----------
reference_spectra : list of SpectralData
Reference spectra representing normal data
threshold : float
Anomaly detection threshold (distance)
"""
self.reference_spectra = reference_spectra
self.threshold = threshold
# Compute the mean of the reference spectra
X_ref = np.array([s.y for s in reference_spectra])
self.reference_mean = np.mean(X_ref, axis=0)
self.reference_std = np.std(X_ref, axis=0)
def on_created(self, event):
"""Callback fired when a file is created."""
if event.is_directory:
return
if event.src_path.endswith('.csv'):
print(f"\nNew file detected: {event.src_path}")
time.sleep(0.5) # Wait for the file write to finish
self.analyze_new_spectrum(event.src_path)
def analyze_new_spectrum(self, filepath: str):
"""Analyze the newly created spectrum."""
try:
# Load the data
loader = SpectralDataLoader()
spectrum = loader.auto_load(filepath)
# Anomaly detection (Euclidean distance)
distance = np.linalg.norm(spectrum.y - self.reference_mean)
if distance > self.threshold:
print(f"[WARN] Anomaly detected! distance: {distance:.3f} "
f"(threshold: {self.threshold:.3f})")
print(f" -> Sending alert: {Path(filepath).name}")
else:
print(f"[OK] Within the normal range. distance: {distance:.3f}")
except Exception as e:
print(f"[ERROR] {e}")
# Usage example (uncomment for real use)
"""
# Prepare the reference data
reference_spectra = [...] # List of normal spectra
# Start the monitoring system
monitor = SpectralMonitor(reference_spectra, threshold=0.3)
observer = Observer()
observer.schedule(monitor, path='./watch_directory', recursive=False)
observer.start()
print("Monitoring started...")
try:
while True:
time.sleep(1)
except KeyboardInterrupt:
observer.stop()
observer.join()
"""
print("Real-time monitoring system template implemented")
Answer: a real-time monitoring system built on the watchdog library.
Store a large collection of spectra in an SQLite database and implement a similar-spectrum search.
Solution:
import json
import sqlite3
import pickle
class SpectralDatabase:
"""
Manager class for a spectral database.
"""
def __init__(self, db_path: str = 'spectra.db'):
"""
Parameters:
-----------
db_path : str
Path to the database file
"""
self.db_path = db_path
self.conn = sqlite3.connect(db_path)
self._create_table()
def _create_table(self):
"""Create the table."""
cursor = self.conn.cursor()
cursor.execute("""
CREATE TABLE IF NOT EXISTS spectra (
id INTEGER PRIMARY KEY AUTOINCREMENT,
filename TEXT,
x_data BLOB,
y_data BLOB,
metadata TEXT,
created_at TIMESTAMP DEFAULT CURRENT_TIMESTAMP
)
""")
self.conn.commit()
def insert_spectrum(self, spectrum: SpectralData):
"""Insert a spectrum."""
cursor = self.conn.cursor()
cursor.execute("""
INSERT INTO spectra (filename, x_data, y_data, metadata)
VALUES (?, ?, ?, ?)
""", (
spectrum.metadata.get('source_file', 'unknown'),
pickle.dumps(spectrum.x),
pickle.dumps(spectrum.y),
json.dumps(spectrum.metadata)
))
self.conn.commit()
return cursor.lastrowid
def search_similar_spectra(self, query_spectrum: SpectralData, top_k: int = 5):
"""
Search for similar spectra.
Parameters:
-----------
query_spectrum : SpectralData
Query spectrum
top_k : int
Number of top matches to return
Returns:
--------
results : list of tuples
(id, filename, similarity_score)
"""
cursor = self.conn.cursor()
cursor.execute("SELECT id, filename, y_data FROM spectra")
results = []
for row in cursor.fetchall():
id, filename, y_blob = row
y_data = pickle.loads(y_blob)
# Cosine similarity
similarity = np.dot(query_spectrum.y, y_data) / \
(np.linalg.norm(query_spectrum.y) * np.linalg.norm(y_data))
results.append((id, filename, similarity))
# Sort by similarity
results.sort(key=lambda x: x[2], reverse=True)
return results[:top_k]
def close(self):
"""Close the database connection."""
self.conn.close()
# Usage example
db = SpectralDatabase('test_spectra.db')
# Insert spectra
for i in range(10):
x = np.linspace(400, 700, 100)
y = np.exp(-((x - (500 + i*10))**2) / (2 * 30**2))
spectrum = SpectralData(x, y, "Wavelength", "Intensity")
spectrum.metadata['source_file'] = f"sample_{i}.csv"
db.insert_spectrum(spectrum)
# Similarity search
query_x = np.linspace(400, 700, 100)
query_y = np.exp(-((query_x - 550)**2) / (2 * 30**2))
query_spectrum = SpectralData(query_x, query_y, "Wavelength", "Intensity")
similar_spectra = db.search_similar_spectra(query_spectrum, top_k=3)
print("Similar spectrum search results:")
for id, filename, score in similar_spectra:
print(f" ID: {id}, file: {filename}, similarity: {score:.3f}")
db.close()
Answer: a spectral database and similarity search built on SQLite.
Use the following list to assess your own progress:
SpectralData classpybaselines.