From 634ba1254285eba9bd9efcf47116a7eb5593f0d0 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Sun, 19 Nov 2023 11:21:20 -0500 Subject: [PATCH 1/2] Add `sdr.design_bandstop_fir()` Closes #149 --- src/sdr/_filter/_design_fir.py | 119 +++++++++++++++++++++++++++++++++ 1 file changed, 119 insertions(+) diff --git a/src/sdr/_filter/_design_fir.py b/src/sdr/_filter/_design_fir.py index 3d693177b..9fc9134f4 100644 --- a/src/sdr/_filter/_design_fir.py +++ b/src/sdr/_filter/_design_fir.py @@ -60,6 +60,16 @@ def _ideal_bandpass(order: int, center_freq: float, bandwidth: float) -> npt.NDA return h_ideal +def _ideal_bandstop(order: int, center_freq: float, bandwidth: float) -> npt.NDArray[np.float_]: + """ + Returns the ideal bandstop filter impulse response. + """ + h_lp = _ideal_lowpass(order, center_freq - bandwidth / 2) + h_hp = _ideal_highpass(order, center_freq + bandwidth / 2) + h_ideal = h_lp + h_hp + return h_ideal + + def _window( order: int, window: Literal["hamming", "hann", "blackman", "blackman-harris", "chebyshev", "kaiser"] @@ -391,3 +401,112 @@ def design_bandpass_fir( h = _normalize_passband(h, center_freq) return h + + +@export +def design_bandstop_fir( + order: int, + center_freq: float, + bandwidth: float, + window: None + | Literal["hamming", "hann", "blackman", "blackman-harris", "chebyshev", "kaiser"] + | npt.ArrayLike = "hamming", + atten: float = 60, +) -> npt.NDArray[np.float_]: + r""" + Designs a bandstop FIR filter impulse response $h[n]$ using the window method. + + Arguments: + order: The filter order $N$. Must be even. + center_freq: The center frequency $f_{center}$, normalized to the Nyquist frequency $f_s / 2$. + bandwidth: The two-sided bandwidth about $f_{center}$, normalized to the Nyquist frequency $f_s / 2$. + window: The time-domain window to use. + + - `None`: No windowing. Equivalently, a length-$N + 1$ vector of ones. + - `"hamming"`: Hamming window, see :func:`scipy.signal.windows.hamming`. + - `"hann"`: Hann window, see :func:`scipy.signal.windows.hann`. + - `"blackman"`: Blackman window, see :func:`scipy.signal.windows.blackman`. + - `"blackman-harris"`: Blackman-Harris window, see :func:`scipy.signal.windows.blackmanharris`. + - `"chebyshev"`: Chebyshev window, see :func:`scipy.signal.windows.chebwin`. + - `"kaiser"`: Kaiser window, see :func:`scipy.signal.windows.kaiser`. + - `npt.ArrayLike`: A custom window. Must be a length-$N + 1$ vector. + + atten: The sidelobe attenuation in dB. Only used if `window` is `"chebyshev"` or `"kaiser"`. + + Returns: + The filter impulse response $h[n]$ with length $N + 1$. The center of the larger passband has 0 dB gain. + + References: + - https://www.mathworks.com/help/dsp/ref/designbandstopfir.html + + Examples: + Design a length-101 bandstop FIR filter with center frequency $f_{center} = 0.4 \cdot f_s / 2$ + and bandwidth $0.75 \cdot f_s / 2$, using a Hamming window. + + .. ipython:: python + + h_hamming = sdr.design_bandstop_fir(100, 0.4, 0.75, window="hamming") + + @savefig sdr_design_bandstop_fir_1.png + plt.figure(figsize=(8, 4)); \ + sdr.plot.impulse_response(h_hamming); + + @savefig sdr_design_bandstop_fir_2.png + plt.figure(figsize=(8, 4)); \ + sdr.plot.magnitude_response(h_hamming); + + Compare filter designs using different windows. + + .. ipython:: python + + h_hann = sdr.design_bandstop_fir(100, 0.4, 0.75, window="hann"); \ + h_blackman = sdr.design_bandstop_fir(100, 0.4, 0.75, window="blackman"); \ + h_blackman_harris = sdr.design_bandstop_fir(100, 0.4, 0.75, window="blackman-harris"); \ + h_chebyshev = sdr.design_bandstop_fir(100, 0.4, 0.75, window="chebyshev"); \ + h_kaiser = sdr.design_bandstop_fir(100, 0.4, 0.75, window="kaiser") + + @savefig sdr_design_bandstop_fir_3.png + plt.figure(figsize=(8, 4)); \ + sdr.plot.magnitude_response(h_hamming, label="Hamming"); \ + sdr.plot.magnitude_response(h_hann, label="Hann"); \ + sdr.plot.magnitude_response(h_blackman, label="Blackman"); \ + sdr.plot.magnitude_response(h_blackman_harris, label="Blackman-Harris"); \ + sdr.plot.magnitude_response(h_chebyshev, label="Chebyshev"); \ + sdr.plot.magnitude_response(h_kaiser, label="Kaiser"); \ + plt.legend(); \ + plt.ylim(-100, 10); + + Group: + dsp-fir-filtering + """ + if not isinstance(order, int): + raise TypeError(f"Argument 'order' must be an integer, not {type(order).__name__}.") + if not order % 2 == 0: + raise ValueError(f"Argument 'order' must be even, not {order}.") + + if not isinstance(center_freq, (int, float)): + raise TypeError(f"Argument 'center_freq' must be a number, not {type(center_freq).__name__}.") + if not 0 <= center_freq <= 1: + raise ValueError(f"Argument 'center_freq' must be between 0 and 1, not {center_freq}.") + + if not isinstance(bandwidth, (int, float)): + raise TypeError(f"Argument 'bandwidth' must be a number, not {type(bandwidth).__name__}.") + if not 0 <= center_freq + bandwidth / 2 <= 1: + raise ValueError(f"Argument 'bandwidth' must be between 0 and 0.5 - 'center_freq', not {bandwidth}.") + if not 0 <= center_freq - bandwidth / 2 <= 1: + raise ValueError(f"Argument 'bandwidth' must be between 0 and 'center_freq', not {bandwidth}.") + + if not isinstance(atten, (int, float)): + raise TypeError(f"Argument 'atten' must be a number, not {type(atten).__name__}.") + if not 0 <= atten: + raise ValueError(f"Argument 'atten' must be non-negative, not {atten}.") + + h_ideal = _ideal_bandstop(order, center_freq, bandwidth) + h_window = _window(order, window, atten=atten) + h = h_ideal * h_window + if center_freq > 0.5: + h = _normalize_passband(h, 0) + else: + h = _normalize_passband(h, 1) + + return h From da0575f8887e10b9c376da4494629d4715810cb5 Mon Sep 17 00:00:00 2001 From: mhostetter Date: Sun, 19 Nov 2023 11:21:48 -0500 Subject: [PATCH 2/2] Add unit tests for `design_bandstop_fir()` --- tests/dsp/fir/test_design_bandstop_fir.py | 345 ++++++++++++++++++++++ 1 file changed, 345 insertions(+) create mode 100644 tests/dsp/fir/test_design_bandstop_fir.py diff --git a/tests/dsp/fir/test_design_bandstop_fir.py b/tests/dsp/fir/test_design_bandstop_fir.py new file mode 100644 index 000000000..f15a13466 --- /dev/null +++ b/tests/dsp/fir/test_design_bandstop_fir.py @@ -0,0 +1,345 @@ +import numpy as np +import pytest + +import sdr + +from .helper import verify_impulse_response + + +def test_exceptions(): + with pytest.raises(TypeError): + sdr.design_bandstop_fir("30", 0.4, 0.25) + with pytest.raises(ValueError): + sdr.design_bandstop_fir(29, 0.4, 0.25) + + with pytest.raises(TypeError): + sdr.design_bandstop_fir(30, "0.4", 0.25) + with pytest.raises(ValueError): + sdr.design_bandstop_fir(30, -0.1, 0.25) + + with pytest.raises(TypeError): + sdr.design_bandstop_fir(30, 0.4, "0.25") + with pytest.raises(ValueError): + sdr.design_bandstop_fir(30, 0.4, 0.61 * 2) + with pytest.raises(ValueError): + sdr.design_bandstop_fir(30, 0.4, 0.41 * 2) + + with pytest.raises(ValueError): + sdr.design_bandstop_fir(30, 0.4, 0.25, window="invalid") + with pytest.raises(ValueError): + sdr.design_bandstop_fir(30, 0.4, 0.25, window=np.ones(10)) + + +def test_custom(): + """ + MATLAB: + >> h = designBandstopFIR(FilterOrder=30, CenterFrequency=0.4, Bandwidth=0.25, Window="custom", CustomWindow=ones(1,31)); + >> transpose(h) + """ + h = sdr.design_bandstop_fir(30, 0.4, 0.25, window=None) + h_truth = np.array( + [ + 0.017963811098946, + 0.010989790375857, + -0.040483655096790, + -0.047470791892434, + 0.018274903941445, + 0.049789192200645, + 0.009251871522192, + 0.000000000000000, + 0.031142203848832, + -0.025642844210333, + -0.130105429160949, + -0.054396687090110, + 0.175429172086089, + 0.201401513132611, + -0.083266843699732, + 0.723490559333834, + -0.083266843699732, + 0.201401513132611, + 0.175429172086089, + -0.054396687090110, + -0.130105429160949, + -0.025642844210333, + 0.031142203848832, + 0.000000000000000, + 0.009251871522192, + 0.049789192200645, + 0.018274903941445, + -0.047470791892434, + -0.040483655096790, + 0.010989790375857, + 0.017963811098946, + ] + ) + verify_impulse_response(h, h_truth, atol=1e-1) + + +def test_hamming(): + """ + MATLAB: + >> h = designBandstopFIR(FilterOrder=30, CenterFrequency=0.4, Bandwidth=0.25, Window="hamming"); + >> transpose(h) + """ + h = sdr.design_bandstop_fir(30, 0.4, 0.25, window="hamming") + h_truth = np.array( + [ + 0.001295920763505, + 0.000892428133759, + -0.004372345234144, + -0.007185275937122, + 0.003826547884793, + 0.013918318028965, + 0.003319260363488, + 0.000000000000000, + 0.016514978737540, + -0.015773739691351, + -0.090339176322020, + -0.041586836521494, + 0.144296905913076, + 0.174392736680514, + -0.074331760946352, + 0.750655412505035, + -0.074331760946352, + 0.174392736680514, + 0.144296905913076, + -0.041586836521494, + -0.090339176322020, + -0.015773739691351, + 0.016514978737540, + 0.000000000000000, + 0.003319260363488, + 0.013918318028965, + 0.003826547884793, + -0.007185275937122, + -0.004372345234144, + 0.000892428133759, + 0.001295920763505, + ] + ) + verify_impulse_response(h, h_truth, atol=1e-3) + + +def test_hann(): + """ + MATLAB: + >> h = designBandstopFIR(FilterOrder=30, CenterFrequency=0.4, Bandwidth=0.25, Window="hann"); + >> transpose(h) + """ + h = sdr.design_bandstop_fir(30, 0.4, 0.25, window="hann") + h_truth = np.array( + [ + 0, + 0.000107362890179, + -0.001564707596353, + -0.004053095706492, + 0.002703194884773, + 0.011129374154929, + 0.002858002590150, + 0.000000000000000, + 0.015377718751261, + -0.015006421380158, + -0.087247368558493, + -0.040590876565245, + 0.141876386618995, + 0.172292816804829, + -0.073637062395384, + 0.752767467273434, + -0.073637062395384, + 0.172292816804829, + 0.141876386618995, + -0.040590876565245, + -0.087247368558493, + -0.015006421380158, + 0.015377718751261, + 0.000000000000000, + 0.002858002590150, + 0.011129374154929, + 0.002703194884773, + -0.004053095706492, + -0.001564707596353, + 0.000107362890179, + 0, + ] + ) + verify_impulse_response(h, h_truth, atol=1e-2) + + +def test_blackman(): + """ + MATLAB: + >> h = designBandstopFIR(FilterOrder=30, CenterFrequency=0.4, Bandwidth=0.25, Window="blackman"); + >> transpose(h) + """ + h = sdr.design_bandstop_fir(30, 0.4, 0.25, window="blackman") + h_truth = np.array( + [ + -0.000000000000000, + 0.000041220652235, + -0.000634590332168, + -0.001785625081886, + 0.001317506937803, + 0.006054485538259, + 0.001737510892287, + 0.000000000000000, + 0.011477784663031, + -0.012227956391224, + -0.076671645271867, + -0.037968922297635, + 0.139356067333892, + 0.175261307377343, + -0.076498342595344, + 0.741352206757141, + -0.076498342595344, + 0.175261307377343, + 0.139356067333892, + -0.037968922297635, + -0.076671645271867, + -0.012227956391224, + 0.011477784663031, + 0.000000000000000, + 0.001737510892287, + 0.006054485538259, + 0.001317506937803, + -0.001785625081886, + -0.000634590332168, + 0.000041220652235, + -0.000000000000000, + ] + ) + verify_impulse_response(h, h_truth, atol=1e-2) + + +def test_blackman_harris(): + """ + MATLAB: + >> h = designBandstopFIR(FilterOrder=30, CenterFrequency=0.4, Bandwidth=0.25, Window="blackman-harris"); + >> transpose(h) + """ + h = sdr.design_bandstop_fir(30, 0.4, 0.25, window="blackman-harris") + h_truth = np.array( + [ + 0.000001060796132, + 0.000008076310132, + -0.000143451692851, + -0.000513101416779, + 0.000480260735282, + 0.002726738086864, + 0.000937988377148, + 0.000000000000000, + 0.008213858159922, + -0.009739012408174, + -0.066659327027275, + -0.035402070593523, + 0.137060855858530, + 0.178985506259001, + -0.079892784955629, + 0.727860137412237, + -0.079892784955629, + 0.178985506259001, + 0.137060855858530, + -0.035402070593523, + -0.066659327027275, + -0.009739012408174, + 0.008213858159922, + 0.000000000000000, + 0.000937988377148, + 0.002726738086864, + 0.000480260735282, + -0.000513101416779, + -0.000143451692851, + 0.000008076310132, + 0.000001060796132, + ] + ) + verify_impulse_response(h, h_truth, atol=1e-1) + + +def test_chebyshev(): + """ + MATLAB: + >> h = designBandstopFIR(FilterOrder=30, CenterFrequency=0.4, Bandwidth=0.25, Window="chebyshev"); + >> transpose(h) + """ + h = sdr.design_bandstop_fir(30, 0.4, 0.25, window="chebyshev") + h_truth = np.array( + [ + 0.000309113441909, + 0.000349144761797, + -0.002354737618114, + -0.004540140197669, + 0.002658020078420, + 0.010356031164874, + 0.002615046535495, + 0.000000000000000, + 0.014316841728137, + -0.014231877095468, + -0.084372501258072, + -0.039973462497828, + 0.141881452704038, + 0.174303969847392, + -0.075032265565100, + 0.747471803375481, + -0.075032265565100, + 0.174303969847392, + 0.141881452704038, + -0.039973462497828, + -0.084372501258072, + -0.014231877095468, + 0.014316841728137, + 0.000000000000000, + 0.002615046535495, + 0.010356031164874, + 0.002658020078420, + -0.004540140197669, + -0.002354737618114, + 0.000349144761797, + 0.000309113441909, + ] + ) + verify_impulse_response(h, h_truth, atol=1e-2) + + +def test_kaiser(): + """ + MATLAB: + >> h = designBandstopFIR(FilterOrder=30, CenterFrequency=0.4, Bandwidth=0.25, Window="kaiser"); + >> transpose(h) + """ + h = sdr.design_bandstop_fir(30, 0.4, 0.25, window="kaiser", atten=20) + h_truth = np.array( + [ + 0.016765450319470, + 0.010339454228757, + -0.038373027043513, + -0.045306490646877, + 0.017552083072319, + 0.048095293336892, + 0.008983538311274, + 0.000000000000000, + 0.030503119574613, + -0.025205246008597, + -0.128266302246497, + -0.053758312599290, + 0.173698261663486, + 0.199683408020401, + -0.082623297525045, + 0.725553680723470, + -0.082623297525045, + 0.199683408020401, + 0.173698261663486, + -0.053758312599290, + -0.128266302246497, + -0.025205246008597, + 0.030503119574613, + 0.000000000000000, + 0.008983538311274, + 0.048095293336892, + 0.017552083072319, + -0.045306490646877, + -0.038373027043513, + 0.010339454228757, + 0.016765450319470, + ] + ) + verify_impulse_response(h, h_truth, atol=1e-1)