diff --git a/docs/api/features/chimeric.md b/docs/api/features/chimeric.md
index 6174b2b..860729d 100644
--- a/docs/api/features/chimeric.md
+++ b/docs/api/features/chimeric.md
@@ -48,6 +48,7 @@ All columns from `FragmentMatchFeatures` with `chimeric_` prefix:
 | `chimeric_max_ion_gap` | Daltons (Da) | Largest gap between matched runner-up ions |
 | `chimeric_b_y_intensity_ratio` | Ratio | Ratio of b-ion to y-ion intensity for runner-up (including isotopic envelopes) |
 | `chimeric_spectral_angle` | Score (0-1) | Normalised spectral angle similarity between runner-up theoretical and observed intensities |
+| `chimeric_xcorr` | Score | SEQUEST fast cross-correlation score for the runner-up peptide. Measures overall agreement between the observed spectrum and the runner-up theoretical spectrum with local background correction. See [FragmentMatchFeatures — Cross-correlation Score](fragment_match.md#cross-correlation-score) for details on the algorithm. |
 
 ```python
 from winnow.calibration.features import ChimericFeatures
diff --git a/docs/api/features/fragment_match.md b/docs/api/features/fragment_match.md
index a15d3c0..a09712a 100644
--- a/docs/api/features/fragment_match.md
+++ b/docs/api/features/fragment_match.md
@@ -57,6 +57,21 @@ For each theoretical peak, we search for the nearest observed peak using binary
 | `max_ion_gap` | Daltons (Da) | Largest m/z difference between two consecutive matched theoretical peaks when sorted by m/z. Large gaps may indicate missing fragmentation coverage. |
 | `b_y_intensity_ratio` | Ratio | Ratio of total matched b-ion intensity to total matched y-ion intensity (including isotopic envelopes). |
 | `spectral_angle` | Score (0-1) | Normalised spectral angle similarity between theoretical and matched observed intensity vectors. A value of 1 indicates perfect correlation, 0 indicates orthogonal vectors. |
+| `xcorr` | Score | SEQUEST fast cross-correlation score. Measures overall agreement between the observed and theoretical spectra with local background correction. Higher values indicate better matches. |
+
+### Cross-correlation Score
+
+The `xcorr` column implements the fast cross-correlation score function from SEQUEST (Eng et al., 2008). Unlike the other features which operate on individually matched peaks, xcorr evaluates the overall pattern of the observed spectrum against the theoretical spectrum using a background-corrected dot product.
+
+The observed spectrum is preprocessed by:
+
+1. **Binning** into near-unit-dalton bins (~1.0005 Da) with square-root intensity compression
+2. **Window normalization** — the maximum intensity within each of 10 equal m/z windows is normalized to a fixed value, making the score robust to different intensity scales
+3. **Background subtraction** — the mean intensity from ±75 neighbouring bins is subtracted from each bin, so that only signal specifically at the correct m/z positions contributes positively
+
+The score is then computed as the dot product of the preprocessed observed spectrum with a binary theoretical spectrum (1.0 at each predicted fragment ion bin position). This background correction naturally penalizes unmatched theoretical ions and discounts noise.
+
+**Reference:** Eng JK, Fischer B, Grossmann J, Maccoss MJ. A fast SEQUEST cross correlation algorithm. *J Proteome Res.* 2008 Oct;7(10):4598-602. doi: 10.1021/pr800420s.
 
 ## Usage
 
diff --git a/tests/calibration/features/test_chimeric.py b/tests/calibration/features/test_chimeric.py
index 548a148..50d9d2e 100644
--- a/tests/calibration/features/test_chimeric.py
+++ b/tests/calibration/features/test_chimeric.py
@@ -80,6 +80,7 @@ def test_properties(self, chimeric_features):
             "chimeric_max_ion_gap",
             "chimeric_b_y_intensity_ratio",
             "chimeric_spectral_angle",
+            "chimeric_xcorr",
             # Missing indicator (learn_from_missing=True by default)
             "is_missing_chimeric_features",
         ]
@@ -158,6 +159,7 @@ def test_learn_from_missing_false_columns_excludes_indicator(self):
             "chimeric_max_ion_gap",
             "chimeric_b_y_intensity_ratio",
             "chimeric_spectral_angle",
+            "chimeric_xcorr",
         ]
 
     # ------------------------------------------------------------------
@@ -232,6 +234,7 @@ def test_compute_preprocessing_pipeline(
         mock_max_ion_gap = 1
         mock_b_y_intensity_ratio = 0.5
         mock_spectral_angle = 0.8
+        mock_xcorr = [1.0, 2.0]
         mock_compute_ions.return_value = IonIdentificationResult(
             mock_match_rate,
             mock_match_intensity,
@@ -241,6 +244,7 @@ def test_compute_preprocessing_pipeline(
             mock_max_ion_gap,
             mock_b_y_intensity_ratio,
             mock_spectral_angle,
+            mock_xcorr,
         )
 
         # Run the compute method
diff --git a/tests/calibration/features/test_fragment_match.py b/tests/calibration/features/test_fragment_match.py
index 786f1c8..066e974 100644
--- a/tests/calibration/features/test_fragment_match.py
+++ b/tests/calibration/features/test_fragment_match.py
@@ -66,6 +66,7 @@ def test_properties(self, prosit_features):
             "max_ion_gap",
             "b_y_intensity_ratio",
             "spectral_angle",
+            "xcorr",
             # Missing indicator (learn_from_missing=True by default)
             "is_missing_fragment_match_features",
         ]
@@ -143,6 +144,7 @@ def test_learn_from_missing_false_columns_excludes_indicator(self):
             "max_ion_gap",
             "b_y_intensity_ratio",
             "spectral_angle",
+            "xcorr",
         ]
 
     def test_learn_from_missing_true_columns_includes_indicator(self):
@@ -256,6 +258,7 @@ def test_compute_with_mock(
         mock_max_ion_gap = 1
         mock_b_y_intensity_ratio = 0.5
         mock_spectral_angle = 0.8
+        mock_xcorr = [1.0, 2.0, 3.0]
         mock_compute_ions.return_value = IonIdentificationResult(
             mock_match_rate,
             mock_match_intensity,
@@ -265,6 +268,7 @@ def test_compute_with_mock(
             mock_max_ion_gap,
             mock_b_y_intensity_ratio,
             mock_spectral_angle,
+            mock_xcorr,
         )
 
         prosit_features.compute(sample_dataset_with_spectra)
diff --git a/tests/calibration/features/test_utils.py b/tests/calibration/features/test_utils.py
index 56e89ac..221def8 100644
--- a/tests/calibration/features/test_utils.py
+++ b/tests/calibration/features/test_utils.py
@@ -17,6 +17,7 @@
     compute_max_ion_gap,
     compute_b_y_intensity_ratio,
     compute_spectral_angle,
+    compute_xcorr,
     parse_mz_tolerance_unit,
     _validate_mz_tolerance,
 )
@@ -969,3 +970,115 @@ def test_compute_spectral_angle_length_mismatch_raises(self):
         observed = [1.0, 2.0]  # Different length
         with pytest.raises(ValueError, match="must align"):
             compute_spectral_angle(theoretical, observed)
+
+
+class TestXcorr:
+    """Test the fast xcorr score function."""
+
+    def test_matching_peaks_positive_score(self):
+        """Peaks exactly at theoretical positions should produce a positive xcorr."""
+        theoretical_mz = [200.0, 400.0, 600.0]
+        observed_mz = [200.0, 400.0, 600.0]
+        observed_intensities = [1000.0, 2000.0, 1500.0]
+
+        score = compute_xcorr(observed_mz, observed_intensities, theoretical_mz)
+        assert score > 0.0
+
+    def test_no_matching_peaks_low_score(self):
+        """Peaks far from theoretical positions should yield a low or negative xcorr."""
+        theoretical_mz = [200.0, 400.0, 600.0]
+        observed_mz = [250.0, 450.0, 650.0]
+        observed_intensities = [1000.0, 2000.0, 1500.0]
+
+        matching_score = compute_xcorr(
+            [200.0, 400.0, 600.0], observed_intensities, [200.0, 400.0, 600.0]
+        )
+        non_matching_score = compute_xcorr(
+            observed_mz, observed_intensities, theoretical_mz
+        )
+        assert non_matching_score < matching_score
+
+    def test_nan_theoretical_returns_zero(self):
+        """NaN theoretical m/z (missing prediction) returns 0.0."""
+        import math
+
+        score = compute_xcorr([100.0], [1000.0], math.nan)
+        assert score == 0.0
+
+    def test_empty_theoretical_returns_zero(self):
+        """Empty theoretical m/z list returns 0.0."""
+        score = compute_xcorr([100.0], [1000.0], [])
+        assert score == 0.0
+
+    def test_empty_observed_returns_zero(self):
+        """Empty observed spectrum returns 0.0."""
+        score = compute_xcorr([], [], [100.0, 200.0])
+        assert score == 0.0
+
+    def test_uniform_spectrum_much_lower_than_matching(self):
+        """A uniform observed spectrum should score much lower than a matching one.
+
+        When every bin has the same intensity, the background subtraction
+        largely cancels the signal. Some residual remains due to edge effects
+        in the background window near spectrum boundaries, but the score
+        should be much lower than when peaks align with theoretical positions.
+        """
+        observed_mz = [float(i) for i in range(50, 500)]
+        observed_intensities = [100.0] * len(observed_mz)
+        theoretical_mz = [100.0, 200.0, 300.0]
+
+        uniform_score = compute_xcorr(observed_mz, observed_intensities, theoretical_mz)
+
+        sparse_mz = [100.0, 200.0, 300.0]
+        sparse_intensities = [1000.0, 1000.0, 1000.0]
+        matching_score = compute_xcorr(sparse_mz, sparse_intensities, theoretical_mz)
+
+        assert uniform_score < matching_score * 0.25
+
+    def test_scale_invariance(self):
+        """Xcorr should be the same regardless of observed intensity scale.
+
+        Window normalization rescales each region to a fixed value, so
+        multiplying all intensities by a constant should not change the score.
+        """
+        theoretical_mz = [200.0, 400.0, 600.0]
+        observed_mz = [200.0, 300.0, 400.0, 500.0, 600.0]
+        intensities_1x = [1000.0, 500.0, 2000.0, 300.0, 1500.0]
+        intensities_100x = [i * 100.0 for i in intensities_1x]
+
+        score_1x = compute_xcorr(observed_mz, intensities_1x, theoretical_mz)
+        score_100x = compute_xcorr(observed_mz, intensities_100x, theoretical_mz)
+        assert score_1x == pytest.approx(score_100x, rel=1e-10)
+
+    def test_more_matches_higher_score(self):
+        """More matching peaks should produce a higher xcorr than fewer matches."""
+        observed_mz = [100.0, 200.0, 300.0, 400.0, 500.0]
+        observed_intensities = [1000.0] * 5
+
+        score_2_matches = compute_xcorr(
+            observed_mz, observed_intensities, [100.0, 200.0]
+        )
+        score_4_matches = compute_xcorr(
+            observed_mz, observed_intensities, [100.0, 200.0, 300.0, 400.0]
+        )
+        assert score_4_matches > score_2_matches
+
+    def test_observed_preprocessing_independent_of_out_of_range_theoretical_ions(self):
+        """Unmatched high-m/z theoretical ions must not change observed normalization."""
+        observed_mz = [166.51, 239.26, 280.07, 391.08, 429.99]
+        observed_intensities = [416.6, 13651.9, 163.4, 806.1, 2155.4]
+        theoretical_mz = [239.26, 429.99, 391.08]
+
+        score = compute_xcorr(observed_mz, observed_intensities, theoretical_mz)
+        score_with_unmatched_ion = compute_xcorr(
+            observed_mz,
+            observed_intensities,
+            theoretical_mz + [2000.0],
+        )
+
+        assert score_with_unmatched_ion == pytest.approx(score)
+
+    def test_single_peak_positive(self):
+        """A single observed peak matching a single theoretical ion should score positive."""
+        score = compute_xcorr([300.0], [5000.0], [300.0])
+        assert score > 0.0
diff --git a/winnow/calibration/features/chimeric.py b/winnow/calibration/features/chimeric.py
index 8256e6f..d364bf0 100644
--- a/winnow/calibration/features/chimeric.py
+++ b/winnow/calibration/features/chimeric.py
@@ -131,6 +131,7 @@ def columns(self) -> List[str]:
             "chimeric_max_ion_gap",
             "chimeric_b_y_intensity_ratio",
             "chimeric_spectral_angle",
+            "chimeric_xcorr",
         ]
         if self.learn_from_missing:
             columns.append("is_missing_chimeric_features")
@@ -338,3 +339,4 @@ def compute(self, dataset: CalibrationDataset) -> None:
             ion_identifications.b_y_intensity_ratio
         )
         dataset.metadata["chimeric_spectral_angle"] = ion_identifications.spectral_angle
+        dataset.metadata["chimeric_xcorr"] = ion_identifications.xcorr
diff --git a/winnow/calibration/features/constants.py b/winnow/calibration/features/constants.py
index cb1a10d..4935401 100644
--- a/winnow/calibration/features/constants.py
+++ b/winnow/calibration/features/constants.py
@@ -7,3 +7,13 @@
 PROTON_MASS = 1.007276
 
 VALID_INTENSITY_MODEL_PROVIDERS = ["prosit", "ms2pip", "alphapeptdeep"]
+# Default bin width matching Comet's fragment_bin_tol for near-unit-dalton bins.
+XCORR_BIN_SIZE = 1.0005079
+# Default bin offset matching Comet's fragment_bin_offset.
+XCORR_BIN_OFFSET = 0.4
+# Number of equal-width m/z windows for intensity normalization (SEQUEST default).
+XCORR_NUM_WINDOWS = 10
+# Maximum m/z offset in bins for background subtraction (SEQUEST default ±75).
+XCORR_MAX_OFFSET = 75
+# Normalization target for the maximum intensity within each window.
+XCORR_WINDOW_NORM_VALUE = 50.0
diff --git a/winnow/calibration/features/fragment_match.py b/winnow/calibration/features/fragment_match.py
index a694717..4d68168 100644
--- a/winnow/calibration/features/fragment_match.py
+++ b/winnow/calibration/features/fragment_match.py
@@ -130,6 +130,7 @@ def columns(self) -> List[str]:
             "max_ion_gap",
             "b_y_intensity_ratio",
             "spectral_angle",
+            "xcorr",
         ]
         if self.learn_from_missing:
             columns.append("is_missing_fragment_match_features")
@@ -313,3 +314,4 @@ def compute(self, dataset: CalibrationDataset) -> None:
             ion_identifications.b_y_intensity_ratio
         )
         dataset.metadata["spectral_angle"] = ion_identifications.spectral_angle
+        dataset.metadata["xcorr"] = ion_identifications.xcorr
diff --git a/winnow/calibration/features/utils.py b/winnow/calibration/features/utils.py
index 5d5f80f..042bd4a 100644
--- a/winnow/calibration/features/utils.py
+++ b/winnow/calibration/features/utils.py
@@ -26,6 +26,11 @@
 from winnow.calibration.features.constants import (
     CARBON_ISOTOPE_MASS_SHIFT,
     VALID_INTENSITY_MODEL_PROVIDERS,
+    XCORR_BIN_SIZE,
+    XCORR_BIN_OFFSET,
+    XCORR_NUM_WINDOWS,
+    XCORR_MAX_OFFSET,
+    XCORR_WINDOW_NORM_VALUE,
 )
 
 ########################################################
@@ -64,6 +69,7 @@ class IonIdentificationResult(NamedTuple):
     max_ion_gap: Tuple[float, ...]
     b_y_intensity_ratio: Tuple[float, ...]
     spectral_angle: Tuple[float, ...]
+    xcorr: Tuple[float, ...]
 
     def __iter__(self):
         raise TypeError(
@@ -72,7 +78,8 @@ def __iter__(self):
             "by accessing result fields by name, e.g. result.ion_match_rate, "
             "result.ion_match_intensity, result.longest_b_series, "
             "result.longest_y_series, result.complementary_ion_count, "
-            "result.max_ion_gap, result.b_y_intensity_ratio, and result.spectral_angle."
+            "result.max_ion_gap, result.b_y_intensity_ratio, "
+            "result.spectral_angle, and result.xcorr."
         )
 
 
@@ -500,7 +507,7 @@ def compute_ion_identifications(
     mz_tolerance_unit: str,
     predictions: Optional[List[str]] = None,
 ) -> IonIdentificationResult:
-    """Computes the ion match rate, match intensity, longest b series, longest y series, complementary ion count, max ion gap and b/y intensity ratio and normalised spectral angle for each spectrum in the dataset.
+    """Computes the ion match rate, match intensity, longest b series, longest y series, complementary ion count, max ion gap, b/y intensity ratio, normalised spectral angle and fast cross-correlation score for each spectrum in the dataset.
 
     Args:
         dataset: DataFrame containing the mass spectrum data.
@@ -521,7 +528,7 @@ def compute_ion_identifications(
     _validate_mz_tolerance(mz_tolerance, mz_tolerance_unit)
 
     per_row_match_results: List[
-        Tuple[float, float, int, int, int, float, float, float]
+        Tuple[float, float, int, int, int, float, float, float, float]
     ] = []
 
     for row_idx, (_, row) in enumerate(dataset.iterrows()):
@@ -562,6 +569,12 @@ def compute_ion_identifications(
         spectral_angle = compute_spectral_angle(
             row[source_intensity_column], match.aligned_m0_intensities
         )
+        # Compute the fast xcorr score
+        xcorr = compute_xcorr(
+            observed_mz=row["mz_array"],
+            observed_intensities=row["intensity_array"],
+            theoretical_mz=row[source_mz_column],
+        )
 
         per_row_match_results.append(
             (
@@ -573,11 +586,12 @@ def compute_ion_identifications(
                 max_ion_gap,
                 b_y_intensity_ratio,
                 spectral_angle,
+                xcorr,
             )
         )
 
     if not per_row_match_results:
-        return IonIdentificationResult((), (), (), (), (), (), (), ())
+        return IonIdentificationResult((), (), (), (), (), (), (), (), ())
 
     return IonIdentificationResult(*zip(*per_row_match_results))
 
@@ -790,3 +804,152 @@ def compute_spectral_angle(
     dot_product = np.clip(dot_product, -1.0, 1.0)
 
     return 1 - (2 * np.arccos(dot_product) / np.pi)
+
+
+########################################################
+# Cross-correlation Score (xcorr)
+########################################################
+
+
+def _bin_observed_spectrum(
+    observed_mz: List[float],
+    observed_intensities: List[float],
+    num_bins: int,
+    bin_size: float,
+    bin_offset: float,
+) -> np.ndarray:
+    """Bin an observed spectrum with sqrt intensity compression.
+
+    Places each peak into a near-unit-dalton bin, taking the square root of
+    its intensity. When multiple peaks land in the same bin, the maximum
+    sqrt-intensity is kept.
+
+    Args:
+        observed_mz: m/z values of the observed spectrum.
+        observed_intensities: Intensities corresponding to observed_mz.
+        num_bins: Total number of bins in the output array.
+        bin_size: Width of each m/z bin in Daltons.
+        bin_offset: Fractional offset for bin assignment.
+
+    Returns:
+        Array of length num_bins with sqrt-compressed, max-per-bin intensities.
+    """
+    binned = np.zeros(num_bins)
+    for mz, intensity in zip(observed_mz, observed_intensities):
+        bin_idx = int(mz / bin_size + bin_offset)
+        if 0 <= bin_idx < num_bins:
+            sqrt_intensity = np.sqrt(max(intensity, 0.0))
+            if sqrt_intensity > binned[bin_idx]:
+                binned[bin_idx] = sqrt_intensity
+    return binned
+
+
+def _normalize_windows(
+    spectrum: np.ndarray,
+    num_windows: int,
+) -> None:
+    """Normalize maximum intensity within equal m/z windows (in-place).
+
+    Divides the spectrum into num_windows equal regions and scales each so
+    that its maximum intensity equals XCORR_WINDOW_NORM_VALUE.
+
+    Args:
+        spectrum: Binned spectrum array, modified in-place.
+        num_windows: Number of equal-width windows.
+    """
+    num_bins = len(spectrum)
+    window_size = max(num_bins // num_windows, 1)
+    for w in range(num_windows):
+        start = w * window_size
+        end = min(start + window_size, num_bins) if w < num_windows - 1 else num_bins
+        window = spectrum[start:end]
+        max_val = window.max()
+        if max_val > 0:
+            spectrum[start:end] = window * (XCORR_WINDOW_NORM_VALUE / max_val)
+
+
+def _subtract_background(
+    spectrum: np.ndarray,
+    max_xcorr_offset: int,
+) -> np.ndarray:
+    """Subtract local background using the fast xcorr preprocessing.
+
+    For each bin i, computes:
+        y'[i] = y[i] - (sum of y[i-offset..i+offset] excluding y[i]) / (2 * offset)
+
+    Uses a cumulative sum for O(n) efficiency.
+
+    Args:
+        spectrum: Window-normalized binned spectrum.
+        max_xcorr_offset: Number of bins in each direction for background.
+
+    Returns:
+        Background-subtracted spectrum.
+    """
+    num_bins = len(spectrum)
+    cumsum = np.concatenate(([0.0], np.cumsum(spectrum)))
+    lo = np.maximum(0, np.arange(num_bins) - max_xcorr_offset)
+    hi = np.minimum(num_bins, np.arange(num_bins) + max_xcorr_offset + 1)
+    window_sums = cumsum[hi] - cumsum[lo]
+    divisor = 2 * max_xcorr_offset
+    return spectrum - (window_sums - spectrum) / divisor
+
+
+def compute_xcorr(
+    observed_mz: List[float],
+    observed_intensities: List[float],
+    theoretical_mz: Union[List[float], float],
+    bin_size: float = XCORR_BIN_SIZE,
+    bin_offset: float = XCORR_BIN_OFFSET,
+    num_windows: int = XCORR_NUM_WINDOWS,
+    max_xcorr_offset: int = XCORR_MAX_OFFSET,
+) -> float:
+    """Compute the SEQUEST fast xcorr score between observed and theoretical spectra.
+
+    Implements the fast cross-correlation score from Eng et al. (2008). The observed
+    spectrum is preprocessed by binning into near-unit-dalton bins with sqrt intensity
+    compression, normalizing peak intensities within equal m/z windows, and subtracting
+    the local background (mean of ±max_xcorr_offset shifted spectra). The xcorr score
+    is the dot product of the preprocessed observed spectrum with a binary theoretical
+    spectrum that has 1.0 at each predicted fragment ion bin.
+
+    Reference:
+        Eng JK, Fischer B, Grossmann J, Maccoss MJ. A fast SEQUEST cross correlation
+        algorithm. J Proteome Res. 2008 Oct;7(10):4598-602.
+
+    Args:
+        observed_mz: m/z values of the observed (acquired) spectrum.
+        observed_intensities: Intensities corresponding to observed_mz.
+        theoretical_mz: m/z values of theoretical fragment ions from intensity
+            prediction. May be NaN (float) for missing predictions.
+        bin_size: Width of each m/z bin in Daltons. Defaults to 1.0005079.
+        bin_offset: Fractional offset for bin assignment. Defaults to 0.4.
+        num_windows: Number of equal-width m/z windows for intensity normalization.
+            Defaults to 10.
+        max_xcorr_offset: Number of bins in each direction for background estimation.
+            The background is divided by 2 * max_xcorr_offset. Defaults to 75.
+
+    Returns:
+        The xcorr score. Returns 0.0 for missing/empty inputs.
+    """
+    if isinstance(theoretical_mz, float):
+        return 0.0
+    if len(observed_mz) == 0 or len(theoretical_mz) == 0:
+        return 0.0
+
+    max_mz = max(observed_mz)
+    num_bins = int(max_mz / bin_size + bin_offset) + 1
+
+    observed_binned = _bin_observed_spectrum(
+        observed_mz, observed_intensities, num_bins, bin_size, bin_offset
+    )
+    _normalize_windows(observed_binned, num_windows)
+    processed = _subtract_background(observed_binned, max_xcorr_offset)
+
+    theoretical_binned = np.zeros(num_bins)
+    for mz in theoretical_mz:
+        bin_idx = int(mz / bin_size + bin_offset)
+        if 0 <= bin_idx < num_bins:
+            theoretical_binned[bin_idx] = 1.0
+
+    return float(np.dot(theoretical_binned, processed))