Release 0.0.23

README.md

@@ -69,7 +69,8 @@ View all available classes and functions in the [API Reference](https://mhostett
   Apply additive white Gaussian noise (AWGN), frequency offset, sample rate offset, IQ imbalance.
 - **Link budgets**: BSC, BEC, AWGN, and BI-AWGN channel capacity, Shannon's limit on $E_b/N_0$ and $S/N$,
   free-space path loss, antenna gain.
-- **Miscellaneous**: Packing and unpacking binary data, hexdump of binary data.
+- **Miscellaneous**: Numerical calculation of PDF of sum and product of random variables. Packing and unpacking binary
+  data, hexdump of binary data.
 - **Plotting**: Time-domain, raster, correlation, stem, discrete Fourier transform (DFT), discrete-time Fourier
   transform (DTFT), periodogram, spectrogram, constellation, symbol map, eye diagram, phase tree,
   bit error rate (BER), symbol error rate (SER), probability of detection, receiver operating characteristic (ROC),

docs/index.rst

@@ -79,7 +79,8 @@ View all available classes and functions in the `API Reference <https://mhostett
   Apply additive white Gaussian noise (AWGN), frequency offset, sample rate offset, IQ imbalance.
 - **Link budgets**: BSC, BEC, AWGN, and BI-AWGN channel capacity, Shannon's limit on $E_b/N_0$ and $S/N$,
   free-space path loss, antenna gain.
-- **Miscellaneous**: Packing and unpacking binary data, hexdump of binary data.
+- **Miscellaneous**: Numerical calculation of PDF of sum and product of random variables. Packing and unpacking binary
+  data, hexdump of binary data.
 - **Plotting**: Time-domain, raster, correlation, stem, discrete Fourier transform (DFT), discrete-time Fourier
   transform (DTFT), periodogram, spectrogram, constellation, symbol map, eye diagram, phase tree,
   bit error rate (BER), symbol error rate (SER), probability of detection, receiver operating characteristic (ROC),

docs/release-notes/v0.0.md

@@ -4,6 +4,29 @@ tocdepth: 2
 
 # v0.0
 
+## v0.0.23
+
+*Released July 28, 2024*
+
+### Changes
+
+- Added calculation of filter noise bandwidth in `sdr.FIR.noise_bandwidth()` and `sdr.IIR.noise_bandwidth()`.
+- Added calculation of threshold level above the noise power in `sdr.threshold_factor()`.
+- Added numerical calculation of the PDF of the sum and product of random variables in `sdr.sum_distribution()`,
+  `sdr.sum_distributions()`, and `sdr.multiply_distributions()`.
+- Renamed `sdr.design_frac_delay_fir()` to `sdr.fractional_delay_fir()`.
+- Renamed `sdr.design_lowpass_fir()` to `sdr.lowpass_fir()`.
+- Renamed `sdr.design_highpass_fir()` to `sdr.highpass_fir()`.
+- Renamed `sdr.design_bandpass_fir()` to `sdr.bandpass_fir()`.
+- Renamed `sdr.design_bandstop_fir()` to `sdr.bandstop_fir()`.
+- Renamed `sdr.design_multirate_fir()` to `sdr.multirate_fir()`.
+- Allowed use of SciPy window definition using `scipy.signal.windows.get_window()` for all filter design and plotting
+  functions.
+
+### Contributors
+
+- Matt Hostetter ([@mhostetter](https://github.com/mhostetter))
+
 ## v0.0.22
 
 *Released July 13, 2024*

docs/requirements.txt

@@ -1,5 +1,5 @@
 sphinx == 7.3.7
-sphinx-immaterial == 0.11.14
+sphinx-immaterial == 0.12.1
 # For some reason, if the below versions aren't specified, the dependency resolution takes 18 minutes in CI.
 myst-parser == 3.0.1
 sphinx-design  == 0.6.0

src/sdr/_data.py

@@ -129,7 +129,7 @@ def hexdump(
 
     Arguments:
         data: The data to display. Each element is considered one byte. Use :func:`sdr.pack()`
-            or :func:`sdr.unpack()` to convert data with variable bits per element.
+            or :func:`sdr.unpack` to convert data with variable bits per element.
         width: The number of bytes per line.
 
     Returns:

src/sdr/_detection/_coherent_integration.py

@@ -18,7 +18,7 @@ def coherent_gain(time_bandwidth: npt.ArrayLike) -> npt.NDArray[np.float64]:
     Computes the SNR improvement by coherent integration.
 
     Arguments:
-        time_bandwidth: The time-bandwidth product $T_C B_n$ in seconds-Hz (unitless). If the noise bandwidth equals
+        time_bandwidth: The time-bandwidth product $T_c B_n$ in seconds-Hz (unitless). If the noise bandwidth equals
             the sample rate, the argument equals the number of samples $N_c$ to coherently integrate.
 
     Returns:
@@ -36,7 +36,7 @@ def coherent_gain(time_bandwidth: npt.ArrayLike) -> npt.NDArray[np.float64]:
 
         The coherent integration gain is the time-bandwidth product
 
-        $$G_c = 10 \log_{10} (T_C B_n) .$$
+        $$G_c = 10 \log_{10} (T_c B_n) .$$
 
         If the noise bandwidth equals the sample rate, the coherent gain is simply
 
@@ -86,8 +86,8 @@ def coherent_gain_loss(
     Computes the coherent gain loss (CGL) given a time or frequency offset.
 
     Arguments:
-        time: The coherent integration time $T_C$ or time offset in $\Delta t$ in seconds.
-        freq: The frequency offset $\Delta f$ or signal bandwidth $B_C$ in Hz.
+        time: The coherent integration time $T_c$ or time offset in $\Delta t$ in seconds.
+        freq: The frequency offset $\Delta f$ or signal bandwidth $B_c$ in Hz.
 
     Returns:
         The coherent gain loss (CGL) in dB.
@@ -96,18 +96,18 @@ def coherent_gain_loss(
         Coherent gain loss is the reduction in SNR due to time or frequency offsets during coherent integration.
         These losses are similar to scalloping loss.
 
-        The coherent gain loss of a signal integrated for $T_C$ seconds with a frequency offset of $\Delta f$ Hz is
+        The coherent gain loss of a signal integrated for $T_c$ seconds with a frequency offset of $\Delta f$ Hz is
 
         $$\text{CGL} = -10 \log_{10} \left( \text{sinc}^2 \left( T_c \Delta f \right) \right) ,$$
 
         where the sinc function is defined as
 
         $$\text{sinc}(x) = \frac{\sin(\pi x)}{\pi x} .$$
 
-        The coherent gain loss of a signal with bandwidth $B_C$ Hz with a detection time offset of $\Delta t$ seconds
+        The coherent gain loss of a signal with bandwidth $B_c$ Hz with a detection time offset of $\Delta t$ seconds
         is
 
-        $$\text{CGL} = -10 \log_{10} \left( \text{sinc}^2 \left( \Delta t B_C \right) \right) .$$
+        $$\text{CGL} = -10 \log_{10} \left( \text{sinc}^2 \left( \Delta t B_c \right) \right) .$$
 
     Examples:
         Compute the coherent gain loss for an integration time of 1 ms and a frequency offset of 235 Hz.

src/sdr/_detection/_theory.py

@@ -11,6 +11,7 @@
 import scipy.stats
 from typing_extensions import Literal
 
+from .._conversion import db as to_db
 from .._conversion import linear
 from .._helper import (
     convert_output,
@@ -21,6 +22,7 @@
     verify_not_specified,
     verify_scalar,
 )
+from .._probability import sum_distribution
 
 
 @export
@@ -242,7 +244,7 @@ def _h0(
             h0 = scipy.stats.norm(0, np.sqrt(sigma2_per))
         elif detector == "linear":
             h0 = scipy.stats.chi(nu, scale=np.sqrt(sigma2_per))
-            h0 = _sum_distribution(h0, n_nc)
+            h0 = sum_distribution(h0, n_nc)
         elif detector == "square-law":
             h0 = scipy.stats.chi2(nu * n_nc, scale=sigma2_per)
 
@@ -488,58 +490,13 @@ def _h1(
             else:
                 # Folded normal distribution has 1 degree of freedom
                 h1 = scipy.stats.foldnorm(np.sqrt(lambda_), scale=np.sqrt(sigma2_per))
-            h1 = _sum_distribution(h1, n_nc)
+            h1 = sum_distribution(h1, n_nc)
         elif detector == "square-law":
             h1 = scipy.stats.ncx2(nu * n_nc, lambda_ * n_nc, scale=sigma2_per)
 
     return h1
 
 
-def _sum_distribution(
-    dist: scipy.stats.rv_continuous, n_nc: int
-) -> scipy.stats.rv_histogram | scipy.stats.rv_continuous:
-    r"""
-    Sums a distribution `n_nc` times.
-
-    Arguments:
-        dist: The distribution to sum.
-        n_nc: The number of times to sum the distribution.
-
-    Returns:
-        The summed distribution.
-    """
-    if n_nc == 1:
-        return dist
-
-    # Determine mean and standard deviation of base distribution
-    mu, sigma2 = dist.stats()
-    sigma = np.sqrt(sigma2)
-
-    # NOTE: I was only able to get this to work with x starting at 0. When the x axis start below zero,
-    # I couldn't get the correct offset for the convolved x axis.
-
-    # Compute the PDF of the base distribution for 10 standard deviations about the mean
-    pdf_x = np.linspace(0, mu + 10 * sigma, 1_001)
-    pdf_y = dist.pdf(pdf_x)
-    dx = np.mean(np.diff(pdf_x))
-
-    # The PDF of the sum of n_nc independent random variables is the convolution of the PDF of the base distribution.
-    # This is efficiently computed in the frequency domain by exponentiating the FFT. The FFT must be zero-padded
-    # enough that the circular convolutions do not wrap around.
-    n_fft = scipy.fft.next_fast_len(pdf_y.size * n_nc)
-    Y = np.fft.fft(pdf_y, n_fft)
-    Y = Y**n_nc
-    y = np.fft.ifft(Y).real
-    y /= y.sum() * dx
-    x = np.arange(y.size) * dx + pdf_x[0]
-
-    # Adjust the histograms bins to be on either side of each point. So there is one extra point added.
-    x = np.append(x, x[-1] + dx)
-    x -= dx / 2
-
-    return scipy.stats.rv_histogram((y, x))
-
-
 @export
 def p_d(
     snr: npt.ArrayLike,
@@ -915,6 +872,82 @@ def _calculate(p_fa, sigma2):
     return convert_output(threshold)
 
 
+@export
+def threshold_factor(
+    p_fa: npt.ArrayLike,
+    detector: Literal["coherent", "linear", "square-law"] = "square-law",
+    complex: bool = True,
+    n_c: int = 1,
+    n_nc: int | None = None,
+    db: bool = False,
+) -> npt.NDArray[np.float64]:
+    r"""
+    Computes the theoretical detection threshold factor $\alpha$.
+
+    Arguments:
+        p_fa: The desired probability of false alarm $P_{fa}$ in $(0, 1)$.
+        detector: The detector type.
+
+            - `"coherent"`: A coherent detector, $$T(x) = \mathrm{Re}\left\{\sum_{i=0}^{N_c-1} x[n-i]\right\} .$$
+            - `"linear"`: A linear detector, $$T(x) = \sum_{j=0}^{N_{nc}-1}\left|\sum_{i=0}^{N_c-1} x[n-i-jN_c]\right| .$$
+            - `"square-law"`: A square-law detector, $$T(x) = \sum_{j=0}^{N_{nc}-1}\left|\sum_{i=0}^{N_c-1} x[n-i-jN_c]\right|^2 .$$
+
+        complex: Indicates whether the input signal is real or complex. This affects how the SNR is converted
+            to noise variance.
+        n_c: The number of samples to coherently integrate $N_c$.
+        n_nc: The number of samples to non-coherently integrate $N_{nc}$. Non-coherent integration is only allowable
+            for linear and square-law detectors.
+        db: Indicates whether to return the detection threshold $\alpha$ in dB.
+
+    Returns:
+        The detection threshold factor $\alpha$.
+
+    Notes:
+        The detection threshold factor $\alpha$ is defined as the ratio of the detection threshold $\gamma$ to the
+        mean of the detector output under the null hypothesis. This is true for linear and square-law detectors.
+
+        $$\alpha = \frac{\gamma}{\frac{1}{N} \sum_{i=1}^{N}\{T(x[i]) \mid \mathcal{H}_0\}}$$
+
+        For coherent detectors, the detection threshold factor $\alpha$ is defined as the ratio of the detection
+        threshold $\gamma$ to the mean of the square of the detector output under the null hypothesis.
+        This is required because the mean of the coherent detector output is zero.
+
+        $$\alpha = \frac{\gamma}{\frac{1}{N} \sum_{i=1}^{N}\{T(x[i])^2 \mid \mathcal{H}_0\}}$$
+
+    Examples:
+        .. ipython:: python
+
+            p_fa = np.logspace(-16, -1, 101)  # Probability of false alarm
+
+            @savefig sdr_threshold_factor_1.png
+            plt.figure(); \
+            plt.semilogx(p_fa, sdr.threshold_factor(p_fa, detector="coherent"), label="Coherent"); \
+            plt.semilogx(p_fa, sdr.threshold_factor(p_fa, detector="linear"), label="Linear"); \
+            plt.semilogx(p_fa, sdr.threshold_factor(p_fa, detector="square-law"), label="Square-Law"); \
+            plt.xlabel("Probability of false alarm, $P_{fa}$"); \
+            plt.ylabel(r"Detection threshold factor, $\alpha$"); \
+            plt.legend(title="Detector"); \
+            plt.title("Detection threshold factor across false alarm rate");
+
+    Group:
+        detection-theory
+    """
+    sigma2 = 1
+    gamma = threshold(p_fa, sigma2, detector, complex, n_c, n_nc)
+    null_hypothesis = h0(sigma2, detector, complex, n_c, n_nc)
+
+    if detector == "coherent":
+        # Since mean is zero, the variance is equivalent to the mean of the square
+        alpha = gamma / null_hypothesis.var()
+    else:
+        alpha = gamma / null_hypothesis.mean()
+
+    if db:
+        alpha = to_db(alpha)
+
+    return alpha
+
+
 @export
 def min_snr(
     p_d: npt.ArrayLike,

src/sdr/_filter/_applications.py

@@ -184,7 +184,7 @@ def __init__(
             window: The SciPy window definition. See :func:`scipy.signal.windows.get_window` for details.
                 If `None`, no window is applied.
             streaming: Indicates whether to use streaming mode. In streaming mode, previous inputs are
-                preserved between calls to :meth:`~sdr.Differentiator.__call__()`.
+                preserved between calls to :meth:`~sdr.Differentiator.__call__`.
 
         Examples:
             See the :ref:`fir-filters` example.
@@ -198,8 +198,7 @@ def __init__(
         h[order // 2] = 0
 
         if window is not None:
-            h_win = scipy.signal.windows.get_window(window, order + 1, fftbins=False)
-            h *= h_win
+            h *= scipy.signal.windows.get_window(window, order + 1, fftbins=False)
 
         super().__init__(h, streaming=streaming)
 

src/sdr/_filter/_delay.py

@@ -13,7 +13,7 @@
 from ._fir import FIR
 
 
-def _ideal_frac_delay(length: int, delay: float) -> npt.NDArray[np.float64]:
+def _ideal_fractional_delay(length: int, delay: float) -> npt.NDArray[np.float64]:
     """
     Returns the ideal fractional delay filter impulse response.
     """
@@ -23,7 +23,7 @@ def _ideal_frac_delay(length: int, delay: float) -> npt.NDArray[np.float64]:
 
 
 @export
-def design_frac_delay_fir(
+def fractional_delay_fir(
     length: int,
     delay: float,
 ) -> npt.NDArray[np.float64]:
@@ -50,38 +50,38 @@ def design_frac_delay_fir(
 
         .. ipython:: python
 
-            h_8 = sdr.design_frac_delay_fir(8, 0.25)
+            h_8 = sdr.fractional_delay_fir(8, 0.25)
 
-            @savefig sdr_design_frac_delay_fir_1.png
+            @savefig sdr_fractional_delay_fir_1.png
             plt.figure(); \
             sdr.plot.impulse_response(h_8);
 
-            @savefig sdr_design_frac_delay_fir_2.png
+            @savefig sdr_fractional_delay_fir_2.png
             plt.figure(); \
             sdr.plot.magnitude_response(h_8); \
             plt.ylim(-4, 1);
 
-            @savefig sdr_design_frac_delay_fir_3.png
+            @savefig sdr_fractional_delay_fir_3.png
             plt.figure(); \
             sdr.plot.group_delay(h_8);
 
         Compare the magnitude response and group delay of filters with different lengths.
 
         .. ipython:: python
 
-            h_16 = sdr.design_frac_delay_fir(16, 0.25); \
-            h_32 = sdr.design_frac_delay_fir(32, 0.25); \
-            h_64 = sdr.design_frac_delay_fir(64, 0.25)
+            h_16 = sdr.fractional_delay_fir(16, 0.25); \
+            h_32 = sdr.fractional_delay_fir(32, 0.25); \
+            h_64 = sdr.fractional_delay_fir(64, 0.25)
 
-            @savefig sdr_design_frac_delay_fir_4.png
+            @savefig sdr_fractional_delay_fir_4.png
             plt.figure(); \
             sdr.plot.magnitude_response(h_8, label="Length 8"); \
             sdr.plot.magnitude_response(h_16, label="Length 16"); \
             sdr.plot.magnitude_response(h_32, label="Length 32"); \
             sdr.plot.magnitude_response(h_64, label="Length 64"); \
             plt.ylim(-4, 1);
 
-            @savefig sdr_design_frac_delay_fir_5.png
+            @savefig sdr_fractional_delay_fir_5.png
             plt.figure(); \
             sdr.plot.group_delay(h_8, label="Length 8"); \
             sdr.plot.group_delay(h_16, label="Length 16"); \
@@ -95,7 +95,7 @@ def design_frac_delay_fir(
     verify_scalar(delay, float=True, inclusive_min=0, inclusive_max=1)
 
     N = length - (length % 2)  # The length guaranteed to be even
-    h_ideal = _ideal_frac_delay(N, delay)
+    h_ideal = _ideal_fractional_delay(N, delay)
 
     if N == 2:
         beta = 0
@@ -181,7 +181,7 @@ def __init__(
         verify_scalar(length, int=True, inclusive_min=2)
         verify_scalar(delay, float=True, inclusive_min=0, inclusive_max=1)
 
-        h = design_frac_delay_fir(length, delay)
+        h = fractional_delay_fir(length, delay)
         super().__init__(h)
 
         self._delay = length // 2 - 1 + delay