diff --git a/.github/workflows/run_unix.yml b/.github/workflows/run_unix.yml index 557001dae..c003901d9 100644 --- a/.github/workflows/run_unix.yml +++ b/.github/workflows/run_unix.yml @@ -354,6 +354,8 @@ jobs: run: sudo -H python3 $GITHUB_WORKSPACE/python_package/examples/tests/transforms.py - name: Downsampling Python run: sudo -H python3 $GITHUB_WORKSPACE/python_package/examples/tests/downsampling.py + - name: Chebyshev Band Filters Python + run: sudo -H python3 $GITHUB_WORKSPACE/python_package/examples/tests/chebyshev_band_filters.py - name: ICA Python run: sudo -H python3 $GITHUB_WORKSPACE/python_package/examples/tests/ica.py - name: CSP Python diff --git a/python_package/examples/tests/chebyshev_band_filters.py b/python_package/examples/tests/chebyshev_band_filters.py new file mode 100644 index 000000000..777071bc9 --- /dev/null +++ b/python_package/examples/tests/chebyshev_band_filters.py @@ -0,0 +1,56 @@ +import numpy as np + +from brainflow.data_filter import DataFilter, FilterTypes + +SAMPLING_RATE = 256 +PASS_FREQ = 10.0 +STOP_FREQ = 50.0 +CENTER_LOW = 5.0 +CENTER_HIGH = 15.0 +ORDER = 4 +RIPPLE = 1.0 + + +def amplitude_at(signal, freq): + spectrum = np.abs(np.fft.rfft(signal)) * 2 / len(signal) + return spectrum[int(round(freq * len(signal) / SAMPLING_RATE))] + + +def main(): + samples = SAMPLING_RATE * 4 + t = np.arange(samples) / SAMPLING_RATE + # One tone inside the 5 to 15 Hz band and one well outside it. + signal = np.sin(2 * np.pi * PASS_FREQ * t) + np.sin(2 * np.pi * STOP_FREQ * t) + + assert np.isclose(amplitude_at(signal, PASS_FREQ), 1.0, atol=0.01) + assert np.isclose(amplitude_at(signal, STOP_FREQ), 1.0, atol=0.01) + + # Chebyshev type 1 band pass and band stop take five design parameters and expect the + # ripple after the band width. Writing it into slot 3 overwrote the band width, which + # produced an unusable design and an all-NaN output. + for name, apply_filter, kept, removed in ( + ('bandpass', DataFilter.perform_bandpass, PASS_FREQ, STOP_FREQ), + ('bandstop', DataFilter.perform_bandstop, STOP_FREQ, PASS_FREQ), + ): + filtered = np.copy(signal) + apply_filter( + filtered, + SAMPLING_RATE, + CENTER_LOW, + CENTER_HIGH, + ORDER, + FilterTypes.CHEBYSHEV_TYPE_1.value, + RIPPLE, + ) + + assert np.all(np.isfinite(filtered)), '%s produced non-finite samples' % name + # The tone the filter is meant to keep survives close to full amplitude. + assert amplitude_at(filtered, kept) > 0.8, (name, amplitude_at(filtered, kept)) + # The tone it is meant to reject is pushed far down. + assert amplitude_at(filtered, removed) < 0.05, (name, amplitude_at(filtered, removed)) + + print('chebyshev band filter regression passed') + + +if __name__ == '__main__': + main() diff --git a/src/data_handler/data_handler.cpp b/src/data_handler/data_handler.cpp index 106c614dd..a52b0db3e 100644 --- a/src/data_handler/data_handler.cpp +++ b/src/data_handler/data_handler.cpp @@ -294,14 +294,17 @@ int perform_bandpass (double *data, int data_len, int sampling_rate, double star } Dsp::Params params; + params.clear (); params[0] = sampling_rate; // sample rate params[1] = order; // order params[2] = center_freq; // center freq - params[3] = band_width; + params[3] = band_width; // band width if ((filter_type == (int)FilterTypes::CHEBYSHEV_TYPE_1) || (filter_type == (int)FilterTypes::CHEBYSHEV_TYPE_1_ZERO_PHASE)) { - params[3] = ripple; // ripple + // band pass and band stop chebyshev designs take 5 params and expect the ripple after + // the band width, unlike the low pass and high pass designs which take 4 + params[4] = ripple; // ripple } f->setParams (params); f->process (data_len, filter_data); @@ -362,14 +365,17 @@ int perform_bandstop (double *data, int data_len, int sampling_rate, double star } Dsp::Params params; + params.clear (); params[0] = sampling_rate; // sample rate params[1] = order; // order params[2] = center_freq; // center freq - params[3] = band_width; + params[3] = band_width; // band width if ((filter_type == (int)FilterTypes::CHEBYSHEV_TYPE_1) || (filter_type == (int)FilterTypes::CHEBYSHEV_TYPE_1_ZERO_PHASE)) { - params[3] = ripple; // ripple + // band pass and band stop chebyshev designs take 5 params and expect the ripple after + // the band width, unlike the low pass and high pass designs which take 4 + params[4] = ripple; // ripple } f->setParams (params); f->process (data_len, filter_data);