From 01c9f96e2d0702ebb26de6d3c42a955706ae2248 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Mon, 14 Aug 2023 20:05:38 -0400 Subject: [PATCH 01/43] Set default `sample_rate` to `None` --- src/sdr/plot/_time_domain.py | 25 ++++++++++++++++--------- 1 file changed, 16 insertions(+), 9 deletions(-) diff --git a/src/sdr/plot/_time_domain.py b/src/sdr/plot/_time_domain.py index da8446ebf..fc10edea7 100644 --- a/src/sdr/plot/_time_domain.py +++ b/src/sdr/plot/_time_domain.py @@ -15,7 +15,7 @@ @export def time_domain( x: npt.ArrayLike, - sample_rate: float = 1.0, + sample_rate: float | None = None, centered: bool = False, offset: float = 0, diff: Literal["color", "line"] = "color", @@ -26,8 +26,8 @@ def time_domain( Arguments: x: The time-domain signal $x[n]$. - sample_rate: The sample rate $f_s$ of the signal in samples/s. If the sample rate is 1, the x-axis will - be label as "Samples". + sample_rate: The sample rate $f_s$ of the signal in samples/s. If `None`, the x-axis will + be labeled as "Samples". centered: Indicates whether to center the x-axis about 0. This argument is mutually exclusive with `offset`. offset: The x-axis offset to apply to the first sample. The units of the offset are $1/f_s$. @@ -76,10 +76,17 @@ def time_domain( Group: plot-time-domain """ - if not isinstance(sample_rate, (int, float)): - raise TypeError(f"Argument 'sample_rate' must be a number, not {type(sample_rate)}.") - x = np.asarray(x) + if not x.ndim == 1: + raise ValueError(f"Argument 'x' must be 1-D, not {x.ndim}-D.") + + if sample_rate is None: + sample_rate_provided = False + sample_rate = 1 + else: + sample_rate_provided = True + if not isinstance(sample_rate, (int, float)): + raise TypeError(f"Argument 'sample_rate' must be a number, not {type(sample_rate)}.") if centered: if x.size % 2 == 0: @@ -112,9 +119,9 @@ def time_domain( if label: plt.legend() - if sample_rate == 1: - plt.xlabel("Samples") - else: + if sample_rate_provided: plt.xlabel("Time (s)") + else: + plt.xlabel("Samples") plt.ylabel("Amplitude") plt.tight_layout() From a2b039231eb1987333ab3fe3333f3ab1f4db65b1 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Mon, 14 Aug 2023 20:29:34 -0400 Subject: [PATCH 02/43] Add `sdr.plot.raster()` --- src/sdr/plot/_time_domain.py | 101 +++++++++++++++++++++++++++++++++++ 1 file changed, 101 insertions(+) diff --git a/src/sdr/plot/_time_domain.py b/src/sdr/plot/_time_domain.py index fc10edea7..c2cb5c07c 100644 --- a/src/sdr/plot/_time_domain.py +++ b/src/sdr/plot/_time_domain.py @@ -6,6 +6,7 @@ import matplotlib.pyplot as plt import numpy as np import numpy.typing as npt +from matplotlib.collections import LineCollection from typing_extensions import Literal from .._helper import export @@ -125,3 +126,103 @@ def time_domain( plt.xlabel("Samples") plt.ylabel("Amplitude") plt.tight_layout() + + +@export +def raster( + x: npt.ArrayLike, + length: int, + stride: int | None = None, + sample_rate: float | None = None, + color: Literal["index"] | str = "index", + colorbar: bool = False, + **kwargs, +): + """ + Plots a raster of the time-domain signal $x[n]$. + + Arguments: + x: The time-domain signal $x[n]$. + length: The length of each raster in samples. + stride: The stride between each raster in samples. If `None`, the stride is set to `length`. + sample_rate: The sample rate $f_s$ of the signal in samples/s. If `None`, the x-axis will + be labeled as "Samples". + color: Indicates how to color the rasters. If `"index"`, the rasters are colored based on their index. + If a valid Matplotlib color, the rasters are all colored with that color. + colorbar: Indicates whether to add a colorbar to the plot. + kwargs: Additional keyword arguments to pass to :obj:`matplotlib.collections.LineCollection`. + The following keyword arguments are set by default. The defaults may be overwritten. + + - `"linewidths"`: `(0.5, 1, 1.5, 2)` + - `"linestyles"`: `"solid"` + + Group: + plot-time-domain + """ + x = np.asarray(x) + if not np.isrealobj(x): + raise TypeError(f"Argument 'x' must be real, not {x.dtype}.") + if not x.ndim == 1: + raise ValueError(f"Argument 'x' must be 1-D, not {x.ndim}-D.") + + if not isinstance(length, int): + raise TypeError(f"Argument 'length' must be an integer, not {type(length)}.") + if not 1 <= length <= x.size: + raise ValueError(f"Argument 'length' must be at least 1 and less than the length of 'x', not {length}.") + + if stride is None: + stride = length + elif not isinstance(stride, int): + raise TypeError(f"Argument 'stride' must be an integer, not {type(stride)}.") + elif not 1 <= stride <= x.size: + raise ValueError(f"Argument 'stride' must be at least 1 and less than the length of 'x', not {stride}.") + + if sample_rate is None: + sample_rate_provided = False + t = np.arange(length) + else: + sample_rate_provided = True + if not isinstance(sample_rate, (int, float)): + raise TypeError(f"Argument 'sample_rate' must be a number, not {type(sample_rate)}.") + t = np.arange(length) / sample_rate + + # Compute the strided data and format into segments for LineCollection + N_rasters = (x.size - length) // stride + 1 + x_strided = np.lib.stride_tricks.as_strided( + x, shape=(N_rasters, length), strides=(x.strides[0] * stride, x.strides[0]), writeable=False + ) + segs = [np.column_stack([t, x_raster]) for x_raster in x_strided] + + # Set the default keyword arguments and override with user-specified keyword arguments + default_kwargs = { + "linewidths": (0.5, 1, 1.5, 2), + "linestyles": "solid", + } + if color == "index": + default_kwargs["array"] = np.arange(N_rasters) + else: + default_kwargs["colors"] = color + kwargs = {**default_kwargs, **kwargs} + + line_collection = LineCollection( + segs, + **kwargs, + ) + + with plt.rc_context(RC_PARAMS): + ax = plt.gca() + ax.add_collection(line_collection) + ax.set_xlim(t.min(), t.max()) + ax.set_ylim(x.min(), x.max()) + + if colorbar: + axcb = plt.colorbar(line_collection) + axcb.set_label("Raster Index") + + plt.grid(True) + if sample_rate_provided: + plt.xlabel("Time (s)") + else: + plt.xlabel("Samples") + plt.ylabel("Amplitude") + plt.tight_layout() From 46ff0f588e9afb80824622bc9d8509015af51942 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Mon, 14 Aug 2023 20:55:38 -0400 Subject: [PATCH 03/43] Add `sdr.plot.eye()` --- src/sdr/plot/_modulation.py | 69 +++++++++++++++++++++++++++++++++++++ 1 file changed, 69 insertions(+) diff --git a/src/sdr/plot/_modulation.py b/src/sdr/plot/_modulation.py index 149956bd4..30d30baee 100644 --- a/src/sdr/plot/_modulation.py +++ b/src/sdr/plot/_modulation.py @@ -10,6 +10,7 @@ from .._helper import export from ._rc_params import RC_PARAMS +from ._time_domain import raster @export @@ -182,6 +183,74 @@ def symbol_map( plt.tight_layout() +@export +def eye( + x: npt.ArrayLike, + sps: int, + span: int = 2, + sample_rate: float | None = None, + color: Literal["index"] | str = "index", + **kwargs, +): + r""" + Plots the eye diagram of the real baseband signal $x[k]$. + + Arguments: + x: The real baseband signal $x[k]$. + sps: The number of samples per symbol. + span: The number of symbols per raster. + sample_rate: The sample rate $f_s$ of the signal in samples/s. If `None`, the x-axis will + be labeled as "Samples". + color: Indicates how to color the rasters. If `"index"`, the rasters are colored based on their index. + If a valid Matplotlib color, the rasters are all colored with that color. + kwargs: Additional keyword arguments to pass to :func:`sdr.plot.raster()`. + + Note: + To plot an eye diagram for I and Q, call this function twice, passing `x.real` and then `x.imag`. + + Example: + Modulate 100 BPSK symbols. + + .. ipython:: python + + psk = sdr.PSK(2); \ + s = np.random.randint(0, psk.order, 100); \ + a = psk.modulate(s) + + Apply a raised cosine pulse shape and examine the eye diagram of the I channel. Since the raised + cosine pulse shape is a Nyquist filter, there is no intersymbol interference (ISI) at the symbol decisions. + + .. ipython:: python + + sps = 25; \ + h = sdr.raised_cosine(0.5, 6, sps); \ + fir = sdr.Interpolator(sps, h); \ + x = fir(a) + + @savefig sdr_plot_eye_1.png + plt.figure(figsize=(8, 4)); \ + sdr.plot.eye(x.real, sps) + + Apply a root raised cosine pulse shape and examine the eye diagram of the I channel. The root raised + cosine filter is not a Nyquist filter, and ISI can be observed. + + .. ipython:: python + + sps = 25; \ + h = sdr.root_raised_cosine(0.5, 6, sps); \ + fir = sdr.Interpolator(sps, h); \ + x = fir(a) + + @savefig sdr_plot_eye_2.png + plt.figure(figsize=(8, 4)); \ + sdr.plot.eye(x.real, sps) + + Group: + plot-modulation + """ + raster(x, span * sps + 1, stride=sps, sample_rate=sample_rate, color=color, **kwargs) + + @export def ber( ebn0: npt.ArrayLike, From b47b2b41ca0442d68bc77a204fd16fcc825a7cc6 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Mon, 14 Aug 2023 21:00:12 -0400 Subject: [PATCH 04/43] Update feature list --- README.md | 5 +++-- docs/index.rst | 5 +++-- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/README.md b/README.md index fa444b36f..82856cd43 100644 --- a/README.md +++ b/README.md @@ -62,8 +62,9 @@ View all available classes and functions in the [API Reference](https://mhostett additive white Gaussian noise (AWGN), frequency offset, sample rate offset, IQ imbalance. - **Link budgets**: Channel capacities, free-space path loss, antenna gains. - **Data manipulation**: Packing and unpacking binary data, hexdumping binary data. -- **Plotting**: Time-domain, periodogram, spectrogram, BER, SER, constellation, symbol map, impulse response, - step response, magnitude response, phase response, phase delay, group delay, and zeros/poles. +- **Plotting**: Time-domain, raster, periodogram, spectrogram, constellation, symbol map, eye diagram, + bit error rate (BER), symbol error rate (SER), impulse response, step response, magnitude response, phase response, + phase delay, group delay, and zeros/poles. ## Examples diff --git a/docs/index.rst b/docs/index.rst index c5d5023cc..42a05eea2 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -53,8 +53,9 @@ View all available classes and functions in the `API Reference Date: Mon, 14 Aug 2023 21:15:20 -0400 Subject: [PATCH 05/43] Fix typo --- src/sdr/plot/_modulation.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/sdr/plot/_modulation.py b/src/sdr/plot/_modulation.py index 30d30baee..14bd3616f 100644 --- a/src/sdr/plot/_modulation.py +++ b/src/sdr/plot/_modulation.py @@ -193,10 +193,10 @@ def eye( **kwargs, ): r""" - Plots the eye diagram of the real baseband signal $x[k]$. + Plots the eye diagram of the real baseband signal $x[n]$. Arguments: - x: The real baseband signal $x[k]$. + x: The real baseband signal $x[n]$. sps: The number of samples per symbol. span: The number of symbols per raster. sample_rate: The sample rate $f_s$ of the signal in samples/s. If `None`, the x-axis will From 1c2a25f595738c63bf68dd79ac7858f433d9de88 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Mon, 14 Aug 2023 22:35:21 -0400 Subject: [PATCH 06/43] Fix raster y limits --- src/sdr/plot/_time_domain.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/sdr/plot/_time_domain.py b/src/sdr/plot/_time_domain.py index c2cb5c07c..d2040046d 100644 --- a/src/sdr/plot/_time_domain.py +++ b/src/sdr/plot/_time_domain.py @@ -213,7 +213,7 @@ def raster( ax = plt.gca() ax.add_collection(line_collection) ax.set_xlim(t.min(), t.max()) - ax.set_ylim(x.min(), x.max()) + ax.set_ylim(-np.abs(x).max(), np.abs(x).max()) if colorbar: axcb = plt.colorbar(line_collection) From 94d9d729a72dda311e4c4cf80f2f5cd8ee12247c Mon Sep 17 00:00:00 2001 From: mhostetter Date: Mon, 14 Aug 2023 22:51:22 -0400 Subject: [PATCH 07/43] Fix line widths on `raster()` plot --- src/sdr/plot/_time_domain.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/sdr/plot/_time_domain.py b/src/sdr/plot/_time_domain.py index d2040046d..4d9afb295 100644 --- a/src/sdr/plot/_time_domain.py +++ b/src/sdr/plot/_time_domain.py @@ -153,7 +153,7 @@ def raster( kwargs: Additional keyword arguments to pass to :obj:`matplotlib.collections.LineCollection`. The following keyword arguments are set by default. The defaults may be overwritten. - - `"linewidths"`: `(0.5, 1, 1.5, 2)` + - `"linewidths"`: 1 - `"linestyles"`: `"solid"` Group: @@ -195,7 +195,7 @@ def raster( # Set the default keyword arguments and override with user-specified keyword arguments default_kwargs = { - "linewidths": (0.5, 1, 1.5, 2), + "linewidths": 1, "linestyles": "solid", } if color == "index": From 2a1e1dcc4f9bd4078777686f3667e5b83ba64375 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Mon, 14 Aug 2023 23:11:04 -0400 Subject: [PATCH 08/43] Support raster of complex data --- src/sdr/plot/_time_domain.py | 19 ++++++++++++++----- 1 file changed, 14 insertions(+), 5 deletions(-) diff --git a/src/sdr/plot/_time_domain.py b/src/sdr/plot/_time_domain.py index 4d9afb295..3908c5a08 100644 --- a/src/sdr/plot/_time_domain.py +++ b/src/sdr/plot/_time_domain.py @@ -142,7 +142,8 @@ def raster( Plots a raster of the time-domain signal $x[n]$. Arguments: - x: The time-domain signal $x[n]$. + x: The time-domain signal $x[n]$. If `x` is complex, the real and imaginary rasters are interleaved. + Time order is preserved. length: The length of each raster in samples. stride: The stride between each raster in samples. If `None`, the stride is set to `length`. sample_rate: The sample rate $f_s$ of the signal in samples/s. If `None`, the x-axis will @@ -160,8 +161,6 @@ def raster( plot-time-domain """ x = np.asarray(x) - if not np.isrealobj(x): - raise TypeError(f"Argument 'x' must be real, not {x.dtype}.") if not x.ndim == 1: raise ValueError(f"Argument 'x' must be 1-D, not {x.ndim}-D.") @@ -191,7 +190,17 @@ def raster( x_strided = np.lib.stride_tricks.as_strided( x, shape=(N_rasters, length), strides=(x.strides[0] * stride, x.strides[0]), writeable=False ) - segs = [np.column_stack([t, x_raster]) for x_raster in x_strided] + + if np.iscomplexobj(x): + segments_real = [np.column_stack([t, x_raster.real]) for x_raster in x_strided] + segments_imag = [np.column_stack([t, x_raster.imag]) for x_raster in x_strided] + + # Interleave the real and imaginary rasters + segments = [None] * (2 * N_rasters) + segments[::2] = segments_real + segments[1::2] = segments_imag + else: + segments = [np.column_stack([t, x_raster]) for x_raster in x_strided] # Set the default keyword arguments and override with user-specified keyword arguments default_kwargs = { @@ -205,7 +214,7 @@ def raster( kwargs = {**default_kwargs, **kwargs} line_collection = LineCollection( - segs, + segments, **kwargs, ) From 38e24f6deda1ff3130002c72621c428412978da5 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Mon, 14 Aug 2023 23:22:00 -0400 Subject: [PATCH 09/43] Support eye diagrams of complex signals --- src/sdr/plot/_modulation.py | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/src/sdr/plot/_modulation.py b/src/sdr/plot/_modulation.py index 14bd3616f..d716d2375 100644 --- a/src/sdr/plot/_modulation.py +++ b/src/sdr/plot/_modulation.py @@ -193,10 +193,11 @@ def eye( **kwargs, ): r""" - Plots the eye diagram of the real baseband signal $x[n]$. + Plots the eye diagram of the baseband modulated signal $x[n]$. Arguments: - x: The real baseband signal $x[n]$. + x: The baseband modulated signal $x[n]$. If `x` is complex, the real and imaginary rasters are interleaved. + Time order is preserved. sps: The number of samples per symbol. span: The number of symbols per raster. sample_rate: The sample rate $f_s$ of the signal in samples/s. If `None`, the x-axis will @@ -205,20 +206,17 @@ def eye( If a valid Matplotlib color, the rasters are all colored with that color. kwargs: Additional keyword arguments to pass to :func:`sdr.plot.raster()`. - Note: - To plot an eye diagram for I and Q, call this function twice, passing `x.real` and then `x.imag`. - Example: - Modulate 100 BPSK symbols. + Modulate 100 QPSK symbols. .. ipython:: python - psk = sdr.PSK(2); \ + psk = sdr.PSK(4, phase_offset=45); \ s = np.random.randint(0, psk.order, 100); \ a = psk.modulate(s) - Apply a raised cosine pulse shape and examine the eye diagram of the I channel. Since the raised - cosine pulse shape is a Nyquist filter, there is no intersymbol interference (ISI) at the symbol decisions. + Apply a raised cosine pulse shape and examine the eye diagram. Since the raised cosine pulse shape + is a Nyquist filter, there is no intersymbol interference (ISI) at the symbol decisions. .. ipython:: python @@ -229,10 +227,12 @@ def eye( @savefig sdr_plot_eye_1.png plt.figure(figsize=(8, 4)); \ - sdr.plot.eye(x.real, sps) + sdr.plot.eye(x, sps) - Apply a root raised cosine pulse shape and examine the eye diagram of the I channel. The root raised - cosine filter is not a Nyquist filter, and ISI can be observed. + Apply a root raised cosine pulse shape and examine the eye diagram. The root raised cosine filter + is not a Nyquist filter, and ISI can be observed. (It should be noted that two cascaded root raised + cosine filters, one for transmit and one for receive, is a Nyquist filter. This is why SRRC pulse shaping + is often used in practice.) .. ipython:: python @@ -243,7 +243,7 @@ def eye( @savefig sdr_plot_eye_2.png plt.figure(figsize=(8, 4)); \ - sdr.plot.eye(x.real, sps) + sdr.plot.eye(x, sps) Group: plot-modulation From cda5009ceb352ea621274c0352c2ced12952fa68 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Mon, 14 Aug 2023 23:22:20 -0400 Subject: [PATCH 10/43] Expand dictionary --- .vscode/settings.json | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/.vscode/settings.json b/.vscode/settings.json index 214339a72..2db195189 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -24,18 +24,23 @@ "downsampling", "exponentials", "Feedforward", + "figsize", "hexdump", + "intersymbol", "lowpass", "matplotlib", "multirate", + "Nyquist", "overdamped", "passband", "periodogram", "periodograms", "polyphase", "pyplot", + "QPSK", "resampler", "resamplers", + "savefig", "scipy", "sinc", "stopband", From 6c5b30300557b5e728631b892e1a744b8e4d13be Mon Sep 17 00:00:00 2001 From: mhostetter Date: Mon, 14 Aug 2023 23:24:52 -0400 Subject: [PATCH 11/43] Fix linter errors --- src/sdr/plot/_time_domain.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/sdr/plot/_time_domain.py b/src/sdr/plot/_time_domain.py index 3908c5a08..55485313d 100644 --- a/src/sdr/plot/_time_domain.py +++ b/src/sdr/plot/_time_domain.py @@ -160,6 +160,7 @@ def raster( Group: plot-time-domain """ + # pylint: disable=too-many-statements x = np.asarray(x) if not x.ndim == 1: raise ValueError(f"Argument 'x' must be 1-D, not {x.ndim}-D.") From d3128292ee7822ac9e5c3ef8692868c1463bb59c Mon Sep 17 00:00:00 2001 From: mhostetter Date: Tue, 15 Aug 2023 08:51:25 -0400 Subject: [PATCH 12/43] Tighten `raster()` y-limits --- src/sdr/plot/_time_domain.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/sdr/plot/_time_domain.py b/src/sdr/plot/_time_domain.py index 55485313d..4da0240c4 100644 --- a/src/sdr/plot/_time_domain.py +++ b/src/sdr/plot/_time_domain.py @@ -223,7 +223,7 @@ def raster( ax = plt.gca() ax.add_collection(line_collection) ax.set_xlim(t.min(), t.max()) - ax.set_ylim(-np.abs(x).max(), np.abs(x).max()) + ax.set_ylim(x.min(), x.max()) if colorbar: axcb = plt.colorbar(line_collection) From ffa8190082d667dc3e65c2fd23c8e4aca16d5119 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Tue, 15 Aug 2023 08:52:40 -0400 Subject: [PATCH 13/43] Make `eye()` y-axis symmetric --- src/sdr/plot/_modulation.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/sdr/plot/_modulation.py b/src/sdr/plot/_modulation.py index d716d2375..7d77807c1 100644 --- a/src/sdr/plot/_modulation.py +++ b/src/sdr/plot/_modulation.py @@ -250,6 +250,12 @@ def eye( """ raster(x, span * sps + 1, stride=sps, sample_rate=sample_rate, color=color, **kwargs) + # Make y-axis symmetric + ax = plt.gca() + ymin, ymax = ax.get_ylim() + ylim = max(np.abs(ymin), np.abs(ymax)) + ax.set_ylim(-ylim, ylim) + @export def ber( From 7ba38818dc5d4504d61044412a32b9aa379ec041 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Tue, 15 Aug 2023 08:55:57 -0400 Subject: [PATCH 14/43] Support 2D input to `raster()` --- src/sdr/plot/_time_domain.py | 56 +++++++++++++++++++++--------------- 1 file changed, 33 insertions(+), 23 deletions(-) diff --git a/src/sdr/plot/_time_domain.py b/src/sdr/plot/_time_domain.py index 4da0240c4..0aa11208e 100644 --- a/src/sdr/plot/_time_domain.py +++ b/src/sdr/plot/_time_domain.py @@ -131,7 +131,7 @@ def time_domain( @export def raster( x: npt.ArrayLike, - length: int, + length: int | None = None, stride: int | None = None, sample_rate: float | None = None, color: Literal["index"] | str = "index", @@ -143,8 +143,9 @@ def raster( Arguments: x: The time-domain signal $x[n]$. If `x` is complex, the real and imaginary rasters are interleaved. - Time order is preserved. - length: The length of each raster in samples. + Time order is preserved. If `x` is 1D, the rastering is determined by `length` and `stride`. + If `x` is 2D, the rows correspond to each raster. + length: The length of each raster in samples. This must be provided if `x` is 1D. stride: The stride between each raster in samples. If `None`, the stride is set to `length`. sample_rate: The sample rate $f_s$ of the signal in samples/s. If `None`, the x-axis will be labeled as "Samples". @@ -162,20 +163,35 @@ def raster( """ # pylint: disable=too-many-statements x = np.asarray(x) - if not x.ndim == 1: - raise ValueError(f"Argument 'x' must be 1-D, not {x.ndim}-D.") - - if not isinstance(length, int): - raise TypeError(f"Argument 'length' must be an integer, not {type(length)}.") - if not 1 <= length <= x.size: - raise ValueError(f"Argument 'length' must be at least 1 and less than the length of 'x', not {length}.") + if not x.ndim in [1, 2]: + raise ValueError(f"Argument 'x' must be 1-D or 2-D, not {x.ndim}-D.") + + if x.ndim == 1: + if not isinstance(length, int): + raise TypeError(f"Argument 'length' must be an integer, not {type(length)}.") + if not 1 <= length <= x.size: + raise ValueError(f"Argument 'length' must be at least 1 and less than the length of 'x', not {length}.") + + if stride is None: + stride = length + elif not isinstance(stride, int): + raise TypeError(f"Argument 'stride' must be an integer, not {type(stride)}.") + elif not 1 <= stride <= x.size: + raise ValueError(f"Argument 'stride' must be at least 1 and less than the length of 'x', not {stride}.") + + # Compute the strided data and format into segments for LineCollection + N_rasters = (x.size - length) // stride + 1 + x_strided = np.lib.stride_tricks.as_strided( + x, shape=(N_rasters, length), strides=(x.strides[0] * stride, x.strides[0]), writeable=False + ) + else: + if not length is None: + raise ValueError("Argument 'length' can not be specified if 'x' is 2-D.") + if not stride is None: + raise ValueError("Argument 'stride' can not be specified if 'x' is 2-D.") - if stride is None: - stride = length - elif not isinstance(stride, int): - raise TypeError(f"Argument 'stride' must be an integer, not {type(stride)}.") - elif not 1 <= stride <= x.size: - raise ValueError(f"Argument 'stride' must be at least 1 and less than the length of 'x', not {stride}.") + N_rasters, length = x.shape + x_strided = x if sample_rate is None: sample_rate_provided = False @@ -186,17 +202,11 @@ def raster( raise TypeError(f"Argument 'sample_rate' must be a number, not {type(sample_rate)}.") t = np.arange(length) / sample_rate - # Compute the strided data and format into segments for LineCollection - N_rasters = (x.size - length) // stride + 1 - x_strided = np.lib.stride_tricks.as_strided( - x, shape=(N_rasters, length), strides=(x.strides[0] * stride, x.strides[0]), writeable=False - ) - + # Interleave the real and imaginary rasters, if necessary if np.iscomplexobj(x): segments_real = [np.column_stack([t, x_raster.real]) for x_raster in x_strided] segments_imag = [np.column_stack([t, x_raster.imag]) for x_raster in x_strided] - # Interleave the real and imaginary rasters segments = [None] * (2 * N_rasters) segments[::2] = segments_real segments[1::2] = segments_imag From 50c23193e0300a48b62297d02289115cdc96699f Mon Sep 17 00:00:00 2001 From: mhostetter Date: Tue, 15 Aug 2023 12:20:22 -0400 Subject: [PATCH 15/43] Change default colormap for `raster()` --- src/sdr/plot/_time_domain.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/sdr/plot/_time_domain.py b/src/sdr/plot/_time_domain.py index 0aa11208e..cf6848833 100644 --- a/src/sdr/plot/_time_domain.py +++ b/src/sdr/plot/_time_domain.py @@ -157,6 +157,7 @@ def raster( - `"linewidths"`: 1 - `"linestyles"`: `"solid"` + - `"cmap"`: `"rainbow"` Group: plot-time-domain @@ -217,6 +218,7 @@ def raster( default_kwargs = { "linewidths": 1, "linestyles": "solid", + "cmap": "rainbow", } if color == "index": default_kwargs["array"] = np.arange(N_rasters) From aa897c8de2b6deefad130672b7e20d465fd1b5b7 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Tue, 15 Aug 2023 12:30:23 -0400 Subject: [PATCH 16/43] Only add colorbar with index-based coloring --- src/sdr/plot/_time_domain.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/sdr/plot/_time_domain.py b/src/sdr/plot/_time_domain.py index cf6848833..80ca3c923 100644 --- a/src/sdr/plot/_time_domain.py +++ b/src/sdr/plot/_time_domain.py @@ -151,7 +151,7 @@ def raster( be labeled as "Samples". color: Indicates how to color the rasters. If `"index"`, the rasters are colored based on their index. If a valid Matplotlib color, the rasters are all colored with that color. - colorbar: Indicates whether to add a colorbar to the plot. + colorbar: Indicates whether to add a colorbar to the plot. This is only added if `color="index"`. kwargs: Additional keyword arguments to pass to :obj:`matplotlib.collections.LineCollection`. The following keyword arguments are set by default. The defaults may be overwritten. @@ -168,6 +168,8 @@ def raster( raise ValueError(f"Argument 'x' must be 1-D or 2-D, not {x.ndim}-D.") if x.ndim == 1: + if not length is not None: + raise ValueError("Argument 'length' must be specified if 'x' is 1-D.") if not isinstance(length, int): raise TypeError(f"Argument 'length' must be an integer, not {type(length)}.") if not 1 <= length <= x.size: @@ -237,7 +239,7 @@ def raster( ax.set_xlim(t.min(), t.max()) ax.set_ylim(x.min(), x.max()) - if colorbar: + if colorbar and color == "index": axcb = plt.colorbar(line_collection) axcb.set_label("Raster Index") From 090f61aadc1c64408d4de3cc3a9bd332d3fecfa2 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Tue, 15 Aug 2023 12:30:37 -0400 Subject: [PATCH 17/43] Enable colorbar by default --- src/sdr/plot/_time_domain.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/sdr/plot/_time_domain.py b/src/sdr/plot/_time_domain.py index 80ca3c923..d5fcb5876 100644 --- a/src/sdr/plot/_time_domain.py +++ b/src/sdr/plot/_time_domain.py @@ -135,7 +135,7 @@ def raster( stride: int | None = None, sample_rate: float | None = None, color: Literal["index"] | str = "index", - colorbar: bool = False, + colorbar: bool = True, **kwargs, ): """ From 9fee01eabab5fe183ed09ee975b31346fe22d77e Mon Sep 17 00:00:00 2001 From: mhostetter Date: Tue, 15 Aug 2023 21:20:18 -0400 Subject: [PATCH 18/43] Use system default float precision --- src/sdr/_farrow.py | 2 +- src/sdr/_filter/_fir.py | 4 ++-- src/sdr/_filter/_iir.py | 4 ++-- src/sdr/_filter/_polyphase.py | 4 ++-- src/sdr/_modulation/_pulse_shapes.py | 18 ++++++------------ src/sdr/_sequence.py | 2 +- src/sdr/plot/_filter.py | 4 ++-- tests/data_manipulation/test_hexdump.py | 2 +- tests/data_manipulation/test_pack.py | 2 +- tests/data_manipulation/test_unpack.py | 2 +- tests/dsp/test_iir.py | 2 +- tests/dsp/test_to_complex_bb.py | 2 +- tests/dsp/test_to_real_pb.py | 2 +- tests/simulation/test_awgn.py | 4 ++-- 14 files changed, 24 insertions(+), 30 deletions(-) diff --git a/src/sdr/_farrow.py b/src/sdr/_farrow.py index 6612bce77..874098d8f 100644 --- a/src/sdr/_farrow.py +++ b/src/sdr/_farrow.py @@ -48,7 +48,7 @@ def __init__(self, streaming: bool = False): [1 / 6, -1, 1 / 2, 1 / 3], [0, 0, 1, 0], ], - dtype=np.float32, + dtype=float, ) self.reset() diff --git a/src/sdr/_filter/_fir.py b/src/sdr/_filter/_fir.py index f0e549f6c..8ec30bde6 100644 --- a/src/sdr/_filter/_fir.py +++ b/src/sdr/_filter/_fir.py @@ -224,7 +224,7 @@ def impulse_response(self, N: int | None = None) -> np.ndarray: raise ValueError("Argument 'N' must be greater than or equal to the filter length.") # Delta impulse function - d = np.zeros(N - self.taps.size + 1, dtype=np.float32) + d = np.zeros(N - self.taps.size + 1, dtype=float) d[0] = 1 h = scipy.signal.convolve(d, self.taps, mode="full") @@ -255,7 +255,7 @@ def step_response(self, N: int | None = None) -> np.ndarray: raise ValueError("Argument 'N' must be greater than or equal to the filter length.") # Unit step function - u = np.ones(N - self.taps.size + 1, dtype=np.float32) + u = np.ones(N - self.taps.size + 1, dtype=float) s = scipy.signal.convolve(u, self.taps, mode="full") diff --git a/src/sdr/_filter/_iir.py b/src/sdr/_filter/_iir.py index 9f8a64813..5ad3999e5 100644 --- a/src/sdr/_filter/_iir.py +++ b/src/sdr/_filter/_iir.py @@ -215,7 +215,7 @@ def impulse_response(self, N: int = 100) -> np.ndarray: See the :ref:`iir-filters` example. """ # Delta impulse function - d = np.zeros(N, dtype=np.float32) + d = np.zeros(N, dtype=float) d[0] = 1 zi = self._state @@ -242,7 +242,7 @@ def step_response(self, N: int = 100) -> np.ndarray: See the :ref:`iir-filters` example. """ # Unit step function - u = np.ones(N, dtype=np.float32) + u = np.ones(N, dtype=float) state = self._state s = self(u) diff --git a/src/sdr/_filter/_polyphase.py b/src/sdr/_filter/_polyphase.py index 19967e1be..ec0a201cd 100644 --- a/src/sdr/_filter/_polyphase.py +++ b/src/sdr/_filter/_polyphase.py @@ -327,12 +327,12 @@ def __init__( taps = multirate_taps(rate, 1) elif taps == "linear": self._method = "linear" - taps = np.zeros(2 * rate, dtype=np.float32) + taps = np.zeros(2 * rate, dtype=float) taps[:rate] = np.arange(0, rate) / rate taps[rate:] = np.arange(rate, 0, -1) / rate elif taps == "zoh": self._method = "zoh" - taps = np.ones(rate, dtype=np.float32) + taps = np.ones(rate, dtype=float) else: raise ValueError(f"Argument 'taps' must be 'kaiser', 'linear', 'zoh', or an array-like, not {taps}.") diff --git a/src/sdr/_modulation/_pulse_shapes.py b/src/sdr/_modulation/_pulse_shapes.py index 6aab9555a..c305b6ec9 100644 --- a/src/sdr/_modulation/_pulse_shapes.py +++ b/src/sdr/_modulation/_pulse_shapes.py @@ -82,9 +82,7 @@ def raised_cosine(alpha: float, span: int, sps: int) -> np.ndarray: if not span * sps % 2 == 0: raise ValueError("The order of the filter (span * sps) must be even.") - # NOTE: Do the math in float64 to avoid numerical issues and then convert to float32 at the end - - t = np.arange(-(span * sps) // 2, (span * sps) // 2 + 1, dtype=np.float64) + t = np.arange(-(span * sps) // 2, (span * sps) // 2 + 1, dtype=float) Ts = sps # Handle special cases where the denominator is zero @@ -102,7 +100,7 @@ def raised_cosine(alpha: float, span: int, sps: int) -> np.ndarray: # Make the filter have unit energy h /= np.sqrt(np.sum(np.abs(h) ** 2)) - return h.astype(np.float32) + return h @export @@ -179,9 +177,7 @@ def root_raised_cosine(alpha: float, span: int, sps: int) -> np.ndarray: if not span * sps % 2 == 0: raise ValueError("The order of the filter (span * sps) must be even.") - # NOTE: Do the math in float64 to avoid numerical issues and then convert to float32 at the end - - t = np.arange(-(span * sps) // 2, (span * sps) // 2 + 1, dtype=np.float64) + t = np.arange(-(span * sps) // 2, (span * sps) // 2 + 1, dtype=float) Ts = sps # Symbol duration (in samples) # Handle special cases where the denominator is zero @@ -201,7 +197,7 @@ def root_raised_cosine(alpha: float, span: int, sps: int) -> np.ndarray: # Make the filter have unit energy h /= np.sqrt(np.sum(np.abs(h) ** 2)) - return h.astype(np.float32) + return h @export @@ -265,9 +261,7 @@ def gaussian(time_bandwidth: float, span: int, sps: int) -> np.ndarray: if not span * sps % 2 == 0: raise ValueError("The order of the filter (span * sps) must be even.") - # NOTE: Do the math in float64 to avoid numerical issues and then convert to float32 at the end - - t = np.arange(-(span * sps) // 2, (span * sps) // 2 + 1, dtype=np.float64) + t = np.arange(-(span * sps) // 2, (span * sps) // 2 + 1, dtype=float) t /= sps # Equation B.2 @@ -280,4 +274,4 @@ def gaussian(time_bandwidth: float, span: int, sps: int) -> np.ndarray: # Normalize coefficients so passband gain is 1 h = h / np.sum(h) - return h.astype(np.float32) + return h diff --git a/src/sdr/_sequence.py b/src/sdr/_sequence.py index 7f870c88d..71d003119 100644 --- a/src/sdr/_sequence.py +++ b/src/sdr/_sequence.py @@ -79,7 +79,7 @@ def barker(length: int, output: Literal["binary", "bipolar"] = "bipolar") -> np. # Map binary Barker code to bipolar. The mapping is equivalent to BPSK: 0 -> 1, 1 -> -1. sequence = 1 - 2 * code - sequence = sequence.astype(np.float64) + sequence = sequence.astype(float) return sequence diff --git a/src/sdr/plot/_filter.py b/src/sdr/plot/_filter.py index c587020a0..ce266d65b 100644 --- a/src/sdr/plot/_filter.py +++ b/src/sdr/plot/_filter.py @@ -91,7 +91,7 @@ def impulse_response( N = 100 # IIR filter # Delta impulse function - d = np.zeros(N, dtype=np.float32) + d = np.zeros(N, dtype=float) d[0] = 1 # Filter the impulse @@ -177,7 +177,7 @@ def step_response( N = 100 # IIR filter # Unit step function - u = np.ones(N, dtype=np.float32) + u = np.ones(N, dtype=float) # Filter the impulse zi = scipy.signal.lfiltic(b, a, y=[], x=[]) diff --git a/tests/data_manipulation/test_hexdump.py b/tests/data_manipulation/test_hexdump.py index 7a6d2decc..1ed3d8c44 100644 --- a/tests/data_manipulation/test_hexdump.py +++ b/tests/data_manipulation/test_hexdump.py @@ -7,7 +7,7 @@ def test_exceptions(): with pytest.raises(ValueError): # Must only have integer arrays - sdr.hexdump(np.array([8, 13, 12], dtype=np.float32)) + sdr.hexdump(np.array([8, 13, 12], dtype=float)) with pytest.raises(ValueError): # Must only have 1D arrays sdr.hexdump(np.array([[8, 13, 12]], dtype=int)) diff --git a/tests/data_manipulation/test_pack.py b/tests/data_manipulation/test_pack.py index d97040951..d054bdb6e 100644 --- a/tests/data_manipulation/test_pack.py +++ b/tests/data_manipulation/test_pack.py @@ -6,7 +6,7 @@ def test_exceptions(): with pytest.raises(ValueError): - x = np.array([1, 0, 0, 0, 1, 1, 0, 1], dtype=np.float32) + x = np.array([1, 0, 0, 0, 1, 1, 0, 1], dtype=float) sdr.pack(x, 2) with pytest.raises(ValueError): x = np.array([1, 0, 0, 0, 1, 1, 0, 1], dtype=int) diff --git a/tests/data_manipulation/test_unpack.py b/tests/data_manipulation/test_unpack.py index 6a3924c62..b59e77e88 100644 --- a/tests/data_manipulation/test_unpack.py +++ b/tests/data_manipulation/test_unpack.py @@ -6,7 +6,7 @@ def test_exceptions(): with pytest.raises(ValueError): - x = np.array([8, 13, 12], dtype=np.float32) + x = np.array([8, 13, 12], dtype=float) sdr.pack(x, 4) with pytest.raises(ValueError): x = np.array([8, 13, 12], dtype=int) diff --git a/tests/dsp/test_iir.py b/tests/dsp/test_iir.py index 99a5ab8cd..c35550155 100644 --- a/tests/dsp/test_iir.py +++ b/tests/dsp/test_iir.py @@ -13,7 +13,7 @@ def test_single_pole_impulse_response(pole): - R. Lyons, Understanding Digital Signal Processing 3rd Edition, Section 6.3.1. """ N = 100 - x = np.zeros(N, dtype=np.float32) + x = np.zeros(N, dtype=float) x[0] = 1 b = np.array([1], dtype=np.complex64) diff --git a/tests/dsp/test_to_complex_bb.py b/tests/dsp/test_to_complex_bb.py index 98c81742a..20f484b2f 100644 --- a/tests/dsp/test_to_complex_bb.py +++ b/tests/dsp/test_to_complex_bb.py @@ -7,7 +7,7 @@ def test_exceptions(): with pytest.raises(ValueError): # x_r must be 1D - x_r = np.zeros((2, 2), dtype=np.float32) + x_r = np.zeros((2, 2), dtype=float) sdr.to_complex_bb(x_r) with pytest.raises(ValueError): # x_r must be real diff --git a/tests/dsp/test_to_real_pb.py b/tests/dsp/test_to_real_pb.py index ab2585cfa..4bb04b7c9 100644 --- a/tests/dsp/test_to_real_pb.py +++ b/tests/dsp/test_to_real_pb.py @@ -11,5 +11,5 @@ def test_exceptions(): sdr.to_real_pb(x_c) with pytest.raises(ValueError): # x_c must be real - x_c = np.zeros(4, dtype=np.float32) + x_c = np.zeros(4, dtype=float) sdr.to_real_pb(x_c) diff --git a/tests/simulation/test_awgn.py b/tests/simulation/test_awgn.py index 9fc628dd0..bd02dd3b0 100644 --- a/tests/simulation/test_awgn.py +++ b/tests/simulation/test_awgn.py @@ -13,7 +13,7 @@ def test_exceptions(): def test_real_from_snr(): A = np.random.uniform(2, 3) # Signal amplitude N = 1000 - x = A * np.ones(N, dtype=np.float32) # Use a constant signal so the power measurement has no uncertainty + x = A * np.ones(N, dtype=float) # Use a constant signal so the power measurement has no uncertainty snr = np.random.uniform(1, 2) # dB snr_linear = 10 ** (snr / 10) @@ -47,7 +47,7 @@ def test_complex_from_snr(): def test_real_from_noise(): A = np.random.uniform(2, 3) # Signal amplitude N = 1000 - x = A * np.ones(N, dtype=np.float32) # Use a constant signal so the power measurement has no uncertainty + x = A * np.ones(N, dtype=float) # Use a constant signal so the power measurement has no uncertainty noise_std = np.random.uniform(1, 2) y = sdr.awgn(x, noise=noise_std**2) From 81a45a0da11687e4414a3d6b2358c9899bfb6172 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Tue, 15 Aug 2023 21:21:52 -0400 Subject: [PATCH 19/43] Use system default complex float precision --- docs/examples/phase-locked-loop.ipynb | 2 +- src/sdr/plot/_time_domain.py | 2 +- tests/dsp/test_decimator.py | 2 +- tests/dsp/test_fir.py | 2 +- tests/dsp/test_iir.py | 4 ++-- tests/dsp/test_interpolator.py | 6 +++--- tests/dsp/test_to_complex_bb.py | 2 +- tests/dsp/test_to_real_pb.py | 2 +- tests/simulation/test_awgn.py | 4 ++-- 9 files changed, 13 insertions(+), 13 deletions(-) diff --git a/docs/examples/phase-locked-loop.ipynb b/docs/examples/phase-locked-loop.ipynb index 3a8fabfa0..d3c6a4131 100644 --- a/docs/examples/phase-locked-loop.ipynb +++ b/docs/examples/phase-locked-loop.ipynb @@ -208,7 +208,7 @@ "N = 75 # samples\n", "theta_0 = 2 * np.pi / 10 # radians/sample\n", "x = np.exp(1j * (theta_0 * np.arange(N) + np.pi)) # Input signal\n", - "y = np.ones(x.size + 1, dtype=np.complex64) # Output signal\n", + "y = np.ones(x.size + 1, dtype=complex) # Output signal\n", "phase_error = np.zeros(x.size)\n", "freq_estimate = np.zeros(x.size)\n", "\n", diff --git a/src/sdr/plot/_time_domain.py b/src/sdr/plot/_time_domain.py index d5fcb5876..f457a150b 100644 --- a/src/sdr/plot/_time_domain.py +++ b/src/sdr/plot/_time_domain.py @@ -60,7 +60,7 @@ def time_domain( .. ipython:: python # Create a QPSK impulse signal - x = np.zeros(1000, dtype=np.complex64); \ + x = np.zeros(1000, dtype=complex); \ symbol_map = np.exp(1j * np.pi / 4) * np.array([1, 1j, -1, -1j]); \ x[::10] = symbol_map[np.random.randint(0, 4, 100)] diff --git a/tests/dsp/test_decimator.py b/tests/dsp/test_decimator.py index 1c0636882..ac6be10a5 100644 --- a/tests/dsp/test_decimator.py +++ b/tests/dsp/test_decimator.py @@ -29,7 +29,7 @@ # fir = sdr.Decimator(r) # y = fir(x, mode) -# xr = np.zeros(N * r, dtype=np.complex64) +# xr = np.zeros(N * r, dtype=complex) # xr[::r] = x[:] # y_truth = scipy.signal.convolve(xr, fir.taps, mode=mode) diff --git a/tests/dsp/test_fir.py b/tests/dsp/test_fir.py index 3a87539da..efa345dc4 100644 --- a/tests/dsp/test_fir.py +++ b/tests/dsp/test_fir.py @@ -27,7 +27,7 @@ def test_streaming(): fir = sdr.FIR(h, streaming=True) d = 5 # Stride - y = np.zeros(N, dtype=np.complex64) + y = np.zeros(N, dtype=complex) for i in range(0, N, d): y[i : i + d] = fir(x[i : i + d]) diff --git a/tests/dsp/test_iir.py b/tests/dsp/test_iir.py index c35550155..bac1ec767 100644 --- a/tests/dsp/test_iir.py +++ b/tests/dsp/test_iir.py @@ -16,8 +16,8 @@ def test_single_pole_impulse_response(pole): x = np.zeros(N, dtype=float) x[0] = 1 - b = np.array([1], dtype=np.complex64) - a = np.array([1, -pole], dtype=np.complex64) + b = np.array([1], dtype=complex) + a = np.array([1, -pole], dtype=complex) iir = sdr.IIR(b, a) h = iir(x) h_truth = pole ** np.arange(N) diff --git a/tests/dsp/test_interpolator.py b/tests/dsp/test_interpolator.py index f9b46fed1..e957deeb1 100644 --- a/tests/dsp/test_interpolator.py +++ b/tests/dsp/test_interpolator.py @@ -30,7 +30,7 @@ def test_non_streaming_full(): fir = sdr.Interpolator(r) y = fir(x, mode) - xr = np.zeros(N * r, dtype=np.complex64) + xr = np.zeros(N * r, dtype=complex) xr[::r] = x[:] y_truth = scipy.signal.convolve(xr, fir.taps, mode=mode) @@ -46,11 +46,11 @@ def test_streaming(): fir = sdr.Interpolator(r, streaming=True) d = 10 # Stride - y = np.zeros(N * r, dtype=np.complex64) + y = np.zeros(N * r, dtype=complex) for i in range(0, N, d): y[i * r : (i + d) * r] = fir(x[i : i + d]) - xr = np.zeros(N * r, dtype=np.complex64) + xr = np.zeros(N * r, dtype=complex) xr[::r] = x[:] y_truth = scipy.signal.convolve(xr, fir.taps, mode="full")[0 : N * r] diff --git a/tests/dsp/test_to_complex_bb.py b/tests/dsp/test_to_complex_bb.py index 20f484b2f..c3da6dafa 100644 --- a/tests/dsp/test_to_complex_bb.py +++ b/tests/dsp/test_to_complex_bb.py @@ -11,5 +11,5 @@ def test_exceptions(): sdr.to_complex_bb(x_r) with pytest.raises(ValueError): # x_r must be real - x_r = np.zeros(4, dtype=np.complex64) + x_r = np.zeros(4, dtype=complex) sdr.to_complex_bb(x_r) diff --git a/tests/dsp/test_to_real_pb.py b/tests/dsp/test_to_real_pb.py index 4bb04b7c9..439f669b3 100644 --- a/tests/dsp/test_to_real_pb.py +++ b/tests/dsp/test_to_real_pb.py @@ -7,7 +7,7 @@ def test_exceptions(): with pytest.raises(ValueError): # x_c must be 1D - x_c = np.zeros((2, 2), dtype=np.complex64) + x_c = np.zeros((2, 2), dtype=complex) sdr.to_real_pb(x_c) with pytest.raises(ValueError): # x_c must be real diff --git a/tests/simulation/test_awgn.py b/tests/simulation/test_awgn.py index bd02dd3b0..89c78b73e 100644 --- a/tests/simulation/test_awgn.py +++ b/tests/simulation/test_awgn.py @@ -30,7 +30,7 @@ def test_real_from_snr(): def test_complex_from_snr(): A = np.random.uniform(2, 3) # Signal amplitude N = 1000 - x = A * np.ones(N, dtype=np.complex64) # Use a constant signal so the power measurement has no uncertainty + x = A * np.ones(N, dtype=complex) # Use a constant signal so the power measurement has no uncertainty snr = np.random.uniform(1, 2) # dB snr_linear = 10 ** (snr / 10) @@ -60,7 +60,7 @@ def test_real_from_noise(): def test_complex_from_noise(): A = np.random.uniform(2, 3) # Signal amplitude N = 1000 - x = A * np.ones(N, dtype=np.complex64) # Use a constant signal so the power measurement has no uncertainty + x = A * np.ones(N, dtype=complex) # Use a constant signal so the power measurement has no uncertainty noise_std = np.random.uniform(1, 2) y = sdr.awgn(x, noise=noise_std**2) From 272624ce2fda0aa9c37577b2c62998d31beefa4c Mon Sep 17 00:00:00 2001 From: mhostetter Date: Wed, 16 Aug 2023 12:42:48 -0400 Subject: [PATCH 20/43] Fix bug in scaling PSDs --- src/sdr/plot/_spectral_estimation.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/sdr/plot/_spectral_estimation.py b/src/sdr/plot/_spectral_estimation.py index 7945efe8e..2754f2d91 100644 --- a/src/sdr/plot/_spectral_estimation.py +++ b/src/sdr/plot/_spectral_estimation.py @@ -63,7 +63,7 @@ def periodogram( return_onesided=x_axis != "two-sided", average=average, ) - Pxx = 10 * np.log10(np.abs(Pxx) ** 2) + Pxx = 10 * np.log10(Pxx) if x_axis == "two-sided": f[f >= 0.5 * sample_rate] -= sample_rate # Wrap frequencies from [0, 1) to [-0.5, 0.5) @@ -143,7 +143,7 @@ def spectrogram( return_onesided=x_axis != "two-sided", mode="psd", ) - Sxx = 10 * np.log10(np.abs(Sxx)) + Sxx = 10 * np.log10(Sxx) if x_axis == "one-sided" and np.iscomplexobj(x): # If complex data, the spectrogram always returns a two-sided spectrum. So we need to remove the second half. From 556ee879dbc6064ef9da44fc5120a48181252cb8 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Wed, 16 Aug 2023 12:51:50 -0400 Subject: [PATCH 21/43] Fix bug in dB to linear conversion --- src/sdr/_conversion.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/sdr/_conversion.py b/src/sdr/_conversion.py index f280f314b..61de42775 100644 --- a/src/sdr/_conversion.py +++ b/src/sdr/_conversion.py @@ -126,8 +126,6 @@ def linear( conversions-decibels """ x = np.asarray(x) - if not np.all(x >= 0): - raise ValueError("Argument 'x' must be non-negative.") if type in ["value", "power"]: return 10 ** (x / 10) From 0027aa49bad7726ae4aba17536d9f3a196b4d825 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Wed, 16 Aug 2023 13:44:16 -0400 Subject: [PATCH 22/43] Fix bug in FSPL in the near-field region --- src/sdr/_link_budget/_path_loss.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/sdr/_link_budget/_path_loss.py b/src/sdr/_link_budget/_path_loss.py index fa260fa71..14d908357 100644 --- a/src/sdr/_link_budget/_path_loss.py +++ b/src/sdr/_link_budget/_path_loss.py @@ -50,5 +50,14 @@ def fspl(d: npt.ArrayLike, f: npt.ArrayLike) -> np.ndarray: """ d = np.asarray(d) f = np.asarray(f) + + # The free-space path loss equation is only valid in the far field. For very small distances, the FSPL + # equation could have a negative result, implying a path gain. This is not possible. So for distances less + # than lambda / (4 * pi), the path loss is set to 0 dB. + lambda_ = scipy.constants.speed_of_light / f + d = np.maximum(d, lambda_ / (4 * np.pi)) + + # The free-space path loss equation loss = 20 * np.log10(4 * np.pi * d * f / scipy.constants.speed_of_light) + return loss From abeda571c7d36dc0e425703178c2b61f3f4378a6 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Wed, 16 Aug 2023 13:53:37 -0400 Subject: [PATCH 23/43] Add an example plot to `fspl()` --- src/sdr/_link_budget/_path_loss.py | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/src/sdr/_link_budget/_path_loss.py b/src/sdr/_link_budget/_path_loss.py index 14d908357..ffa38a002 100644 --- a/src/sdr/_link_budget/_path_loss.py +++ b/src/sdr/_link_budget/_path_loss.py @@ -24,6 +24,11 @@ def fspl(d: npt.ArrayLike, f: npt.ArrayLike) -> np.ndarray: Returns: The free-space path loss (FSPL) in dB. + Note: + The free-space path loss equation is only valid in the far field. For $d < \lambda / 4 \pi$, the FSPL + equation has a negative result, implying a path gain. This is not possible. So these path losses + are set to 0 dB. + Examples: Compute the free-space path loss for a 1 km link at 1 GHz. @@ -45,6 +50,21 @@ def fspl(d: npt.ArrayLike, f: npt.ArrayLike) -> np.ndarray: sdr.fspl(1e3, 2e9) + Plot the free-space path loss at 1 GHz for distances up to 1 km. + + .. ipython:: python + + d = np.linspace(0, 1_000, 1000) # m + fspl = sdr.fspl(d, 1e9) # dB + + @savefig sdr_fspl_1.png + plt.figure(figsize=(8, 4)); \ + plt.plot(d, fspl); \ + plt.xlabel('Distance (m)'); \ + plt.ylabel('Free-space path loss (dB)'); \ + plt.grid(True); \ + plt.tight_layout() + Group: link-budget-path-losses """ From d3771204e16fda7501a4198c1b3c46889c432aa2 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Wed, 16 Aug 2023 15:41:02 -0400 Subject: [PATCH 24/43] Test FSPL in near field --- tests/link_budgets/test_fspl.py | 41 +++++++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) diff --git a/tests/link_budgets/test_fspl.py b/tests/link_budgets/test_fspl.py index 8be39ba46..bfc4e0d3a 100644 --- a/tests/link_budgets/test_fspl.py +++ b/tests/link_budgets/test_fspl.py @@ -125,3 +125,44 @@ def test_28p25_GHz(): fspl = sdr.fspl(d[i], 28.25e9) assert isinstance(fspl, float) assert fspl == pytest.approx(fspl_truth[i]) + + +def test_300_MHz_near_field(): + """ + Matlab: + >> d = linspace(0, 0.25, 20); d' + >> L = fspl(d, physconst('LightSpeed')/300e6); L' + """ + d = np.linspace(0, 0.25, 20) + fspl = sdr.fspl(d, 300e6) + fspl_truth = np.array( + [ + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 1.275897270945932, + 2.435736210499666, + 3.458786659447293, + 4.373936470660794, + 5.201790173825296, + 5.957561391613292, + 6.652803516797531, + 7.296497184225555, + 7.895761651774420, + 8.456336123779291, + 8.982914898226273, + 9.479386572726916, + 9.949008489717375, + ] + ) + assert isinstance(fspl, np.ndarray) + assert np.allclose(fspl, fspl_truth) + + i = np.random.randint(0, d.size) + fspl = sdr.fspl(d[i], 300e6) + assert isinstance(fspl, float) + assert fspl == pytest.approx(fspl_truth[i]) From c91f5ef6c772298c426970bb6c2e16923a11733d Mon Sep 17 00:00:00 2001 From: mhostetter Date: Wed, 16 Aug 2023 15:54:19 -0400 Subject: [PATCH 25/43] Allow for returning power measurements in dB --- src/sdr/_measurement/_power.py | 25 ++++++++++++++++++------- 1 file changed, 18 insertions(+), 7 deletions(-) diff --git a/src/sdr/_measurement/_power.py b/src/sdr/_measurement/_power.py index 98fcc4221..0ac2de231 100644 --- a/src/sdr/_measurement/_power.py +++ b/src/sdr/_measurement/_power.py @@ -6,31 +6,37 @@ import numpy as np import numpy.typing as npt +from .._conversion import db as to_db from .._helper import export @export -def peak_power(x: npt.ArrayLike) -> float: +def peak_power(x: npt.ArrayLike, db: bool = False) -> float: r""" Measures the peak power of a time-domain signal $x[n]$. - $$P_{\text{peak}} = \max \left( \left| x[n] \right|^2 \right)$$ + $$P_{\text{peak}} = \max \left| x[n] \right|^2$$ Arguments: x: The time-domain signal $x[n]$ to measure. + db: Indicates whether to return the result in dB. Returns: - The peak power of $x[n]$ in units^2. + The peak power. If `db=False`, $P_{\text{peak}}$ is returned. + If `db=True`, $10 \log_{10} P_{\text{peak}}$ is returned. Group: measurement-power """ x = np.asarray(x) - return np.max(np.abs(x) ** 2) + P_peak = np.max(np.abs(x) ** 2) + if db: + P_peak = to_db(P_peak, type="power") + return P_peak @export -def average_power(x: npt.ArrayLike) -> float: +def average_power(x: npt.ArrayLike, db: bool = False) -> float: r""" Measures the average power of a time-domain signal $x[n]$. @@ -38,15 +44,20 @@ def average_power(x: npt.ArrayLike) -> float: Arguments: x: The time-domain signal $x[n]$ to measure. + db: Indicates whether to return the result in dB. Returns: - The average power of $x[n]$ in units^2. + The average power. If `db=False`, $P_{\text{avg}}$ is returned. + If `db=True`, $10 \log_{10} P_{\text{avg}}$ is returned. Group: measurement-power """ x = np.asarray(x) - return np.mean(np.abs(x) ** 2) + P_avg = np.mean(np.abs(x) ** 2) + if db: + P_avg = to_db(P_avg, type="power") + return P_avg @export From e096808ada97b006971452f625d303cf29a08465 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Wed, 16 Aug 2023 15:59:04 -0400 Subject: [PATCH 26/43] Allow for returning voltage measurements in dB --- src/sdr/_measurement/_voltage.py | 25 ++++++++++++++++++------- 1 file changed, 18 insertions(+), 7 deletions(-) diff --git a/src/sdr/_measurement/_voltage.py b/src/sdr/_measurement/_voltage.py index 710a63235..dcaeeca58 100644 --- a/src/sdr/_measurement/_voltage.py +++ b/src/sdr/_measurement/_voltage.py @@ -6,31 +6,37 @@ import numpy as np import numpy.typing as npt +from .._conversion import db as to_db from .._helper import export @export -def peak_voltage(x: npt.ArrayLike) -> float: +def peak_voltage(x: npt.ArrayLike, db: bool = False) -> float: r""" Measures the peak voltage of a time-domain signal $x[n]$. - $$V_{\text{peak}} = \max \left( \left| x[n] \right| \right)$$ + $$V_{\text{peak}} = \max \left| x[n] \right|$$ Arguments: x: The time-domain signal $x[n]$ to measure. + db: Indicates whether to return the result in dB. Returns: - The peak voltage of $x[n]$ in units. + The peak voltage. If `db=False`, $V_{\text{peak}}$ is returned. + If `db=True`, $20 \log_{10} V_{\text{peak}}$ is returned. Group: measurement-voltage """ x = np.asarray(x) - return np.max(np.abs(x)) + V_peak = np.max(np.abs(x)) + if db: + V_peak = to_db(V_peak, type="voltage") + return V_peak @export -def rms_voltage(x: npt.ArrayLike) -> float: +def rms_voltage(x: npt.ArrayLike, db: bool = False) -> float: r""" Measures the root-mean-square (RMS) voltage of a time-domain signal $x[n]$. @@ -38,15 +44,20 @@ def rms_voltage(x: npt.ArrayLike) -> float: Arguments: x: The time-domain signal $x[n]$ to measure. + db: Indicates whether to return the result in dB. Returns: - The RMS voltage of $x[n]$ in units. + The root-mean-square voltage. If `db=False`, $V_{\text{rms}}$ is returned. + If `db=True`, $20 \log_{10} V_{\text{rms}}$ is returned. Group: measurement-voltage """ x = np.asarray(x) - return np.sqrt(np.mean(np.abs(x) ** 2)) + V_rms = np.sqrt(np.mean(np.abs(x) ** 2)) + if db: + V_rms = to_db(V_rms, type="voltage") + return V_rms @export From a514d9ad459565f2ba5b0ea7c63c3b1af0a06a87 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Fri, 18 Aug 2023 17:08:52 -0400 Subject: [PATCH 27/43] Support real mixing --- src/sdr/_signal.py | 29 ++++++++++++++++++++++------- 1 file changed, 22 insertions(+), 7 deletions(-) diff --git a/src/sdr/_signal.py b/src/sdr/_signal.py index b8a21f316..cb5461cf0 100644 --- a/src/sdr/_signal.py +++ b/src/sdr/_signal.py @@ -11,18 +11,26 @@ @export -def mix(x: npt.ArrayLike, freq: float = 0, phase: float = 0, sample_rate: float = 1) -> np.ndarray: +def mix( + x: npt.ArrayLike, + freq: float = 0, + phase: float = 0, + sample_rate: float = 1, + complex: bool = True, # pylint: disable=redefined-builtin +) -> np.ndarray: r""" - Mixes the time-domain signal $x[n]$ with a complex exponential. - - $$y[n] = x[n] \cdot \exp \left[ j \left( \frac{2 \pi f}{f_s} n + \phi \right) \right]$$ + Mixes the time-domain signal $x[n]$ with a complex exponential or real sinusoid. Arguments: x: The time-domain signal $x[n]$. - freq: The frequency $f$ of the complex exponential in Hz (or 1/samples if `sample_rate=1`). + freq: The frequency $f$ of the sinusoid in Hz (or 1/samples if `sample_rate=1`). The frequency must satisfy $-f_s/2 \le f \le f_s/2$. - phase: The phase $\phi$ of the complex exponential in degrees. + phase: The phase $\phi$ of the sinusoid in degrees. sample_rate: The sample rate $f_s$ of the signal. + complex: Indicates whether to mix by a complex exponential or real sinusoid. + + - `True`: $y[n] = x[n] \cdot \exp \left[ j \left( \frac{2 \pi f}{f_s} n + \phi \right) \right]$ + - `False`: $y[n] = x[n] \cdot \cos \left( \frac{2 \pi f}{f_s} n + \phi \right)$ Returns: The mixed signal $y[n]$. @@ -65,6 +73,8 @@ def mix(x: npt.ArrayLike, freq: float = 0, phase: float = 0, sample_rate: float raise TypeError(f"Argument 'phase' must be a number, not {type(phase)}.") if not isinstance(sample_rate, (int, float)): raise TypeError(f"Argument 'sample_rate' must be a number, not {type(sample_rate)}.") + if not isinstance(complex, bool): + raise TypeError(f"Argument 'complex' must be a bool, not {type(complex)}.") if not -sample_rate / 2 <= freq <= sample_rate / 2: raise ValueError(f"Argument 'freq' must be in the range [{-sample_rate/2}, {sample_rate/2}], not {freq}.") @@ -72,7 +82,12 @@ def mix(x: npt.ArrayLike, freq: float = 0, phase: float = 0, sample_rate: float raise ValueError(f"Argument 'sample_rate' must be positive, not {sample_rate}.") t = np.arange(len(x)) / sample_rate # Time vector in seconds - y = x * np.exp(1j * (2 * np.pi * freq * t + np.deg2rad(phase))) + if complex: + lo = np.exp(1j * (2 * np.pi * freq * t + np.deg2rad(phase))) + else: + lo = np.cos(2 * np.pi * freq * t + np.deg2rad(phase)) + + y = x * lo return y From e5cff0326f179c325a3d959952e73e94789c70c1 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Fri, 18 Aug 2023 18:57:44 -0400 Subject: [PATCH 28/43] Add `sdr.upsample()` Implements #99 --- src/sdr/_signal.py | 54 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 54 insertions(+) diff --git a/src/sdr/_signal.py b/src/sdr/_signal.py index cb5461cf0..8fdf24bc1 100644 --- a/src/sdr/_signal.py +++ b/src/sdr/_signal.py @@ -258,3 +258,57 @@ def to_real_pb(x_c: npt.ArrayLike) -> np.ndarray: x_r = x_c.real return x_r + + +@export +def upsample(x: npt.ArrayLike, rate: int) -> np.ndarray: + r""" + Upsamples the time-domain signal $x[n]$ by the factor $r$. + + Warning: + This function does not perform any anti-alias filtering. The upsampled signal $y[n]$ will have + frequency components above the Nyquist frequency of the original signal $x[n]$. + + Arguments: + x: The time-domain signal $x[n]$ with sample rate $f_s$. + rate: The upsampling factor $r$. + + Returns: + The upsampled signal $y[n]$ with sample rate $f_s r$. + + Examples: + Upsample a complex exponential by a factor of 4. + + .. ipython:: python + + x = np.exp(1j*2*np.pi/16*np.arange(20)); x + y = sdr.upsample(x, 4); y + + @savefig sdr_upsample_1.png + plt.figure(figsize=(8, 4)); \ + sdr.plot.time_domain(x, sample_rate=1, label="$x[n]$"); \ + sdr.plot.time_domain(y, sample_rate=4, label="$y[n]$"); + + The spectrum of $y[n]$ has 3 additional copies of the spectrum of $x[n]$. + + .. ipython:: python + + @savefig sdr_upsample_2.png + plt.figure(figsize=(8, 4)); \ + sdr.plot.periodogram(x, fft=2048, sample_rate=1, label="$x[n]$"); \ + sdr.plot.periodogram(y, fft=2048, sample_rate=4, label="$y[n]$"); + + Group: + dsp-signal-manipulation + """ + x = np.asarray(x) + + if not isinstance(rate, int): + raise TypeError(f"Argument 'rate' must be an int, not {type(rate)}.") + if not rate > 0: + raise ValueError(f"Argument 'rate' must be positive, not {rate}.") + + y = np.zeros(x.size * rate, dtype=x.dtype) + y[::rate] = x + + return y From 956c57254c62a5d0b5ed7c3412f8465120759196 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Fri, 18 Aug 2023 18:57:58 -0400 Subject: [PATCH 29/43] Add unit tests for `sdr.upsample()` --- tests/dsp/test_upsample.py | 70 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 70 insertions(+) create mode 100644 tests/dsp/test_upsample.py diff --git a/tests/dsp/test_upsample.py b/tests/dsp/test_upsample.py new file mode 100644 index 000000000..13964ccc3 --- /dev/null +++ b/tests/dsp/test_upsample.py @@ -0,0 +1,70 @@ +import numpy as np +import pytest + +import sdr + + +def test_exceptions(): + x = np.random.randn(10) + + with pytest.raises(TypeError): + # Rate must be an integer + sdr.upsample(x, 4.0) + with pytest.raises(ValueError): + # Rate must be positive + sdr.upsample(x, 0) + + +def test_1(): + """ + Matlab: + >> x = 0:9; + >> upsample(x, 4)' + """ + x = np.arange(10) + y = sdr.upsample(x, 4) + y_truth = np.array( + [ + 0, + 0, + 0, + 0, + 1, + 0, + 0, + 0, + 2, + 0, + 0, + 0, + 3, + 0, + 0, + 0, + 4, + 0, + 0, + 0, + 5, + 0, + 0, + 0, + 6, + 0, + 0, + 0, + 7, + 0, + 0, + 0, + 8, + 0, + 0, + 0, + 9, + 0, + 0, + 0, + ] + ) + assert np.array_equal(y, y_truth) From 1173ed2184e8e3ccf72285fd669d073d4f1467f9 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Fri, 18 Aug 2023 19:02:35 -0400 Subject: [PATCH 30/43] Improve warning statement --- src/sdr/_signal.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/sdr/_signal.py b/src/sdr/_signal.py index 8fdf24bc1..b1f651d50 100644 --- a/src/sdr/_signal.py +++ b/src/sdr/_signal.py @@ -266,8 +266,9 @@ def upsample(x: npt.ArrayLike, rate: int) -> np.ndarray: Upsamples the time-domain signal $x[n]$ by the factor $r$. Warning: - This function does not perform any anti-alias filtering. The upsampled signal $y[n]$ will have - frequency components above the Nyquist frequency of the original signal $x[n]$. + This function does not perform any anti-aliasing filtering. The upsampled signal $y[n]$ will have + frequency components above the Nyquist frequency of the original signal $x[n]$. For efficient + polyphase interpolation (with anti-aliasing filtering), see :class:`sdr.Interpolator`. Arguments: x: The time-domain signal $x[n]$ with sample rate $f_s$. @@ -276,13 +277,16 @@ def upsample(x: npt.ArrayLike, rate: int) -> np.ndarray: Returns: The upsampled signal $y[n]$ with sample rate $f_s r$. + See Also: + sdr.Interpolator + Examples: Upsample a complex exponential by a factor of 4. .. ipython:: python - x = np.exp(1j*2*np.pi/16*np.arange(20)); x - y = sdr.upsample(x, 4); y + x = np.exp(1j * 2 * np.pi / 16 * np.arange(20)) + y = sdr.upsample(x, 4) @savefig sdr_upsample_1.png plt.figure(figsize=(8, 4)); \ From db66961cee5aae00f8f0c53bd857c2f5b5772ee6 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Fri, 18 Aug 2023 20:40:20 -0400 Subject: [PATCH 31/43] Add `sdr.downsample()` Implements #100 --- src/sdr/_signal.py | 76 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 76 insertions(+) diff --git a/src/sdr/_signal.py b/src/sdr/_signal.py index b1f651d50..4829e5c7c 100644 --- a/src/sdr/_signal.py +++ b/src/sdr/_signal.py @@ -316,3 +316,79 @@ def upsample(x: npt.ArrayLike, rate: int) -> np.ndarray: y[::rate] = x return y + + +@export +def downsample(x: npt.ArrayLike, rate: int) -> np.ndarray: + r""" + Downsamples the time-domain signal $x[n]$ by the factor $r$. + + Warning: + This function does not perform any anti-aliasing filtering. The downsampled signal $y[n]$ will have + spectral aliasing. For efficient polyphase decimation (with anti-aliasing filtering), see + :class:`sdr.Decimator`. + + Arguments: + x: The time-domain signal $x[n]$ with sample rate $f_s$. + rate: The downsampling factor $r$. + + Returns: + The downsampled signal $y[n]$ with sample rate $f_s / r$. + + Examples: + Downsample a complex exponential by a factor of 4. + + .. ipython:: python + + sample_rate = 400; \ + x1 = np.exp(1j * 2 * np.pi * 0 / sample_rate * np.arange(200)); \ + x2 = np.exp(1j * 2 * np.pi * 130 / sample_rate * np.arange(200)); \ + x3 = np.exp(1j * 2 * np.pi * -140 / sample_rate * np.arange(200)); \ + x = x1 + x2 + x3 + y = sdr.downsample(x, 4) + + @savefig sdr_downsample_1.png + plt.figure(figsize=(8, 4)); \ + sdr.plot.time_domain(x, sample_rate=sample_rate); \ + plt.title("Input signal $x[n]$"); \ + plt.tight_layout(); + + @savefig sdr_downsample_2.png + plt.figure(figsize=(8, 4)); \ + sdr.plot.time_domain(y, sample_rate=sample_rate/4); \ + plt.title("Downsampled signal $y[n]$"); \ + plt.tight_layout(); + + The spectrum of $x[n]$ has aliased. Any spectral content above the Nyquist frequency of $f_s / 2$ + will *fold* into the spectrum of $y[n]$. The CW at 0 Hz remains at 0 Hz (unaliased). + The CW at 130 Hz folds into 30 Hz. The CW at -140 Hz folds into -40 Hz. + + .. ipython:: python + + @savefig sdr_downsample_3.png + plt.figure(figsize=(8, 4)); \ + sdr.plot.periodogram(x, fft=2048, sample_rate=sample_rate); \ + plt.xlim(-sample_rate/2, sample_rate/2); \ + plt.title("Input signal $x[n]$"); \ + plt.tight_layout(); + + @savefig sdr_downsample_4.png + plt.figure(figsize=(8, 4)); \ + sdr.plot.periodogram(y, fft=2048, sample_rate=sample_rate/4); \ + plt.xlim(-sample_rate/2, sample_rate/2); \ + plt.title("Downsampled signal $y[n]$"); \ + plt.tight_layout(); + + Group: + dsp-signal-manipulation + """ + x = np.asarray(x) + + if not isinstance(rate, int): + raise TypeError(f"Argument 'rate' must be an int, not {type(rate)}.") + if not rate > 0: + raise ValueError(f"Argument 'rate' must be positive, not {rate}.") + + y = x[::rate] + + return y From 808e2af2e4c75b02c6d31c8ec70e17d9367a6fe0 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Fri, 18 Aug 2023 20:41:27 -0400 Subject: [PATCH 32/43] Add unit tests for `sdr.downsample()` --- tests/dsp/test_downsample.py | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) create mode 100644 tests/dsp/test_downsample.py diff --git a/tests/dsp/test_downsample.py b/tests/dsp/test_downsample.py new file mode 100644 index 000000000..dbeb0f265 --- /dev/null +++ b/tests/dsp/test_downsample.py @@ -0,0 +1,27 @@ +import numpy as np +import pytest + +import sdr + + +def test_exceptions(): + x = np.random.randn(40) + + with pytest.raises(TypeError): + # Rate must be an integer + sdr.downsample(x, 4.0) + with pytest.raises(ValueError): + # Rate must be positive + sdr.downsample(x, 0) + + +def test_1(): + """ + Matlab: + >> x = 0:39; + >> downsample(x, 4)' + """ + x = np.arange(40) + y = sdr.downsample(x, 4) + y_truth = np.array([0, 4, 8, 12, 16, 20, 24, 28, 32, 36]) + assert np.array_equal(y, y_truth) From 565cb1f1ce85042344065541b24f8d0beeab56cd Mon Sep 17 00:00:00 2001 From: mhostetter Date: Fri, 18 Aug 2023 20:58:06 -0400 Subject: [PATCH 33/43] Improve example plots --- src/sdr/_signal.py | 29 +++++++++++++++++++++++------ 1 file changed, 23 insertions(+), 6 deletions(-) diff --git a/src/sdr/_signal.py b/src/sdr/_signal.py index 4829e5c7c..0be006ec0 100644 --- a/src/sdr/_signal.py +++ b/src/sdr/_signal.py @@ -285,22 +285,39 @@ def upsample(x: npt.ArrayLike, rate: int) -> np.ndarray: .. ipython:: python - x = np.exp(1j * 2 * np.pi / 16 * np.arange(20)) + sample_rate = 100; \ + x = np.exp(1j * 2 * np.pi * 15 / sample_rate * np.arange(20)) y = sdr.upsample(x, 4) @savefig sdr_upsample_1.png plt.figure(figsize=(8, 4)); \ - sdr.plot.time_domain(x, sample_rate=1, label="$x[n]$"); \ - sdr.plot.time_domain(y, sample_rate=4, label="$y[n]$"); + sdr.plot.time_domain(x, sample_rate=sample_rate); \ + plt.title("Input signal $x[n]$"); \ + plt.tight_layout(); + + @savefig sdr_upsample_2.png + plt.figure(figsize=(8, 4)); \ + sdr.plot.time_domain(y, sample_rate=sample_rate*4); \ + plt.title("Upsampled signal $y[n]$"); \ + plt.tight_layout(); The spectrum of $y[n]$ has 3 additional copies of the spectrum of $x[n]$. .. ipython:: python - @savefig sdr_upsample_2.png + @savefig sdr_upsample_3.png plt.figure(figsize=(8, 4)); \ - sdr.plot.periodogram(x, fft=2048, sample_rate=1, label="$x[n]$"); \ - sdr.plot.periodogram(y, fft=2048, sample_rate=4, label="$y[n]$"); + sdr.plot.periodogram(x, fft=2048, sample_rate=sample_rate); \ + plt.xlim(-sample_rate*2, sample_rate*2); \ + plt.title("Input signal $x[n]$"); \ + plt.tight_layout(); + + @savefig sdr_upsample_4.png + plt.figure(figsize=(8, 4)); \ + sdr.plot.periodogram(y, fft=2048, sample_rate=sample_rate*4); \ + plt.xlim(-sample_rate*2, sample_rate*2); \ + plt.title("Upsampled signal $y[n]$"); \ + plt.tight_layout(); Group: dsp-signal-manipulation From 2a6c6df5fec85ae7e0f0842b28d5e9cdf8ed4d06 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Fri, 18 Aug 2023 21:15:06 -0400 Subject: [PATCH 34/43] Remove unnecessary conversion to scalar --- src/sdr/_link_budget/_capacity.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/sdr/_link_budget/_capacity.py b/src/sdr/_link_budget/_capacity.py index a9041e292..fb0ccf018 100644 --- a/src/sdr/_link_budget/_capacity.py +++ b/src/sdr/_link_budget/_capacity.py @@ -68,7 +68,7 @@ def bsc_capacity(p: npt.ArrayLike) -> np.ndarray: C = 1 - Hb(p) - return C if C.ndim > 0 else C.item() + return C @export @@ -116,7 +116,7 @@ def bec_capacity(p: npt.ArrayLike) -> np.ndarray: C = 1 - p - return C if C.ndim > 0 else C.item() + return C @export @@ -195,4 +195,4 @@ def awgn_capacity(snr: npt.ArrayLike, bandwidth: float | None = None) -> np.ndar else: C = np.log2(1 + snr_linear) # bits/2D - return C if C.ndim > 0 else C.item() + return C From f418e5088b0293668b78b24774c44a995d908463 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Fri, 18 Aug 2023 21:21:01 -0400 Subject: [PATCH 35/43] Suppress `log(0)` warning in entropy function Fixes #76 --- src/sdr/_link_budget/_capacity.py | 22 ++++++---------------- 1 file changed, 6 insertions(+), 16 deletions(-) diff --git a/src/sdr/_link_budget/_capacity.py b/src/sdr/_link_budget/_capacity.py index fb0ccf018..864545c43 100644 --- a/src/sdr/_link_budget/_capacity.py +++ b/src/sdr/_link_budget/_capacity.py @@ -5,6 +5,7 @@ import numpy as np import numpy.typing as npt +import scipy.stats from .._helper import export @@ -14,13 +15,8 @@ def Hb(x: npt.ArrayLike) -> np.ndarray: Computes the binary entropy function $H_b(x)$. """ x = np.asarray(x) - H = -x * np.log2(x) - (1 - x) * np.log2(1 - x) - # When x is 0 or 1, the binary entropy function is undefined. However, the limit of the binary entropy function - # as x approaches 0 or 1 is 0. Therefore, we can replace the undefined values with 0. - H = np.nan_to_num(H, nan=0) - - return H + return scipy.stats.entropy([x, 1 - x], base=2) @export @@ -66,9 +62,7 @@ def bsc_capacity(p: npt.ArrayLike) -> np.ndarray: if not (np.all(0 <= p) and np.all(p <= 1)): raise ValueError(f"Argument 'p' must be between 0 and 1, not {p}.") - C = 1 - Hb(p) - - return C + return 1 - Hb(p) @export @@ -114,9 +108,7 @@ def bec_capacity(p: npt.ArrayLike) -> np.ndarray: if not (np.all(0 <= p) and np.all(p <= 1)): raise ValueError(f"Argument 'p' must be between 0 and 1, not {p}.") - C = 1 - p - - return C + return 1 - p @export @@ -191,8 +183,6 @@ def awgn_capacity(snr: npt.ArrayLike, bandwidth: float | None = None) -> np.ndar snr_linear = 10 ** (snr / 10) if bandwidth: - C = bandwidth * np.log2(1 + snr_linear) # bits/s - else: - C = np.log2(1 + snr_linear) # bits/2D + return bandwidth * np.log2(1 + snr_linear) # bits/s - return C + return np.log2(1 + snr_linear) # bits/2D From e435fba8454fd9548b4f481e57b3a20539337341 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Fri, 18 Aug 2023 21:27:49 -0400 Subject: [PATCH 36/43] Format `.toml` file --- pyproject.toml | 24 +++++++----------------- 1 file changed, 7 insertions(+), 17 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 1d38b45e1..669ec1dbf 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,18 +1,12 @@ [build-system] -requires = [ - "setuptools >= 62", - "wheel", - "setuptools_scm[toml] >= 6.2" -] +requires = ["setuptools >= 62", "wheel", "setuptools_scm[toml] >= 6.2"] [project] name = "sdr" -authors = [ - {name = "Matt Hostetter", email = "matthostetter@gmail.com"}, -] +authors = [{ name = "Matt Hostetter", email = "matthostetter@gmail.com" }] description = "A Python package for software-defined radio" readme = "README.md" -license = {text = "MIT"} +license = { text = "MIT" } keywords = [ "software-defined radio", "sdr", @@ -54,7 +48,7 @@ dependencies = [ "scipy", "matplotlib", # "numba >= 0.55, < 0.58", # v0.55 is needed for support of NumPy 1.21 - "typing_extensions >= 4.0.0", # v4.0.0 is needed for use of Self (Python 3.11+) and Literal (Python 3.8+) + "typing_extensions >= 4.0.0", # v4.0.0 is needed for use of Self (Python 3.11+) and Literal (Python 3.8+) ] dynamic = ["version"] @@ -93,11 +87,9 @@ where = ["src"] universal = false [tool.pylint] -ignore-paths = [ - "src/sdr/_version.py", -] +ignore-paths = ["src/sdr/_version.py"] disable = [ - "comparison-with-callable", # pylint doesn't understand metaclass properties + "comparison-with-callable", # pylint doesn't understand metaclass properties "fixme", "global-statement", "invalid-name", @@ -132,9 +124,7 @@ profile = "black" [tool.pytest.ini_options] minversion = "6.2" addopts = "-s --showlocals" -testpaths = [ - "tests" -] +testpaths = ["tests"] [tool.coverage.report] exclude_lines = [ From 18f2e16b40519d8c768fd3f1160ae78ede6a085e Mon Sep 17 00:00:00 2001 From: mhostetter Date: Fri, 18 Aug 2023 21:28:22 -0400 Subject: [PATCH 37/43] Exclude plot functions from coverage Once the library is more mature, unit tests against images can be added. It will be too cumbersome at this stage. --- pyproject.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/pyproject.toml b/pyproject.toml index 669ec1dbf..70140761b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -135,3 +135,4 @@ exclude_lines = [ "raise NotImplementedError", "raise RuntimeError", ] +omit = "src/sdr/plot/*" From ab2f229a26b3f7b4ac33df76e11495bce57ea8ed Mon Sep 17 00:00:00 2001 From: mhostetter Date: Fri, 18 Aug 2023 21:28:54 -0400 Subject: [PATCH 38/43] Expand dictionary --- .vscode/settings.json | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/.vscode/settings.json b/.vscode/settings.json index 2db195189..88ff65fdc 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -20,10 +20,13 @@ ], "cSpell.words": [ "awgn", - "Baseband", + "baseband", + "downsample", + "downsampled", + "downsamples", "downsampling", "exponentials", - "Feedforward", + "feedforward", "figsize", "hexdump", "intersymbol", @@ -46,6 +49,7 @@ "stopband", "underdamped", "upsampled", + "upsamples", "upsampling" ] } \ No newline at end of file From 8ab16cfd747806952fc56ab8513d9644408e5373 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Fri, 18 Aug 2023 21:38:13 -0400 Subject: [PATCH 39/43] Fix `.toml` syntax --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 70140761b..9d241c648 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -135,4 +135,4 @@ exclude_lines = [ "raise NotImplementedError", "raise RuntimeError", ] -omit = "src/sdr/plot/*" +omit = ["src/sdr/plot/*"] From 574d3508ccd0f3134a6caac12ecdd002c8e1dd65 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Fri, 18 Aug 2023 22:01:21 -0400 Subject: [PATCH 40/43] Fix codecov path --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 9d241c648..659eff017 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -135,4 +135,4 @@ exclude_lines = [ "raise NotImplementedError", "raise RuntimeError", ] -omit = ["src/sdr/plot/*"] +omit = ["*/plot/*"] From 0a7db7d4dc91e798e5a6181ff71e882bb5cae886 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Fri, 18 Aug 2023 22:18:15 -0400 Subject: [PATCH 41/43] Use `linear()` inside library --- docs/examples/pulse-shapes.ipynb | 16 ++++++++-------- src/sdr/_conversion.py | 14 +++++++------- src/sdr/_link_budget/_antenna.py | 3 ++- src/sdr/_link_budget/_capacity.py | 3 ++- src/sdr/_measurement/_power.py | 2 +- src/sdr/_modulation/_psk.py | 10 +++++----- src/sdr/_pll.py | 3 ++- src/sdr/_simulation/_impairment.py | 7 ++++--- 8 files changed, 31 insertions(+), 27 deletions(-) diff --git a/docs/examples/pulse-shapes.ipynb b/docs/examples/pulse-shapes.ipynb index d369f4a91..1b4f58ef3 100644 --- a/docs/examples/pulse-shapes.ipynb +++ b/docs/examples/pulse-shapes.ipynb @@ -218,10 +218,10 @@ "w, H_rc_0p9 = scipy.signal.freqz(rc_0p9, 1, worN=1024, whole=False, fs=sps)\n", "\n", "# Compute the relative power in the main lobe of the pulses\n", - "P_rect = 10 * np.log10(np.cumsum(np.abs(H_rect) ** 2) / np.sum(np.abs(H_rect) ** 2))\n", - "P_rc_0p1 = 10 * np.log10(np.cumsum(np.abs(H_rc_0p1) ** 2) / np.sum(np.abs(H_rc_0p1) ** 2))\n", - "P_rc_0p5 = 10 * np.log10(np.cumsum(np.abs(H_rc_0p5) ** 2) / np.sum(np.abs(H_rc_0p5) ** 2))\n", - "P_rc_0p9 = 10 * np.log10(np.cumsum(np.abs(H_rc_0p9) ** 2) / np.sum(np.abs(H_rc_0p9) ** 2))\n", + "P_rect = sdr.db(np.cumsum(np.abs(H_rect) ** 2) / np.sum(np.abs(H_rect) ** 2))\n", + "P_rc_0p1 = sdr.db(np.cumsum(np.abs(H_rc_0p1) ** 2) / np.sum(np.abs(H_rc_0p1) ** 2))\n", + "P_rc_0p5 = sdr.db(np.cumsum(np.abs(H_rc_0p5) ** 2) / np.sum(np.abs(H_rc_0p5) ** 2))\n", + "P_rc_0p9 = sdr.db(np.cumsum(np.abs(H_rc_0p9) ** 2) / np.sum(np.abs(H_rc_0p9) ** 2))\n", "\n", "plt.figure(figsize=(10, 5))\n", "plt.plot(w, P_rect, color=\"k\", linestyle=\":\", label=\"Rectangular\")\n", @@ -406,10 +406,10 @@ "w, H_srrc_0p9 = scipy.signal.freqz(srrc_0p9, 1, worN=1024, whole=False, fs=sps)\n", "\n", "# Compute the relative power in the main lobe of the pulses\n", - "P_rect = 10 * np.log10(np.cumsum(np.abs(H_rect) ** 2) / np.sum(np.abs(H_rect) ** 2))\n", - "P_srrc_0p1 = 10 * np.log10(np.cumsum(np.abs(H_srrc_0p1) ** 2) / np.sum(np.abs(H_srrc_0p1) ** 2))\n", - "P_srrc_0p5 = 10 * np.log10(np.cumsum(np.abs(H_srrc_0p5) ** 2) / np.sum(np.abs(H_srrc_0p5) ** 2))\n", - "P_srrc_0p9 = 10 * np.log10(np.cumsum(np.abs(H_srrc_0p9) ** 2) / np.sum(np.abs(H_srrc_0p9) ** 2))\n", + "P_rect = sdr.db(np.cumsum(np.abs(H_rect) ** 2) / np.sum(np.abs(H_rect) ** 2))\n", + "P_srrc_0p1 = sdr.db(np.cumsum(np.abs(H_srrc_0p1) ** 2) / np.sum(np.abs(H_srrc_0p1) ** 2))\n", + "P_srrc_0p5 = sdr.db(np.cumsum(np.abs(H_srrc_0p5) ** 2) / np.sum(np.abs(H_srrc_0p5) ** 2))\n", + "P_srrc_0p9 = sdr.db(np.cumsum(np.abs(H_srrc_0p9) ** 2) / np.sum(np.abs(H_srrc_0p9) ** 2))\n", "\n", "plt.figure(figsize=(10, 5))\n", "plt.plot(w, P_rect, color=\"k\", linestyle=\":\", label=\"Rectangular\")\n", diff --git a/src/sdr/_conversion.py b/src/sdr/_conversion.py index 61de42775..2789045c3 100644 --- a/src/sdr/_conversion.py +++ b/src/sdr/_conversion.py @@ -176,8 +176,8 @@ def ebn0_to_esn0(ebn0: npt.ArrayLike, bps: int, rate: int = 1) -> np.ndarray: conversions-from-ebn0 """ ebn0 = np.asarray(ebn0) # Energy per information bit - ecn0 = ebn0 + 10 * np.log10(rate) # Energy per coded bit - esn0 = ecn0 + 10 * np.log10(bps) # Energy per symbol + ecn0 = ebn0 + db(rate) # Energy per coded bit + esn0 = ecn0 + db(bps) # Energy per symbol return esn0 @@ -217,7 +217,7 @@ def ebn0_to_snr(ebn0: npt.ArrayLike, bps: int, rate: int = 1, sps: int = 1) -> n conversions-from-ebn0 """ esn0 = ebn0_to_esn0(ebn0, bps, rate=rate) # SNR per symbol - snr = esn0 - 10 * np.log10(sps) # SNR per sample + snr = esn0 - db(sps) # SNR per sample return snr @@ -261,8 +261,8 @@ def esn0_to_ebn0(esn0: npt.ArrayLike, bps: int, rate: int = 1) -> np.ndarray: conversions-from-esn0 """ esn0 = np.asarray(esn0) - ecn0 = esn0 - 10 * np.log10(bps) # Energy per coded bit - ebn0 = ecn0 - 10 * np.log10(rate) # Energy per information bit + ecn0 = esn0 - db(bps) # Energy per coded bit + ebn0 = ecn0 - db(rate) # Energy per information bit return ebn0 @@ -300,7 +300,7 @@ def esn0_to_snr(esn0: npt.ArrayLike, sps: int = 1) -> np.ndarray: conversions-from-esn0 """ esn0 = np.asarray(esn0) # SNR per symbol - snr = esn0 - 10 * np.log10(sps) # SNR per sample + snr = esn0 - db(sps) # SNR per sample return snr @@ -384,5 +384,5 @@ def snr_to_esn0(snr: npt.ArrayLike, sps: int = 1) -> np.ndarray: conversions-from-snr """ snr = np.asarray(snr) - esn0 = snr + 10 * np.log10(sps) + esn0 = snr + db(sps) return esn0 diff --git a/src/sdr/_link_budget/_antenna.py b/src/sdr/_link_budget/_antenna.py index 537245088..e2235fb64 100644 --- a/src/sdr/_link_budget/_antenna.py +++ b/src/sdr/_link_budget/_antenna.py @@ -7,6 +7,7 @@ import numpy.typing as npt import scipy.constants +from .._conversion import db from .._helper import export @@ -63,7 +64,7 @@ def parabolic_antenna( lambda_ = scipy.constants.speed_of_light / freq # Wavelength in meters G = (np.pi * diameter / lambda_) ** 2 * efficiency # Gain in linear units - G = 10 * np.log10(G) # Gain in dBi + G = db(G) # Gain in dBi theta = np.arcsin(3.83 * lambda_ / (np.pi * diameter)) # Beamwidth in radians theta = np.rad2deg(theta) # Beamwidth in degrees diff --git a/src/sdr/_link_budget/_capacity.py b/src/sdr/_link_budget/_capacity.py index 864545c43..e2504f947 100644 --- a/src/sdr/_link_budget/_capacity.py +++ b/src/sdr/_link_budget/_capacity.py @@ -7,6 +7,7 @@ import numpy.typing as npt import scipy.stats +from .._conversion import linear from .._helper import export @@ -180,7 +181,7 @@ def awgn_capacity(snr: npt.ArrayLike, bandwidth: float | None = None) -> np.ndar link-budget-channel-capacity """ snr = np.asarray(snr) - snr_linear = 10 ** (snr / 10) + snr_linear = linear(snr) if bandwidth: return bandwidth * np.log2(1 + snr_linear) # bits/s diff --git a/src/sdr/_measurement/_power.py b/src/sdr/_measurement/_power.py index 0ac2de231..ad23712aa 100644 --- a/src/sdr/_measurement/_power.py +++ b/src/sdr/_measurement/_power.py @@ -83,4 +83,4 @@ def papr(x: npt.ArrayLike) -> float: measurement-power """ x = np.asarray(x) - return 10 * np.log10(peak_power(x) / average_power(x)) + return to_db(peak_power(x) / average_power(x)) diff --git a/src/sdr/_modulation/_psk.py b/src/sdr/_modulation/_psk.py index 7577d87a3..35a813efc 100644 --- a/src/sdr/_modulation/_psk.py +++ b/src/sdr/_modulation/_psk.py @@ -9,7 +9,7 @@ import scipy.special from typing_extensions import Literal -from .._conversion import ebn0_to_esn0, esn0_to_ebn0 +from .._conversion import ebn0_to_esn0, esn0_to_ebn0, linear from .._data import unpack from .._helper import export, extend_docstring from .._probability import Q @@ -130,9 +130,9 @@ def ber(self, ebn0: npt.ArrayLike | None = None, diff_encoded: bool = False) -> M = self.order k = self.bps ebn0 = np.asarray(ebn0) - ebn0_linear = 10 ** (ebn0 / 10) + ebn0_linear = linear(ebn0) esn0 = ebn0_to_esn0(ebn0, k) - esn0_linear = 10 ** (esn0 / 10) + esn0_linear = linear(esn0) if not diff_encoded: if M == 2: @@ -219,9 +219,9 @@ def ser(self, esn0: npt.ArrayLike | None = None, diff_encoded: bool = False) -> M = self.order k = self.bps esn0 = np.asarray(esn0) - esn0_linear = 10 ** (esn0 / 10) + esn0_linear = linear(esn0) ebn0 = esn0_to_ebn0(esn0, k) - ebn0_linear = 10 ** (ebn0 / 10) + ebn0_linear = linear(ebn0) if not diff_encoded: if M == 2: diff --git a/src/sdr/_pll.py b/src/sdr/_pll.py index cbdab4757..793210519 100644 --- a/src/sdr/_pll.py +++ b/src/sdr/_pll.py @@ -5,6 +5,7 @@ import numpy as np +from ._conversion import linear from ._filter import IIR from ._helper import export from ._loop_filter import LoopFilter @@ -189,7 +190,7 @@ def phase_error_variance(self, cn0: float) -> float: Examples: See the :ref:`phase-locked-loop` example. """ - cn0_linear = 10 ** (cn0 / 10) + cn0_linear = linear(cn0) return self.Bn / cn0_linear @property diff --git a/src/sdr/_simulation/_impairment.py b/src/sdr/_simulation/_impairment.py index 6c617cd35..8c97af423 100644 --- a/src/sdr/_simulation/_impairment.py +++ b/src/sdr/_simulation/_impairment.py @@ -6,6 +6,7 @@ import numpy as np import numpy.typing as npt +from .._conversion import linear from .._farrow import FarrowResampler from .._helper import export from .._measurement import average_power @@ -86,7 +87,7 @@ def awgn( """ x = np.asarray(x) if snr is not None: - snr_linear = 10 ** (snr / 10) + snr_linear = linear(snr) signal_power = average_power(x) noise_power = signal_power / snr_linear elif noise is not None: @@ -181,8 +182,8 @@ def iq_imbalance(x: npt.ArrayLike, amplitude: float, phase: float = 0) -> np.nda phase = np.deg2rad(phase) # TODO: Should the phase be negative for I? - gain_i = 10 ** (0.5 * amplitude / 20) * np.exp(1j * -0.5 * phase) - gain_q = 10 ** (-0.5 * amplitude / 20) * np.exp(1j * 0.5 * phase) + gain_i = linear(0.5 * amplitude, type="voltage") * np.exp(1j * -0.5 * phase) + gain_q = linear(-0.5 * amplitude, type="voltage") * np.exp(1j * 0.5 * phase) y = gain_i * x.real + 1j * gain_q * x.imag From c3b39c306326d84213b584534508a2b472cd8581 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Fri, 18 Aug 2023 22:28:53 -0400 Subject: [PATCH 42/43] Add `sdr.wavelength()` --- src/sdr/_link_budget/_antenna.py | 29 ++++++++++++++++++++++++++++- 1 file changed, 28 insertions(+), 1 deletion(-) diff --git a/src/sdr/_link_budget/_antenna.py b/src/sdr/_link_budget/_antenna.py index e2235fb64..27a802db1 100644 --- a/src/sdr/_link_budget/_antenna.py +++ b/src/sdr/_link_budget/_antenna.py @@ -11,6 +11,33 @@ from .._helper import export +@export +def wavelength(freq: float) -> float: + r""" + Calculates the wavelength $\lambda$ of a electromagnetic wave with frequency $f$. + + $$\lambda = \frac{c}{f}$$ + + Arguments: + freq: The frequency $f$ in Hz of the signal. + + Returns: + The wavelength $\lambda$ in meters. + + Examples: + The wavelength of a 1 GHz signal is 0.3 meters. + + .. ipython:: python + + sdr.wavelength(1e9) + + Group: + link-budget-antennas + """ + freq = np.asarray(freq) + return scipy.constants.speed_of_light / freq + + @export def parabolic_antenna( freq: float, @@ -62,7 +89,7 @@ def parabolic_antenna( if not np.all((0 <= efficiency) & (efficiency <= 1)): raise ValueError("Argument 'efficiency' must be between 0 and 1.") - lambda_ = scipy.constants.speed_of_light / freq # Wavelength in meters + lambda_ = wavelength(freq) # Wavelength in meters G = (np.pi * diameter / lambda_) ** 2 * efficiency # Gain in linear units G = db(G) # Gain in dBi From 3a794e1bf357246a66df57188ad70dafdcd45fd4 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Sun, 20 Aug 2023 20:29:08 -0400 Subject: [PATCH 43/43] Add release notes for v0.0.6 --- docs/release-notes/v0.0.md | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/docs/release-notes/v0.0.md b/docs/release-notes/v0.0.md index 8fa1d1c4e..1ff2896d5 100644 --- a/docs/release-notes/v0.0.md +++ b/docs/release-notes/v0.0.md @@ -4,6 +4,24 @@ tocdepth: 2 # v0.0 +## v0.0.6 + +*Released August 20, 2023* + +### Changes + +- Added raster plots in `sdr.plot.raster()`. +- Added eye diagrams in `sdr.plot.eye()`. +- Added upsampling (without anti-aliasing filtering) in `sdr.upsample()`. +- Added downsampling (without anti-aliasing filtering) in `sdr.downsample()`. +- Added wavelength calculation in `sdr.wavelength()`. +- Supported real sinusoid mixing. +- Supported returning measurements in dB. + +### Contributors + +- Matt Hostetter ([@mhostetter](https://github.com/mhostetter)) + ## v0.0.5 *Released August 13, 2023*