diff --git a/README.rst b/README.rst
index 683b209..9a26a98 100644
--- a/README.rst
+++ b/README.rst
@@ -172,6 +172,7 @@ read off computational-mechanics quantities:
 
    eps = EpsilonMachine.from_hmm(gm)          # minimize an HMM presentation
    # eps = EpsilonMachine.from_sequence(data, method="cssr", Lmax=4)  # infer
+   # eps = EpsilonMachine.from_sequence(data, method="spectral", prefix_length=3, rank=2)
 
    eps.statistical_complexity()   # 0.9183 bits  (C_mu)
    eps.entropy_rate()             # 0.6667 bits/symbol  (h_mu)
diff --git a/docs/generators/epsilon_inference.rst b/docs/generators/epsilon_inference.rst
index 6af90eb..8586acd 100644
--- a/docs/generators/epsilon_inference.rst
+++ b/docs/generators/epsilon_inference.rst
@@ -9,10 +9,11 @@ Sample-based reconstruction of ε-machines from observed symbol sequences.
 This complements the **oracle** path :func:`~sofic.generators.epsilon_construction.build_epsilon_machine`,
 which merges probabilistically equivalent states in a *given* generator.
 
-Two algorithms are implemented:
+Three algorithms are implemented:
 
 * **CSSR** (Causal-State Splitting Reconstruction) :cite:`Shalizi2002,Shalizi2004`
 * **Subtree merging** (Crutchfield--Young batch reconstruction) :cite:`CrutchfieldYoung1989,Crutchfield1994`
+* **Spectral** (Hankel SVD, then mixed-state extraction) :cite:`Balle2014,Hsu2012,Ellison2009`
 
 Quick start
 ===========
@@ -55,16 +56,50 @@ numerical tolerance for finite-sample estimates.
 
 .. autofunction:: subtree_merge
 
+Spectral reconstruction
+=======================
+
+Spectral reconstruction learns a weighted finite automaton from the Hankel matrix
+of block probabilities :cite:`Balle2014,Hsu2012`, then extracts causal states as
+mixed states of the learned operators :cite:`Ellison2009`. When the operators
+are non-negative in the learned basis this is a Mealy projection followed by
+:meth:`~sofic.generators.epsilon_machine.EpsilonMachine.from_hmm`; signed
+operators use mixed-state enumeration of the observable operators (not a
+clustering heuristic). The same extraction is
+:func:`~sofic.inference.spectral.project_to_epsilon_machine`.
+
+Pass ``rank`` when the model order is known (two for the golden mean and even
+process). Otherwise the Hankel singular-value gap selects the rank.
+
+.. code-block:: python
+
+   inferred = EpsilonMachine.from_sequence(
+       observations, method="spectral", prefix_length=3, rank=2
+   )
+
+.. autofunction:: spectral
+
 Unified entry point
 ===================
 
 Use :meth:`~sofic.generators.epsilon_machine.EpsilonMachine.from_sequence` to
-dispatch to CSSR or subtree merging (see :doc:`epsilon_machine`).
+dispatch to CSSR, subtree merging, or spectral reconstruction
+(see :doc:`epsilon_machine`).
+
+Related inference methods
+=========================
+
+**transCSSR** — input/output ε-transducers; see :doc:`epsilon_transducer_inference`.
 
-Related inference methods (not yet implemented)
-===============================================
+**Bayesian structural inference** — conjugate Dirichlet–multinomial evidence over
+candidate unifilar topologies :cite:`Strelioff2014`; see :doc:`../inference/epsilon`.
+This is not a Gibbs clustering heuristic over history labels.
 
-**transCSSR** — input/output ε-transducers (Darmon, 2014).
+Subtree merging is the literature reconstruction by morph clustering; a separate
+agglomerative "causal-state merging" procedure is not provided. k-means on
+history-morph embeddings (sometimes labelled "neural state discovery" despite
+involving no neural network) is also not provided — mixed-state extraction is
+the causal-state construction used after spectral learning.
 
 **VLMC / context algorithm** — sparse Markov trees; not causally minimal in general.
 
@@ -76,4 +111,5 @@ Related inference methods (not yet implemented)
 
 **RKHS ε-machines** — continuous-time extension (arXiv:2011.14821).
 
-See also :doc:`epsilon_machine`, :doc:`constructions`, and :doc:`hmm_inference`.
+See also :doc:`epsilon_machine`, :doc:`constructions`, :doc:`hmm_inference`,
+and :doc:`../inference/spectral`.
diff --git a/docs/inference/inference.rst b/docs/inference/inference.rst
index 791bf5e..b228b4f 100644
--- a/docs/inference/inference.rst
+++ b/docs/inference/inference.rst
@@ -31,7 +31,8 @@ The historical names ``InferMC`` and ``InferEM`` are retained as aliases for
 .. note::
 
    For *non-Bayesian* reconstruction — Causal-State Splitting Reconstruction
-   (CSSR) and subtree merging — see the point-estimate routines in
+   (CSSR), subtree merging, and spectral mixed-state extraction — see the
+   point-estimate routines in
    :doc:`../generators/epsilon_inference`,
    :doc:`../generators/hmm_inference`, and
    :doc:`../generators/stack_inference`.
diff --git a/docs/inference/spectral.rst b/docs/inference/spectral.rst
index d5dc458..478542f 100644
--- a/docs/inference/spectral.rst
+++ b/docs/inference/spectral.rst
@@ -61,6 +61,25 @@ in the learned basis and raises :class:`SpectralInferenceError` otherwise.
 
    machine = project_to_nmachine(model)   # observable-operator generator with a graph
 
+Projection to an ε-machine
+==========================
+
+:func:`project_to_epsilon_machine` extracts the causal presentation: a
+non-negative Mealy projection when one exists in the learned basis, otherwise
+mixed-state enumeration of the observable operators
+:cite:`Ellison2009`. The same path is
+:func:`~sofic.generators.epsilon_inference.spectral` /
+``EpsilonMachine.from_sequence(..., method="spectral")``.
+
+.. code-block:: python
+
+   from sofic.generators.epsilon_inference import spectral
+   from sofic.examples import golden_mean
+
+   process = golden_mean(0.5)
+   eps = spectral(word_probability=process.word_probability, alphabet=(0, 1), prefix_length=3, rank=2)
+   len(list(eps.states()))  # 2
+
 API
 ===
 
@@ -74,4 +93,6 @@ API
 
 .. autofunction:: project_to_mealy
 
+.. autofunction:: project_to_epsilon_machine
+
 .. autoexception:: SpectralInferenceError
diff --git a/sofic/generators/__init__.py b/sofic/generators/__init__.py
index 7d78ad4..07d2d97 100644
--- a/sofic/generators/__init__.py
+++ b/sofic/generators/__init__.py
@@ -16,7 +16,7 @@
     synergistic_information_flow,
     transfer_entropy,
 )
-from sofic.generators.epsilon_inference import cssr, subtree_merge
+from sofic.generators.epsilon_inference import cssr, spectral, subtree_merge
 from sofic.generators.epsilon_machine import EpsilonMachine
 from sofic.generators.epsilon_transducer import EpsilonTransducer
 from sofic.generators.lumping import LumpabilityError, is_lumpable, lump, normalize_partition
@@ -71,6 +71,7 @@
     "is_lumpable",
     "lump",
     "normalize_partition",
+    "spectral",
     "subtree_merge",
     "fit_stack_hmm_mle",
     "learn_stack_hmm_papni",
diff --git a/sofic/generators/epsilon_inference.py b/sofic/generators/epsilon_inference.py
index 0d4ff4e..4cc97f3 100644
--- a/sofic/generators/epsilon_inference.py
+++ b/sofic/generators/epsilon_inference.py
@@ -1,13 +1,15 @@
-"""Sample-based ε-machine reconstruction (CSSR and subtree merging).
+"""Sample-based ε-machine reconstruction (CSSR, subtree merging, and spectral).
 
 CSSR follows Shalizi, Shalizi & Crutchfield (arXiv:cs/0210025). Subtree merging
-follows Crutchfield & Young (PRL 1989; PRE 1994).
+follows Crutchfield & Young (PRL 1989; PRE 1994). Spectral reconstruction learns
+a weighted finite automaton by Hankel SVD :cite:`Balle2014,Hsu2012` and extracts
+causal states as mixed states of the learned operators :cite:`Ellison2009`.
 """
 
 from __future__ import annotations
 
 from collections import Counter, defaultdict
-from collections.abc import Callable, Sequence
+from collections.abc import Callable, Iterable, Sequence
 from dataclasses import dataclass, field
 from typing import Any, ClassVar, Literal
 
@@ -641,3 +643,61 @@ def subtree_merge(
         history: state_id for state_id, histories_in_state in states.items() for history in histories_in_state
     }
     return _counts_to_mealy(states, counts, history_to_state, seq, length=L)
+
+
+def spectral(
+    sequences: Iterable[Any] | None = None,
+    *,
+    word_probability: Callable[[Sequence[Any]], float] | None = None,
+    alphabet: Sequence[Any] | None = None,
+    rank: int | None = None,
+    prefix_length: int = 3,
+    suffix_length: int | None = None,
+    singular_value_threshold: float = 1e-3,
+    min_singular_value: float = 1e-12,
+    max_states: int = 10_000,
+) -> EpsilonMachine:
+    """Reconstruct an ε-machine by spectral learning then mixed-state extraction.
+
+    Learns a weighted finite automaton / observable-operator model from block
+    statistics :cite:`Balle2014,Hsu2012`, then extracts causal states as the
+    mixed states of those operators :cite:`Ellison2009`. When the learned
+    operators are non-negative this is a Mealy projection followed by
+    :meth:`~sofic.generators.epsilon_machine.EpsilonMachine.from_hmm`; signed
+    operators use mixed-state enumeration rather than a clustering heuristic.
+
+    Parameters
+    ----------
+    sequences
+        A single observed realization or an iterable of realizations. Ignored
+        when ``word_probability`` is given.
+    word_probability
+        Optional exact block-probability function ``f(word) -> float``.
+        ``alphabet`` is then required.
+    alphabet
+        Observation alphabet. Inferred from ``sequences`` when omitted.
+    rank
+        Number of latent states. When ``None`` the rank is chosen from the
+        Hankel singular-value spectrum.
+    prefix_length, suffix_length
+        Maximum lengths of the prefix and suffix bases. ``suffix_length``
+        defaults to ``prefix_length``.
+    singular_value_threshold, min_singular_value
+        Cutoffs for automatic rank selection; see
+        :func:`~sofic.inference.spectral.learn_spectral_wfa`.
+    max_states
+        Safety cap on enumerated mixed states.
+    """
+    from sofic.inference.spectral import learn_spectral_wfa, project_to_epsilon_machine
+
+    model = learn_spectral_wfa(
+        sequences,
+        word_probability=word_probability,
+        alphabet=alphabet,
+        rank=rank,
+        prefix_length=prefix_length,
+        suffix_length=suffix_length,
+        singular_value_threshold=singular_value_threshold,
+        min_singular_value=min_singular_value,
+    )
+    return project_to_epsilon_machine(model, max_states=max_states)
diff --git a/sofic/generators/epsilon_machine.py b/sofic/generators/epsilon_machine.py
index b35ac94..dcd0709 100644
--- a/sofic/generators/epsilon_machine.py
+++ b/sofic/generators/epsilon_machine.py
@@ -51,7 +51,7 @@ def from_sequence(
         cls,
         sequence: Sequence[Any],
         *,
-        method: Literal["cssr", "subtree"] = "cssr",
+        method: Literal["cssr", "subtree", "spectral"] = "cssr",
         **kwargs: Any,
     ) -> EpsilonMachine:
         """Reconstruct an ε-machine from an observed symbol sequence.
@@ -61,11 +61,14 @@ def from_sequence(
         sequence
             Observed process realization.
         method
-            ``"cssr"`` for Causal-State Splitting Reconstruction, or
-            ``"subtree"`` for depth-``L`` subtree merging (pass ``L=...``).
+            ``"cssr"`` for Causal-State Splitting Reconstruction,
+            ``"subtree"`` for depth-``L`` subtree merging (pass ``L=...``), or
+            ``"spectral"`` for Hankel-SVD learning followed by mixed-state
+            extraction.
         **kwargs
-            Forwarded to :func:`~sofic.generators.epsilon_inference.cssr` or
-            :func:`~sofic.generators.epsilon_inference.subtree_merge`.
+            Forwarded to :func:`~sofic.generators.epsilon_inference.cssr`,
+            :func:`~sofic.generators.epsilon_inference.subtree_merge`, or
+            :func:`~sofic.generators.epsilon_inference.spectral`.
         """
         if method == "cssr":
             from sofic.generators.epsilon_inference import cssr
@@ -75,6 +78,10 @@ def from_sequence(
             from sofic.generators.epsilon_inference import subtree_merge
 
             return subtree_merge(sequence, **kwargs)
+        if method == "spectral":
+            from sofic.generators.epsilon_inference import spectral
+
+            return spectral(sequence, **kwargs)
         raise ValueError(f"unknown inference method {method!r}")
 
     def copy(self) -> Self:
diff --git a/sofic/inference/__init__.py b/sofic/inference/__init__.py
index 6d7d984..05d7664 100644
--- a/sofic/inference/__init__.py
+++ b/sofic/inference/__init__.py
@@ -17,6 +17,7 @@
     SpectralInferenceError,
     hankel_matrices,
     learn_spectral_wfa,
+    project_to_epsilon_machine,
     project_to_mealy,
     project_to_nmachine,
     spectral_singular_values,
@@ -37,6 +38,7 @@
     "SpectralInferenceError",
     "hankel_matrices",
     "learn_spectral_wfa",
+    "project_to_epsilon_machine",
     "project_to_mealy",
     "project_to_nmachine",
     "spectral_singular_values",
diff --git a/sofic/inference/spectral.py b/sofic/inference/spectral.py
index 1737b92..abf9b23 100644
--- a/sofic/inference/spectral.py
+++ b/sofic/inference/spectral.py
@@ -43,11 +43,15 @@
     "SpectralInferenceError",
     "hankel_matrices",
     "learn_spectral_wfa",
+    "project_to_epsilon_machine",
     "project_to_mealy",
     "project_to_nmachine",
     "spectral_singular_values",
 ]
 
+_BELIEF_DECIMALS = 6
+_MASS_ATOL = 1e-12
+
 
 class SpectralInferenceError(ValueError):
     """Raised when spectral learning or projection cannot proceed."""
@@ -443,3 +447,118 @@ def project_to_mealy(qr: QuasiRealization, *, tol: float = 1e-8, validate: bool
     if validate:
         machine.validate()
     return machine
+
+
+def project_to_epsilon_machine(
+    qr: QuasiRealization,
+    *,
+    tol: float = 1e-8,
+    max_states: int = 10_000,
+) -> Any:
+    """Extract an ε-machine from a learned spectral model.
+
+    When the observable operators admit a non-negative realization in the
+    learned basis, this is :func:`project_to_mealy` followed by
+    :meth:`~sofic.generators.epsilon_machine.EpsilonMachine.from_hmm`. Signed
+    operators are converted by enumerating mixed states of the observable
+    operators (belief updates ``b A_x / (b A_x τ)``) and merging
+    predictively equivalent recurrent states :cite:`Ellison2009`. This is the
+    computational-mechanics extraction, not a clustering heuristic.
+
+    Raises
+    ------
+    SpectralInferenceError
+        If mixed-state enumeration exceeds ``max_states`` or the initial
+        vector is degenerate.
+    """
+    from sofic.generators.epsilon_machine import EpsilonMachine
+
+    try:
+        mealy = project_to_mealy(qr, tol=tol, validate=True)
+    except SpectralInferenceError:
+        mealy = _mealy_from_operator_mixed_states(qr, max_states=max_states)
+    return EpsilonMachine.from_hmm(mealy)
+
+
+def _mealy_from_operator_mixed_states(
+    qr: QuasiRealization,
+    *,
+    max_states: int = 10_000,
+    decimals: int = _BELIEF_DECIMALS,
+) -> Any:
+    """Build a unifilar Mealy HMM whose states are mixed states of ``qr``."""
+    from collections import deque
+
+    from sofic.generators.mealy import MealyHMM
+    from sofic.generators.mixed_state import MixedState
+    from sofic.graph import ATTR_EMISSION, ATTR_PROB, TransitionGraph
+
+    maps = qr.symbol_maps
+    tau = np.asarray(qr.tau, dtype=float)
+    symbols = tuple(sorted(maps, key=repr))
+    eta0 = MixedState.from_vector(qr.pi, decimals=decimals)
+    if eta0 is None:
+        raise SpectralInferenceError("degenerate initial vector; cannot extract mixed states")
+
+    graph = TransitionGraph()
+    discovered: dict[MixedState, MixedState] = {}
+    queue: deque[MixedState] = deque()
+
+    def register(state: MixedState) -> MixedState:
+        existing = discovered.get(state)
+        if existing is not None:
+            return existing
+        atol = 10 ** (-decimals)
+        for known in discovered:
+            if all(np.isclose(a, b, rtol=0.0, atol=atol) for a, b in zip(known.belief, state.belief, strict=True)):
+                discovered[state] = known
+                return known
+        if len(discovered) >= max_states:
+            raise SpectralInferenceError(
+                f"mixed-state extraction exceeded max_states={max_states}; "
+                "use project_to_nmachine for the signed observable-operator model"
+            )
+        discovered[state] = state
+        graph.add_state(state)
+        queue.append(state)
+        return state
+
+    register(eta0)
+    while queue:
+        eta = queue.popleft()
+        row = eta.as_array()
+        emissions: list[tuple[Any, MixedState, float]] = []
+        for symbol in symbols:
+            nxt = row @ maps[symbol]
+            prob = float(nxt @ tau)
+            if prob <= _MASS_ATOL:
+                continue
+            successor = MixedState.from_vector(nxt, decimals=decimals)
+            if successor is None:
+                continue
+            emissions.append((symbol, register(successor), prob))
+        total = sum(prob for _symbol, _successor, prob in emissions)
+        if total <= _MASS_ATOL:
+            continue
+        for symbol, successor, prob in emissions:
+            graph.add_transition(
+                eta,
+                successor,
+                **{ATTR_PROB: float(prob / total), ATTR_EMISSION: symbol},
+            )
+
+    keep = graph.terminal_recurrent_states()
+    if not keep:
+        keep = frozenset(discovered.values())
+    recurrent = TransitionGraph()
+    for state in keep:
+        recurrent.add_state(state)
+        for transition in graph.out_transitions(state):
+            if transition.target in keep:
+                recurrent.add_transition(transition.source, transition.target, **dict(transition.data))
+    initial = {eta0: 1.0} if eta0 in keep else {}
+    return MealyHMM(
+        graph=recurrent,
+        initial_distribution=initial,
+        observation_alphabet=frozenset(symbols),
+    )
diff --git a/tests/test_epsilon_inference.py b/tests/test_epsilon_inference.py
index 00c4ca3..b025d24 100644
--- a/tests/test_epsilon_inference.py
+++ b/tests/test_epsilon_inference.py
@@ -10,7 +10,7 @@
 import pytest
 
 from sofic.examples.epsilon_machines import bernoulli, even_process, golden_mean
-from sofic.generators.epsilon_inference import cssr, subtree_merge
+from sofic.generators.epsilon_inference import cssr, spectral, subtree_merge
 from sofic.generators.epsilon_machine import EpsilonMachine
 from sofic.generators.hmm_inference import sample
 from sofic.graph import ATTR_EMISSION, ATTR_PROB
@@ -160,3 +160,47 @@ def test_cssr_short_sequence_raises():
 def test_subtree_merge_short_sequence_raises():
     with pytest.raises(ValueError, match="at least two"):
         subtree_merge([1], L=1)
+
+
+def test_spectral_bernoulli_single_state():
+    oracle = bernoulli(0.3)
+    alphabet = sorted(oracle.observation_alphabet, key=repr)
+    inferred = spectral(word_probability=oracle.word_probability, alphabet=alphabet, prefix_length=2, rank=1)
+    inferred.validate()
+    assert len(list(inferred.states())) == 1
+    assert inferred.entropy_rate() == pytest.approx(oracle.entropy_rate(), abs=1e-9)
+
+
+def test_spectral_golden_mean_recovers_two_states():
+    oracle = golden_mean(0.5)
+    alphabet = sorted(oracle.observation_alphabet, key=repr)
+    inferred = spectral(word_probability=oracle.word_probability, alphabet=alphabet, prefix_length=3, rank=2)
+    inferred.validate()
+    assert len(list(inferred.states())) == 2
+    assert inferred.entropy_rate() == pytest.approx(oracle.entropy_rate(), abs=1e-6)
+    assert inferred.statistical_complexity() == pytest.approx(oracle.statistical_complexity(), abs=1e-6)
+    assert _signatures_isomorphic(inferred, oracle, prob_tol=0.05)
+
+
+def test_spectral_even_process_recovers_two_states():
+    oracle = even_process(0.5)
+    alphabet = sorted(oracle.observation_alphabet, key=repr)
+    inferred = spectral(word_probability=oracle.word_probability, alphabet=alphabet, prefix_length=3, rank=2)
+    inferred.validate()
+    assert len(list(inferred.states())) == 2
+    assert inferred.entropy_rate() == pytest.approx(oracle.entropy_rate(), abs=1e-6)
+    assert inferred.statistical_complexity() == pytest.approx(oracle.statistical_complexity(), abs=1e-6)
+    assert _signatures_isomorphic(inferred, oracle, prob_tol=0.05)
+
+
+def test_from_sequence_spectral_dispatch(rng: np.random.Generator):
+    oracle = bernoulli(0.4)
+    observations, _ = sample(oracle, 400, rng)
+    inferred = EpsilonMachine.from_sequence(observations, method="spectral", prefix_length=2, rank=1)
+    inferred.validate()
+    assert len(list(inferred.states())) == 1
+
+
+def test_from_sequence_unknown_method():
+    with pytest.raises(ValueError, match="unknown inference method"):
+        EpsilonMachine.from_sequence([0, 1, 0], method="nsd")
diff --git a/tests/test_spectral_inference.py b/tests/test_spectral_inference.py
index 98a2e65..4350830 100644
--- a/tests/test_spectral_inference.py
+++ b/tests/test_spectral_inference.py
@@ -14,6 +14,7 @@
 from sofic.inference.spectral import (
     SpectralInferenceError,
     learn_spectral_wfa,
+    project_to_epsilon_machine,
     project_to_mealy,
     project_to_nmachine,
     spectral_singular_values,
@@ -134,6 +135,27 @@ def test_project_to_mealy_rejects_signed_realization():
         project_to_mealy(learned)
 
 
+def test_project_to_epsilon_machine_from_signed_golden_mean():
+    from sofic.generators.epsilon_machine import EpsilonMachine
+
+    model = golden_mean(0.5)
+    alphabet = sorted(model.observation_alphabet, key=repr)
+    learned = learn_spectral_wfa(word_probability=model.word_probability, alphabet=alphabet, prefix_length=3, rank=2)
+    machine = project_to_epsilon_machine(learned)
+    assert isinstance(machine, EpsilonMachine)
+    machine.validate()
+    assert len(list(machine.states())) == 2
+    assert machine.entropy_rate() == pytest.approx(model.entropy_rate(), abs=1e-6)
+
+
+def test_project_to_epsilon_machine_respects_max_states():
+    model = golden_mean(0.5)
+    alphabet = sorted(model.observation_alphabet, key=repr)
+    learned = learn_spectral_wfa(word_probability=model.word_probability, alphabet=alphabet, prefix_length=3, rank=2)
+    with pytest.raises(SpectralInferenceError, match="max_states"):
+        project_to_epsilon_machine(learned, max_states=1)
+
+
 # --- Input validation ------------------------------------------------------