plasma-signal-processor

A comprehensive signal processing library for scientific instrumentation, plasma diagnostics, and experimental physics data analysis. Built from real-world experience in fusion energy research.

Version 0.2.0: a faster FFT that handles any input length, a working inverse FFT, Hilbert envelopes, Welch power spectral density, and correctness fixes checked against numpy and scipy. See the changelog.

Known Limitations (v0.2.0)

• Savitzky-Golay filter: Only supports window size 5 (throws for other sizes)
• Peak prominence: Uses a simplified approximation, not true topographic prominence
• FFT scaling: fft() and computeFFT() return raw, unnormalized DFT values (for an amplitude spectrum, divide by N/2, or by N at DC and Nyquist)
• Arbitrary-length FFTs (Bluestein's algorithm) are several times slower than power-of-2 lengths; pad with zeroPad() when the exact length doesn't matter
• Median filter: O(n·w·log w) complexity, which may be slow for large windows or real-time use

For production-critical applications, validate results against established tools (scipy, MATLAB, etc.).

Features

• Noise Analysis & Reduction: SNR calculation, outlier detection, wavelet denoising, baseline correction
• Digital Filtering: Moving average, Savitzky-Golay, Butterworth, median filters, and more
• Frequency Analysis: FFT and inverse FFT of any length, periodogram and Welch PSD, analytic signal (Hilbert transform), spectrograms, dominant frequency detection
• Signal Analysis: Peak detection, cross-correlation, phase shift, zero-crossing detection, Hilbert and RMS envelopes
• Statistical Tools: Comprehensive statistics, interpolation, window functions
• TypeScript: Full type definitions for excellent IDE support

Installation

npm install plasma-signal-processor

Or with yarn:

yarn add plasma-signal-processor

Quick Start

import {
  NoiseAnalysis,
  DigitalFilters,
  FrequencyAnalysis,
  SignalAnalysis,
} from 'plasma-signal-processor';

// Your noisy signal data
const signal = [/* your data */];
const sampleRate = 10000; // Hz

// 1. Analyze noise characteristics
const noise = NoiseAnalysis.estimateNoise(signal, { method: 'mad' });
console.log('Noise level:', noise.stdDev);

// 2. Remove baseline drift
const corrected = NoiseAnalysis.removeBaseline(signal, 'median');

// 3. Apply smoothing filter
const filtered = DigitalFilters.savitzkyGolay(corrected, 5);

// 4. Calculate SNR
const snr = NoiseAnalysis.calculateSNR(filtered, noise.stdDev);
console.log('SNR:', snr, 'dB');

// 5. Find peaks
const peaks = SignalAnalysis.findPeaks(filtered, {
  minHeight: 10,
  minDistance: 50,
});

// 6. Frequency analysis
const psd = FrequencyAnalysis.welch(filtered, sampleRate, { segmentLength: 1024 });
const dominant = FrequencyAnalysis.findDominantFrequency(filtered, sampleRate);
console.log('Dominant frequency:', dominant.frequency, 'Hz');

Core Modules

1. Noise Analysis

Handle noisy scientific data with robust noise estimation and reduction techniques.

import { NoiseAnalysis } from 'plasma-signal-processor';

// Estimate noise using MAD (robust to outliers)
const noise = NoiseAnalysis.estimateNoise(data, { method: 'mad' });

// Calculate Signal-to-Noise Ratio
const snr = NoiseAnalysis.calculateSNR(signal, noiseLevel);

// Detect and remove outliers
const outliers = NoiseAnalysis.detectOutliers(data, 1.5);
const cleaned = NoiseAnalysis.removeOutliers(data, 1.5);

// Noise gating (threshold below which signal is zeroed)
const gated = NoiseAnalysis.noiseGate(data, 0.1, {
  relative: true,
  smoothTransition: true,
});

// Baseline correction
const baselineCorrected = NoiseAnalysis.removeBaseline(data, 'rolling', 100);

// Wavelet denoising
const denoised = NoiseAnalysis.waveletDenoise(data, 1.0);

2. Digital Filters

Various filtering techniques optimized for scientific instrumentation.

import { DigitalFilters } from 'plasma-signal-processor';

// Moving average (simple smoothing)
const smoothed = DigitalFilters.movingAverage(data, 10);

// Savitzky-Golay filter (preserves peak shapes)
const sgFiltered = DigitalFilters.savitzkyGolay(data, 5);

// Exponential moving average (good for real-time)
const ema = DigitalFilters.exponentialMovingAverage(data, 0.1);

// Lowpass filter (remove high-frequency noise)
const lowpassed = DigitalFilters.lowpass(data, {
  cutoffFrequency: 1000,
  sampleRate: 10000,
});

// Highpass filter (remove DC offset and drift)
const highpassed = DigitalFilters.highpass(data, {
  cutoffFrequency: 100,
  sampleRate: 10000,
});

// Bandpass filter
const bandpassed = DigitalFilters.bandpass(data, {
  lowCutoff: 100,
  highCutoff: 1000,
  sampleRate: 10000,
});

// Median filter (excellent for spike removal)
const medianFiltered = DigitalFilters.median(data, 5);

// Butterworth filter (sharper rolloff)
const butterworth = DigitalFilters.butterworthLowpass(data, {
  cutoffFrequency: 1000,
  sampleRate: 10000,
});

// Multi-pass filtering for sharper response
const sharpFiltered = DigitalFilters.multipass(
  data,
  (d) => DigitalFilters.lowpass(d, { cutoffFrequency: 1000, sampleRate: 10000 }),
  3 // number of passes
);

3. Frequency Analysis (FFT)

Comprehensive Fourier analysis tools. See How It Works for the math.

import { FrequencyAnalysis } from 'plasma-signal-processor';

// Forward and inverse FFT of any length (same convention as numpy.fft)
const spectrum = FrequencyAnalysis.fft(data); // Complex[]
const restored = FrequencyAnalysis.ifft(spectrum).map((c) => c.real);

// FFT with full spectral information (zero-padded to a power of 2)
const fft = FrequencyAnalysis.computeFFT(data, sampleRate);
// Returns: { frequencies, magnitudes, phases, powerSpectrum, real, imaginary }

// The window function is applied before zero-padding and is configurable
// (default 'hanning'; use 'rectangular' to disable windowing)
const raw = FrequencyAnalysis.computeFFT(data, sampleRate, 'rectangular');

// Power spectral density (units²/Hz) from a single periodogram...
const psd = FrequencyAnalysis.powerSpectralDensity(data, sampleRate);

// ...or with Welch's method: averaged overlapping segments, far less variance
const welch = FrequencyAnalysis.welch(data, sampleRate, {
  segmentLength: 1024, // default 256
  overlap: 512, // default segmentLength / 2
  windowFunction: 'hanning',
  detrend: 'constant', // subtract each segment's mean
});
// Returns: { frequencies, psd, segmentCount }

// Analytic signal x + i·H{x}: instantaneous amplitude and phase
const { real, imag } = FrequencyAnalysis.analyticSignal(data);
const amplitude = real.map((re, i) => Math.hypot(re, imag[i]));
const phase = real.map((re, i) => Math.atan2(imag[i], re));

// Spectrogram (time-frequency analysis)
const spectrogram = FrequencyAnalysis.spectrogram(data, sampleRate, {
  windowSize: 256,
  hopSize: 128,
  windowFunction: 'hanning',
});

// Find dominant frequency
const dominant = FrequencyAnalysis.findDominantFrequency(data, sampleRate);
console.log(dominant.frequency, dominant.magnitude);

// Total harmonic distortion
const thd = FrequencyAnalysis.totalHarmonicDistortion(
  data,
  sampleRate,
  fundamentalFrequency
);

4. Signal Analysis

Advanced signal processing and feature extraction.

import { SignalAnalysis } from 'plasma-signal-processor';

// Peak detection with prominence
const peaks = SignalAnalysis.findPeaks(data, {
  minHeight: 10,
  minDistance: 50,
  threshold: 2, // minimum prominence
});

// Cross-correlation between signals
const correlation = SignalAnalysis.crossCorrelation(signal1, signal2, {
  normalize: true,
  maxLag: 100,
  scaling: 'biased', // default; 'unbiased' divides each lag by its overlap
});

// Auto-correlation
const autoCorr = SignalAnalysis.autoCorrelation(signal);

// Zero crossing detection
const crossings = SignalAnalysis.findZeroCrossings(data, 'rising');

// Signal envelope: Hilbert (instantaneous amplitude) or moving RMS
const hilbertEnvelope = SignalAnalysis.envelope(data, 'hilbert');
const rmsEnvelope = SignalAnalysis.envelope(data, 'rms', 50);

// Derivative and integral
const derivative = SignalAnalysis.derivative(data, sampleRate);
const integral = SignalAnalysis.integrate(data, sampleRate);

// Edge detection
const edges = SignalAnalysis.detectEdges(data, 0.5, 'both');

// Phase shift in radians (positive when signal2 lags signal1). Pass
// `frequency` as well when it is known, e.g. a modulation frequency
const phaseShift = SignalAnalysis.phaseShift(signal1, signal2, { sampleRate });

5. Statistics & Utilities

import { Statistics, WindowFunctions, Interpolation } from 'plasma-signal-processor';

// Comprehensive statistics
const stats = Statistics.summary(data);
// Returns: { mean, median, stdDev, variance, min, max, range, skewness, kurtosis }

// Individual metrics
const mean = Statistics.mean(data);
const median = Statistics.median(data);
const stdDev = Statistics.stdDev(data);
const rms = Statistics.rms(data);
const percentile95 = Statistics.percentile(data, 95);

// Window functions
const windowed = WindowFunctions.applyWindow(data, 'hanning');

// Interpolation
const interpolated = Interpolation.linear(xValues, yValues, newX);

How It Works

The transforms below are unit-tested against numpy and scipy reference values, agreeing to about 1e-12.

Fourier transforms

fft() and ifft() use the same convention as numpy.fft:

Power-of-2 lengths run an iterative, in-place radix-2 FFT on Float64Arrays with cached sine and cosine tables. Other lengths use Bluestein's algorithm, which rewrites the transform as a chirp-modulated convolution and evaluates it with power-of-2 FFTs of length M ≥ 2N − 1:

Both paths are O(N log N), and ifft(fft(x)) returns x for any length.

| N | v0.1.0 recursive | v0.2.0 iterative | Speedup | v0.2.0 N + 1 (Bluestein) | |---:|---:|---:|---:|---:| | 1,024 | 0.31 ms | 0.10 ms | 3.0× | 0.95 ms | | 4,096 | 1.68 ms | 0.17 ms | 10.0× | 2.41 ms | | 16,384 | 8.54 ms | 0.89 ms | 9.6× | 11.84 ms | | 65,536 | 45.90 ms | 4.35 ms | 10.5× | 66.24 ms | | 262,144 | 290.03 ms | 45.66 ms | 6.4× | 437.04 ms | | 1,048,576 | 1855.67 ms | 534.91 ms | 3.5× | 3518.97 ms |

Timings come from npm run docs:figures and depend on the machine. For large inputs, a good share of the time goes to building the Complex[] output rather than to the transform itself.

Welch power spectral density

welch() splits the record into L segments of M samples, advancing D = M − overlap samples each time. Each segment has its mean removed and is windowed and transformed, and the periodograms are averaged:

The DC bin, and the Nyquist bin when Nfft is even, are not doubled. Averaging trades frequency resolution (fs/M) for variance, which lets weak tones rise out of the noise:

Analytic signal and envelope

analyticSignal() removes the negative frequencies and doubles the positive ones, matching scipy.signal.hilbert. Its magnitude is the instantaneous amplitude that envelope(x, 'hilbert') returns:

The envelope is exact when the record holds whole cycles; otherwise expect some ripple near the ends.

Real-World Use Cases

Thomson Scattering Diagnostics

// Process Thomson scattering data from plasma diagnostic
const raw = loadDiagnosticData();

// Remove baseline
const corrected = NoiseAnalysis.removeBaseline(raw, 'median');

// Smooth while preserving peak shape
const filtered = DigitalFilters.savitzkyGolay(corrected, 5);

// Find scattering peaks
const peaks = SignalAnalysis.findPeaks(filtered, {
  minHeight: threshold,
  minDistance: 100,
});

// Calculate electron temperature from peak width
const temperature = calculateTemperature(peaks);

XUV Spectroscopy

// Analyze X-ray/UV spectroscopy data
const spectrum = loadXUVData();

// Noise reduction
const denoised = NoiseAnalysis.waveletDenoise(spectrum, 2.0);

// Baseline subtraction (rolling window for varying background)
const baselineRemoved = NoiseAnalysis.removeBaseline(denoised, 'rolling', 200);

// Peak detection for spectral lines
const spectralLines = SignalAnalysis.findPeaks(baselineRemoved, {
  minHeight: 100,
  threshold: 5, // prominence
});

Time-Series Analysis

// Analyze plasma instabilities over time
const timeSeries = loadTimeSeriesData();
const sampleRate = 100000; // 100 kHz

// Compute spectrogram to see frequency evolution
const spectrogram = FrequencyAnalysis.spectrogram(timeSeries, sampleRate, {
  windowSize: 1024,
  hopSize: 256,
  windowFunction: 'hamming',
});

// Low-variance spectrum of the whole shot
const spectrum = FrequencyAnalysis.welch(timeSeries, sampleRate, {
  segmentLength: 4096,
});

// Mode amplitude over time from a magnetic probe signal
const modeAmplitude = SignalAnalysis.envelope(timeSeries, 'hilbert');

// Detect mode transitions
const edges = SignalAnalysis.detectEdges(timeSeries, 0.1, 'both');

API Reference

Full API documentation is available in the TypeScript definitions. Your IDE will provide autocomplete and inline documentation.

Examples

See the examples/ directory for complete working examples:

• thomson-scattering.ts - Processing Thomson scattering diagnostic data
• real-time-processing.ts - Real-time signal processing patterns
• spectrogram-analysis.ts - Time-frequency analysis of plasma instabilities

Development

# Install dependencies
npm install

# Build
npm run build

# Run tests
npm test

# Run tests in watch mode
npm run test:watch

# Generate coverage report
npm run test:coverage

# Lint
npm run lint

# Format code
npm run format

# Regenerate the README figures and equations
# (needs Python with numpy and matplotlib, and LaTeX with dvisvgm)
npm run docs:figures

Roadmap

Planned for v0.3.0:

• Zero-phase filtering (filtfilt) and general Butterworth design (highpass, bandpass, higher orders)
• Savitzky-Golay filters for any window size and polynomial order
• Polynomial fitting for calibration curves, and polynomial or asymmetric-least-squares baseline removal
• FFT-based cross-correlation in O(N log N)
• Normalized amplitude and power spectrum helpers
• Float64Array inputs and an ESM build
• True topographic peak prominence and peak widths
• Faster median filter and stateful real-time filter classes

Background

This library was developed from scientific instrumentation workflows, with particular focus on plasma diagnostics for fusion energy research. It incorporates techniques commonly used in:

• Plasma diagnostics (Thomson scattering, XUV spectroscopy, polarimetry)
• High-noise experimental environments
• Real-time data acquisition systems
• Scientific instrumentation

Note: While the algorithms are grounded in physics research, this is a general-purpose DSP library suitable for any scientific signal processing application.

Contributing

Contributions are welcome! Please feel free to submit a Pull Request.

License

Apache License 2.0 – see the LICENSE file for details.

Author

Developed from scientific instrumentation workflows in experimental physics research.

Acknowledgments

Built from techniques used in plasma physics research and diagnostic systems for fusion energy experiments.