Add sdr.preferred_pairs()

Fixes #384

src/sdr/_sequence/_maximum.py

@@ -4,7 +4,8 @@
 
 from __future__ import annotations
 
-from typing import Any, overload
+import itertools
+from typing import Any, Iterator, overload
 
 import galois
 import numpy as np
@@ -164,3 +165,128 @@ def m_sequence(
  if not q == 2:
  raise ValueError(f"Argument 'output' can only be 'bipolar' when q = 2, not {q}.")
  return code_to_sequence(field_to_code(code))
+
+
+@export
+def preferred_pairs(
+ degree: int,
+ poly1: PolyLike | None = None,
+) -> Iterator[tuple[Poly, Poly]]:
+ r"""
+ Generates primitive polynomials of degree $m$ that produce preferred pair $m$-sequences.
+
+ Arguments:
+ degree: The degree $m$ of the $m$-sequences.
+ poly1: The first polynomial $f(x)$ in the preferred pair. If `None`, all primitive polynomials of degree $m$
+ that yield preferred pair $m$-sequences are returned.
+
+ Returns:
+ An iterator of primitive polynomials $(f(x), g(x))$ of degree $m$ that produce preferred pair $m$-sequences.
+
+ See Also:
+ sdr.gold_code
+
+ Notes:
+ A preferred pair of primitive polynomials of degree $m$ are two polynomials $f(x)$ and $g(x)$ such that the
+ periodic cross-correlation of the two $m$-sequences generated by $f(x)$ and $g(x)$ have only the values in
+ $\{-1, -t(m), t(m) - 2\}$.
+
+ $$
+ t(m) = \begin{cases}
+ 2^{(m+1)/2} + 1 & \text{if $m$ is odd} \\
+ 2^{(m+2)/2} + 1 & \text{if $m$ is even}
+ \end{cases}
+ $$
+
+ There are no preferred pairs for degrees divisible by 4.
+
+ References:
+ - John Proakis, *Digital Communications*, Chapter 12.2-5: Generation of PN Sequences.
+
+ Examples:
+ Generate one preferred pair with $f(x) = x^5 + x^3 + 1$.
+
+ .. ipython:: python
+
+ next(sdr.preferred_pairs(5, poly1="x^5 + x^3 + 1"))
+
+ Generate all preferred pairs with $f(x) = x^5 + x^3 + 1$.
+
+ .. ipython:: python
+
+ list(sdr.preferred_pairs(5, poly1="x^5 + x^3 + 1"))
+
+ Generate all preferred pairs with degree 5.
+
+ .. ipython:: python
+
+ list(sdr.preferred_pairs(5))
+
+ Note, there are no preferred pairs for degrees divisible by 4.
+
+ .. ipython:: python
+
+ list(sdr.preferred_pairs(4))
+ list(sdr.preferred_pairs(8))
+
+ Group:
+ sequences-maximum-length
+ """
+ if not isinstance(degree, int):
+ raise TypeError(f"Argument 'degree' must be an integer, not {type(degree)}.")
+ if not degree > 0:
+ raise ValueError(f"Argument 'degree' must be positive, not {degree}.")
+
+ if degree % 4 == 0:
+ # There are no preferred pairs for degrees divisible by 4
+ return
+
+ # Compute t(m) for degree m, Equation 12.2-73
+ if degree % 2 == 1:
+ t_m = 2 ** ((degree + 1) // 2) + 1
+ else:
+ t_m = 2 ** ((degree + 2) // 2) + 1
+
+ # Determine the valid cross-correlation values for preferred pairs, Page 799
+ valid_values = [-1, -t_m, t_m - 2]
+
+ if poly1 is None:
+ # Find all combinations of primitive polynomials of degree m
+ for poly1, poly2 in itertools.combinations(galois.primitive_polys(2, degree), 2):
+ # Create first m-sequence with the first polynomial
+ u = m_sequence(degree, poly1, output="bipolar")
+ U = np.fft.fft(u)
+
+ # Create second m-sequence with the second polynomial
+ v = m_sequence(degree, poly2, output="bipolar")
+ V = np.fft.fft(v)
+
+ # Compute the periodic cross-correlation of the two m-sequences
+ pccf = np.fft.ifft(U * V.conj()).real # The inputs are real, so the output is real
+ pccf = np.around(pccf).astype(int) # Round to the nearest integer so isin() works
+
+ if np.all(np.isin(pccf, valid_values)):
+ yield poly1, poly2
+ else:
+ # Find all combinations of the first polynomial with all primitive polynomials of degree m
+ poly1 = Poly._PolyLike(poly1, field=galois.GF(2))
+ if not poly1.degree == degree:
+ raise ValueError(f"Argument 'poly1' must be a polynomial of degree {degree}, not {poly1.degree}.")
+ if not poly1.is_primitive():
+ raise ValueError(f"Argument 'poly1' must be a primitive polynomial, {poly1} is not.")
+
+ # Create first m-sequence with the first polynomial
+ u = m_sequence(degree, poly1, output="bipolar")
+ U = np.fft.fft(u)
+
+ for poly2 in galois.primitive_polys(2, degree):
+ # Create second m-sequence with the second polynomial
+ v = m_sequence(degree, poly2, output="bipolar")
+ V = np.fft.fft(v)
+
+ # Compute the periodic cross-correlation of the two m-sequences
+ pccf = np.fft.ifft(U * V.conj()).real # The inputs are real, so the output is real
+ pccf = np.around(pccf).astype(int) # Round to the nearest integer so isin() works
+
+ if np.all(np.isin(pccf, valid_values)):
+ yield poly1, poly2