"""Utility functions for filtering."""
import numpy as np
from scipy.signal import freqz, sosfreqz
from neurodsp.utils.decorators import multidim
from neurodsp.filt.checks import check_filter_definition
###################################################################################################
###################################################################################################
def infer_passtype(f_range):
"""Given frequency definition of a filter, infer the passtype.
Parameters
----------
f_range : tuple of (float, float)
Cutoff frequency(ies) used for filter, specified as f_lo & f_hi.
Returns
-------
pass_type : str
Which kind of filter pass_type is consistent with the frequency definition provided.
Notes
-----
Assumes that a definition with two frequencies is a 'bandpass' (not 'bandstop').
"""
if f_range[0] is None:
pass_type = 'lowpass'
elif f_range[1] is None:
pass_type = 'highpass'
else:
pass_type = 'bandpass'
# Check the inferred passtype & frequency definition is valid
_ = check_filter_definition(pass_type, f_range)
return pass_type
[docs]def compute_frequency_response(filter_coefs, a_vals, fs):
"""Compute the frequency response of a filter.
Parameters
----------
filter_coefs : 1d or 2d array
If 1d, interpreted as the B-value filter coefficients.
If 2d, interpreted as the second-order (sos) filter coefficients.
a_vals : 1d array or None
The A-value filter coefficients for a filter.
If second-order filter coefficients are provided in `filter_coefs`, must be None.
fs : float
Sampling rate, in Hz.
Returns
-------
f_db : 1d array
Frequency vector corresponding to attenuation decibels, in Hz.
db : 1d array
Degree of attenuation for each frequency specified in `f_db`, in dB.
Examples
--------
Compute the frequency response for an FIR filter:
>>> from neurodsp.filt.fir import design_fir_filter
>>> filter_coefs = design_fir_filter(fs=500, pass_type='bandpass', f_range=(8, 12))
>>> f_db, db = compute_frequency_response(filter_coefs, 1, fs=500)
Compute the frequency response for an IIR filter, which uses SOS coefficients:
>>> from neurodsp.filt.iir import design_iir_filter
>>> sos_coefs = design_iir_filter(fs=500, pass_type='bandpass',
... f_range=(8, 12), butterworth_order=3)
>>> f_db, db = compute_frequency_response(sos_coefs, None, fs=500)
"""
if filter_coefs.ndim == 1 and a_vals is not None:
# Compute response for B & A value filter coefficient inputs
w_vals, h_vals = freqz(filter_coefs, a_vals, worN=int(fs * 2))
elif filter_coefs.ndim == 2 and a_vals is None:
# Compute response for sos filter coefficient inputs
w_vals, h_vals = sosfreqz(filter_coefs, worN=int(fs * 2))
else:
raise ValueError("The organization of the filter coefficient inputs is not understood.")
f_db = w_vals * fs / (2. * np.pi)
db = 20 * np.log10(abs(h_vals))
return f_db, db
[docs]def compute_pass_band(fs, pass_type, f_range):
"""Compute the pass bandwidth of a filter.
Parameters
----------
fs : float
Sampling rate, in Hz.
pass_type : {'bandpass', 'bandstop', 'lowpass', 'highpass'}
Which kind of filter to apply:
* 'bandpass': apply a bandpass filter
* 'bandstop': apply a bandstop (notch) filter
* 'lowpass': apply a lowpass filter
* 'highpass' : apply a highpass filter
f_range : tuple of (float, float) or float
Cutoff frequency(ies) used for filter, specified as f_lo & f_hi.
For 'bandpass' & 'bandstop', must be a tuple.
For 'lowpass' or 'highpass', can be a float that specifies pass frequency, or can be
a tuple and is assumed to be (None, f_hi) for 'lowpass', and (f_lo, None) for 'highpass'.
Returns
-------
pass_bw : float
The pass bandwidth of the filter, in Hz.
Examples
--------
Compute the bandwidth of a bandpass filter:
>>> compute_pass_band(fs=500, pass_type='bandpass', f_range=(5, 25))
20.0
"""
f_lo, f_hi = check_filter_definition(pass_type, f_range)
if pass_type in ['bandpass', 'bandstop']:
pass_bw = f_hi - f_lo
elif pass_type == 'highpass':
pass_bw = compute_nyquist(fs) - f_lo
elif pass_type == 'lowpass':
pass_bw = f_hi
return pass_bw
[docs]def compute_transition_band(f_db, db, low=-20, high=-3):
"""Compute transition bandwidth of a filter.
Parameters
----------
f_db : 1d array
Frequency vector corresponding to attenuation decibels, in Hz.
db : 1d array
Degree of attenuation for each frequency specified in `f_db`, in dB.
low : float, optional, default: -20
The lower limit that defines the transition band, in dB.
high : float, optional, default: -3
The upper limit that defines the transition band, in dB.
Returns
-------
transition_band : float
The transition bandwidth of the filter, in Hz.
Examples
--------
Compute the transition band of an FIR filter, using the computed frequency response:
>>> from neurodsp.filt.fir import design_fir_filter
>>> filter_coefs = design_fir_filter(fs=500, pass_type='bandpass', f_range=(1, 25))
>>> f_db, db = compute_frequency_response(filter_coefs, 1, fs=500)
>>> compute_transition_band(f_db, db, low=-20, high=-3)
0.5
Compute the transition band of an IIR filter, using the computed frequency response:
>>> from neurodsp.filt.iir import design_iir_filter
>>> sos = design_iir_filter(fs=500, pass_type='bandstop',
... f_range=(10, 20), butterworth_order=7)
>>> f_db, db = compute_frequency_response(sos, None, fs=500)
>>> compute_transition_band(f_db, db, low=-20, high=-3)
2.0
"""
# This gets the indices of transitions to the values in searched for range
inds = np.where(np.diff(np.logical_and(db > low, db < high)))[0]
if len(inds) == 0:
raise ValueError("Invalid transition band. Redefine f_range.")
else:
# This steps through the indices, in pairs, selecting from the vector to select from
transition_band = np.max([(b - a) for a, b in zip(f_db[inds[0::2]], f_db[inds[1::2]])])
return transition_band
[docs]def compute_nyquist(fs):
"""Compute the Nyquist frequency.
Parameters
----------
fs : float
Sampling rate, in Hz.
Returns
-------
float
The Nyquist frequency of a signal with the given sampling rate, in Hz.
Examples
--------
Compute the Nyquist frequency for a 500 Hz sampling rate:
>>> compute_nyquist(fs=500)
250.0
"""
return fs / 2.
[docs]@multidim()
def remove_filter_edges(sig, filt_len):
"""Drop the edges, by making NaN, from a filtered signal, to avoid edge artifacts.
Parameters
----------
sig : 1d array
Filtered signal to have edge artifacts removed from.
filt_len : int
Length of the filter that was applied.
Returns
-------
sig : 1d array
Filter signal with edge artifacts switched to NaNs.
Examples
--------
Apply a filter and remove the filter edges of the filtered signal:
>>> from neurodsp.filt.fir import design_fir_filter, apply_fir_filter
>>> from neurodsp.sim import sim_combined
>>> sig = sim_combined(n_seconds=10, fs=500,
... components={'sim_powerlaw': {}, 'sim_oscillation' : {'freq': 10}})
>>> filter_coefs = design_fir_filter(fs=500, pass_type='bandpass', f_range=(1, 25))
>>> filt_sig = apply_fir_filter(sig, filter_coefs)
>>> filt_sig_no_edges = remove_filter_edges(filt_sig, filt_len=len(filter_coefs))
"""
n_rmv = int(np.ceil(filt_len / 2))
sig[:n_rmv] = np.nan
sig[-n_rmv:] = np.nan
return sig