diff --git a/src/mcltspice/output_format.py b/src/mcltspice/output_format.py
new file mode 100644
index 0000000..314d772
--- /dev/null
+++ b/src/mcltspice/output_format.py
@@ -0,0 +1,81 @@
+"""Output-formatting helpers for compact JSON responses.
+
+Full float64 precision (~17 digits) bloats JSON output; these round to an
+appropriate precision and downsample long arrays before serialization.
+"""
+
+import math
+
+import numpy as np
+
+
+def round_sig(x: float, sig: int = 6) -> float:
+    """Round to N significant figures for compact JSON output."""
+    if x == 0 or not math.isfinite(x):
+        return float(x)
+    digits = sig - 1 - int(math.floor(math.log10(abs(x))))
+    return round(x, digits)
+
+
+def compact_list(
+    values, decimal_places: int | None = None, sig_figs: int | None = None
+) -> list[float]:
+    """Round a list of floats for compact JSON output.
+
+    Full float64 precision (~17 digits) causes massive JSON bloat.
+    Rounding to appropriate precision cuts output size 50-60%.
+    """
+    if decimal_places is not None:
+        return [round(float(v), decimal_places) for v in values]
+    if sig_figs is not None:
+        return [round_sig(float(v), sig_figs) for v in values]
+    return [float(v) for v in values]
+
+
+def downsample_indices(
+    n_total: int,
+    max_points: int,
+    signals: list[np.ndarray] | None = None,
+) -> np.ndarray:
+    """Compute indices for downsampling, preserving peaks in AC data.
+
+    For AC data (when complex signal arrays are provided), uses
+    peak-preserving bucket selection: divides the data into max_points
+    buckets and keeps the point in each bucket whose magnitude deviates
+    most from the bucket mean. This ensures narrow resonance peaks
+    (e.g. high-Q bandpass filters) aren't lost by blind stride sampling.
+
+    For transient data (no signals provided), uses simple stride.
+    """
+    if n_total <= max_points:
+        return np.arange(n_total)
+
+    if signals is None or len(signals) == 0:
+        step = max(1, n_total // max_points)
+        return np.arange(0, n_total, step)
+
+    # Pre-compute magnitude_dB for each signal
+    mag_db_list = []
+    for sig in signals:
+        mag = np.abs(sig)
+        with np.errstate(divide="ignore"):
+            db = np.where(mag > 0, 20 * np.log10(mag), -200.0)
+        mag_db_list.append(db)
+
+    # Bucket-based peak-preserving selection
+    indices = []
+    bucket_size = n_total / max_points
+
+    for i in range(max_points):
+        start = int(i * bucket_size)
+        end = min(int((i + 1) * bucket_size), n_total)
+        if start >= end:
+            continue
+        # Combined importance: sum of |deviation from bucket mean| across signals
+        combined = np.zeros(end - start)
+        for db in mag_db_list:
+            bucket = db[start:end]
+            combined += np.abs(bucket - np.mean(bucket))
+        indices.append(start + int(np.argmax(combined)))
+
+    return np.array(indices)
diff --git a/src/mcltspice/server.py b/src/mcltspice/server.py
index 426998d..616142e 100644
--- a/src/mcltspice/server.py
+++ b/src/mcltspice/server.py
@@ -49,6 +49,7 @@ from .optimizer import (
     format_engineering,
     optimize_component_values,
 )
+from .output_format import compact_list
 from .plotting import build_waveform_svg
 from .power_analysis import compute_efficiency, compute_power_metrics
 from .raw_parser import parse_raw_file
@@ -77,13 +78,8 @@ from .tuning import (
 from .waveform_expr import WaveformCalculator
 from .waveform_math import (
     compute_bandwidth,
-    compute_fft,
-    compute_peak_to_peak,
-    compute_rise_time,
-    compute_rms,
-    compute_settling_time,
-    compute_thd,
 )
+from .waveform_query import analyze_signal, extract_waveform
 
 mcp = FastMCP(
     name="mcltspice",
@@ -125,83 +121,6 @@ mcp = FastMCP(
 )
 
 
-# ============================================================================
-# OUTPUT HELPERS
-# ============================================================================
-
-
-def _round_sig(x: float, sig: int = 6) -> float:
-    """Round to N significant figures for compact JSON output."""
-    if x == 0 or not math.isfinite(x):
-        return float(x)
-    digits = sig - 1 - int(math.floor(math.log10(abs(x))))
-    return round(x, digits)
-
-
-def _compact_list(
-    values, decimal_places: int | None = None, sig_figs: int | None = None
-) -> list[float]:
-    """Round a list of floats for compact JSON output.
-
-    Full float64 precision (~17 digits) causes massive JSON bloat.
-    Rounding to appropriate precision cuts output size 50-60%.
-    """
-    if decimal_places is not None:
-        return [round(float(v), decimal_places) for v in values]
-    if sig_figs is not None:
-        return [_round_sig(float(v), sig_figs) for v in values]
-    return [float(v) for v in values]
-
-
-def _downsample_indices(
-    n_total: int,
-    max_points: int,
-    signals: list[np.ndarray] | None = None,
-) -> np.ndarray:
-    """Compute indices for downsampling, preserving peaks in AC data.
-
-    For AC data (when complex signal arrays are provided), uses
-    peak-preserving bucket selection: divides the data into max_points
-    buckets and keeps the point in each bucket whose magnitude deviates
-    most from the bucket mean. This ensures narrow resonance peaks
-    (e.g. high-Q bandpass filters) aren't lost by blind stride sampling.
-
-    For transient data (no signals provided), uses simple stride.
-    """
-    if n_total <= max_points:
-        return np.arange(n_total)
-
-    if signals is None or len(signals) == 0:
-        step = max(1, n_total // max_points)
-        return np.arange(0, n_total, step)
-
-    # Pre-compute magnitude_dB for each signal
-    mag_db_list = []
-    for sig in signals:
-        mag = np.abs(sig)
-        with np.errstate(divide="ignore"):
-            db = np.where(mag > 0, 20 * np.log10(mag), -200.0)
-        mag_db_list.append(db)
-
-    # Bucket-based peak-preserving selection
-    indices = []
-    bucket_size = n_total / max_points
-
-    for i in range(max_points):
-        start = int(i * bucket_size)
-        end = min(int((i + 1) * bucket_size), n_total)
-        if start >= end:
-            continue
-        # Combined importance: sum of |deviation from bucket mean| across signals
-        combined = np.zeros(end - start)
-        for db in mag_db_list:
-            bucket = db[start:end]
-            combined += np.abs(bucket - np.mean(bucket))
-        indices.append(start + int(np.argmax(combined)))
-
-    return np.array(indices)
-
-
 # ============================================================================
 # SIMULATION TOOLS
 # ============================================================================
@@ -329,89 +248,14 @@ def get_waveform(
         x_max: Maximum x-axis value (frequency Hz or time s) to include
     """
     raw = parse_raw_file(raw_file_path)
-
-    # Extract specific run if requested
-    if run is not None:
-        if not raw.is_stepped:
-            return {"error": "Not a stepped simulation - no multiple runs available"}
-        if run < 1 or run > raw.n_runs:
-            return {"error": f"Run {run} out of range (1..{raw.n_runs})"}
-        raw = raw.get_run_data(run)
-
-    x_axis = raw.get_time()
-    x_name = "time"
-    if x_axis is None:
-        x_axis = raw.get_frequency()
-        x_name = "frequency"
-
-    if x_axis is None:
-        return {"error": "No time or frequency axis found in raw file"}
-
-    x_real = x_axis.real if np.iscomplexobj(x_axis) else x_axis
-    is_complex = np.iscomplexobj(raw.data)
-
-    # Apply x-axis range filter
-    mask = np.ones(len(x_real), dtype=bool)
-    if x_min is not None:
-        mask &= x_real >= x_min
-    if x_max is not None:
-        mask &= x_real <= x_max
-    filtered_indices = np.where(mask)[0]
-
-    if len(filtered_indices) == 0:
-        return {
-            "error": f"No data points in {x_name} range [{x_min}, {x_max}]",
-            "total_points": len(x_real),
-            f"{x_name}_range": [float(x_real[0]), float(x_real[-1])],
-        }
-
-    # Collect signal data for peak-preserving downsampling (AC only)
-    ac_signals = []
-    if is_complex:
-        for name in signal_names:
-            data = raw.get_variable(name)
-            if data is not None:
-                ac_signals.append(data[filtered_indices])
-
-    # Downsample: peak-preserving for AC, stride for transient
-    ds_indices = _downsample_indices(
-        len(filtered_indices),
-        max_points,
-        signals=ac_signals if ac_signals else None,
+    return extract_waveform(
+        raw,
+        signal_names,
+        max_points=max_points,
+        run=run,
+        x_min=x_min,
+        x_max=x_max,
     )
-    sample_indices = filtered_indices[ds_indices]
-
-    result = {
-        "x_axis_name": x_name,
-        "x_axis_data": _compact_list(x_real[sample_indices], sig_figs=6),
-        "signals": {},
-        "total_points": len(x_real),
-        "returned_points": len(sample_indices),
-        "is_stepped": raw.is_stepped,
-        "n_runs": raw.n_runs,
-    }
-
-    for name in signal_names:
-        data = raw.get_variable(name)
-        if data is not None:
-            sampled = data[sample_indices]
-            if np.iscomplexobj(sampled):
-                result["signals"][name] = {
-                    "magnitude_db": _compact_list(
-                        [20 * math.log10(abs(x)) if abs(x) > 0 else -200 for x in sampled],
-                        decimal_places=2,
-                    ),
-                    "phase_degrees": _compact_list(
-                        [math.degrees(math.atan2(x.imag, x.real)) for x in sampled],
-                        decimal_places=2,
-                    ),
-                }
-            else:
-                result["signals"][name] = {
-                    "values": _compact_list(sampled, sig_figs=6)
-                }
-
-    return result
 
 
 @mcp.tool()
@@ -484,66 +328,17 @@ def analyze_waveform(
         thd_n_harmonics: Number of harmonics for THD calculation
     """
     raw = parse_raw_file(raw_file_path)
-
-    time = raw.get_time()
-    signal = raw.get_variable(signal_name)
-
-    if signal is None:
-        return {
-            "error": f"Signal '{signal_name}' not found. Available: "
-            f"{[v.name for v in raw.variables]}"
-        }
-
-    # Use real parts for time-domain analysis
-    if np.iscomplexobj(time):
-        time = time.real
-    if np.iscomplexobj(signal):
-        signal = np.abs(signal)
-
-    results = {"signal": signal_name}
-
-    for analysis in analyses:
-        if analysis == "rms":
-            results["rms"] = compute_rms(signal)
-
-        elif analysis == "peak_to_peak":
-            results["peak_to_peak"] = compute_peak_to_peak(signal)
-
-        elif analysis == "settling_time":
-            if time is not None:
-                results["settling_time"] = compute_settling_time(
-                    time,
-                    signal,
-                    final_value=settling_final_value,
-                    tolerance_percent=settling_tolerance_pct,
-                )
-
-        elif analysis == "rise_time":
-            if time is not None:
-                results["rise_time"] = compute_rise_time(
-                    time,
-                    signal,
-                    low_pct=rise_low_pct,
-                    high_pct=rise_high_pct,
-                )
-
-        elif analysis == "fft":
-            if time is not None:
-                results["fft"] = compute_fft(
-                    time,
-                    signal,
-                    max_harmonics=fft_max_harmonics,
-                )
-
-        elif analysis == "thd":
-            if time is not None:
-                results["thd"] = compute_thd(
-                    time,
-                    signal,
-                    n_harmonics=thd_n_harmonics,
-                )
-
-    return results
+    return analyze_signal(
+        raw,
+        signal_name,
+        analyses,
+        settling_tolerance_pct=settling_tolerance_pct,
+        settling_final_value=settling_final_value,
+        rise_low_pct=rise_low_pct,
+        rise_high_pct=rise_high_pct,
+        fft_max_harmonics=fft_max_harmonics,
+        thd_n_harmonics=thd_n_harmonics,
+    )
 
 
 @mcp.tool()
@@ -1052,9 +847,9 @@ def evaluate_waveform_expression(
     if x_axis is not None:
         x_real = x_axis.real if np.iscomplexobj(x_axis) else x_axis
         response["x_axis_name"] = x_name
-        response["x_axis_data"] = _compact_list(x_real[sample_indices], sig_figs=6)
+        response["x_axis_data"] = compact_list(x_real[sample_indices], sig_figs=6)
 
-    response["values"] = _compact_list(result[sample_indices], sig_figs=6)
+    response["values"] = compact_list(result[sample_indices], sig_figs=6)
     response["available_signals"] = calc.available_signals()
 
     return response
diff --git a/src/mcltspice/waveform_query.py b/src/mcltspice/waveform_query.py
new file mode 100644
index 0000000..aab81c7
--- /dev/null
+++ b/src/mcltspice/waveform_query.py
@@ -0,0 +1,197 @@
+"""Read-side queries over a parsed RawFile: data extraction and analysis.
+
+Both functions take an already-parsed RawFile and return a JSON-ready dict, so
+they can be unit-tested without an LTspice run. The get_waveform / analyze_waveform
+MCP tools wrap these with the .raw file read.
+"""
+
+import math
+
+import numpy as np
+
+from .output_format import compact_list, downsample_indices
+from .raw_parser import RawFile
+from .waveform_math import (
+    compute_fft,
+    compute_peak_to_peak,
+    compute_rise_time,
+    compute_rms,
+    compute_settling_time,
+    compute_thd,
+)
+
+
+def extract_waveform(
+    raw: RawFile,
+    signal_names: list[str],
+    max_points: int = 1000,
+    run: int | None = None,
+    x_min: float | None = None,
+    x_max: float | None = None,
+) -> dict:
+    """Extract sampled signal data from a RawFile.
+
+    Returns the x-axis plus per-signal values (magnitude_db/phase for AC,
+    values for transient), or an {"error": ...} dict for a bad run number or an
+    empty x-range. AC data is downsampled peak-preservingly; transient by stride.
+    """
+    # Extract specific run if requested
+    if run is not None:
+        if not raw.is_stepped:
+            return {"error": "Not a stepped simulation - no multiple runs available"}
+        if run < 1 or run > raw.n_runs:
+            return {"error": f"Run {run} out of range (1..{raw.n_runs})"}
+        raw = raw.get_run_data(run)
+
+    x_axis = raw.get_time()
+    x_name = "time"
+    if x_axis is None:
+        x_axis = raw.get_frequency()
+        x_name = "frequency"
+
+    if x_axis is None:
+        return {"error": "No time or frequency axis found in raw file"}
+
+    x_real = x_axis.real if np.iscomplexobj(x_axis) else x_axis
+    is_complex = np.iscomplexobj(raw.data)
+
+    # Apply x-axis range filter
+    mask = np.ones(len(x_real), dtype=bool)
+    if x_min is not None:
+        mask &= x_real >= x_min
+    if x_max is not None:
+        mask &= x_real <= x_max
+    filtered_indices = np.where(mask)[0]
+
+    if len(filtered_indices) == 0:
+        return {
+            "error": f"No data points in {x_name} range [{x_min}, {x_max}]",
+            "total_points": len(x_real),
+            f"{x_name}_range": [float(x_real[0]), float(x_real[-1])],
+        }
+
+    # Collect signal data for peak-preserving downsampling (AC only)
+    ac_signals = []
+    if is_complex:
+        for name in signal_names:
+            data = raw.get_variable(name)
+            if data is not None:
+                ac_signals.append(data[filtered_indices])
+
+    # Downsample: peak-preserving for AC, stride for transient
+    ds_indices = downsample_indices(
+        len(filtered_indices),
+        max_points,
+        signals=ac_signals if ac_signals else None,
+    )
+    sample_indices = filtered_indices[ds_indices]
+
+    result = {
+        "x_axis_name": x_name,
+        "x_axis_data": compact_list(x_real[sample_indices], sig_figs=6),
+        "signals": {},
+        "total_points": len(x_real),
+        "returned_points": len(sample_indices),
+        "is_stepped": raw.is_stepped,
+        "n_runs": raw.n_runs,
+    }
+
+    for name in signal_names:
+        data = raw.get_variable(name)
+        if data is not None:
+            sampled = data[sample_indices]
+            if np.iscomplexobj(sampled):
+                result["signals"][name] = {
+                    "magnitude_db": compact_list(
+                        [20 * math.log10(abs(x)) if abs(x) > 0 else -200 for x in sampled],
+                        decimal_places=2,
+                    ),
+                    "phase_degrees": compact_list(
+                        [math.degrees(math.atan2(x.imag, x.real)) for x in sampled],
+                        decimal_places=2,
+                    ),
+                }
+            else:
+                result["signals"][name] = {
+                    "values": compact_list(sampled, sig_figs=6)
+                }
+
+    return result
+
+
+def analyze_signal(
+    raw: RawFile,
+    signal_name: str,
+    analyses: list[str],
+    settling_tolerance_pct: float = 2.0,
+    settling_final_value: float | None = None,
+    rise_low_pct: float = 10.0,
+    rise_high_pct: float = 90.0,
+    fft_max_harmonics: int = 50,
+    thd_n_harmonics: int = 10,
+) -> dict:
+    """Run one or more analyses on a signal from a RawFile.
+
+    Supported analyses: rms, peak_to_peak, settling_time, rise_time, fft, thd.
+    Returns a dict keyed by analysis name, or {"error": ...} if the signal is
+    absent. Frequency-domain analyses are skipped when there's no time axis.
+    """
+    time = raw.get_time()
+    signal = raw.get_variable(signal_name)
+
+    if signal is None:
+        return {
+            "error": f"Signal '{signal_name}' not found. Available: "
+            f"{[v.name for v in raw.variables]}"
+        }
+
+    # Use real parts for time-domain analysis
+    if np.iscomplexobj(time):
+        time = time.real
+    if np.iscomplexobj(signal):
+        signal = np.abs(signal)
+
+    results = {"signal": signal_name}
+
+    for analysis in analyses:
+        if analysis == "rms":
+            results["rms"] = compute_rms(signal)
+
+        elif analysis == "peak_to_peak":
+            results["peak_to_peak"] = compute_peak_to_peak(signal)
+
+        elif analysis == "settling_time":
+            if time is not None:
+                results["settling_time"] = compute_settling_time(
+                    time,
+                    signal,
+                    final_value=settling_final_value,
+                    tolerance_percent=settling_tolerance_pct,
+                )
+
+        elif analysis == "rise_time":
+            if time is not None:
+                results["rise_time"] = compute_rise_time(
+                    time,
+                    signal,
+                    low_pct=rise_low_pct,
+                    high_pct=rise_high_pct,
+                )
+
+        elif analysis == "fft":
+            if time is not None:
+                results["fft"] = compute_fft(
+                    time,
+                    signal,
+                    max_harmonics=fft_max_harmonics,
+                )
+
+        elif analysis == "thd":
+            if time is not None:
+                results["thd"] = compute_thd(
+                    time,
+                    signal,
+                    n_harmonics=thd_n_harmonics,
+                )
+
+    return results