Filtering#

The filtering module holds the functions for frequency manipulation :

nap.apply_bandstop_filter
nap.apply_lowpass_filter
nap.apply_highpass_filter
nap.apply_bandpass_filter

The functions have similar calling signatures. For example, to filter a 1000 Hz signal between 10 and 20 Hz using a Butterworth filter:

>>> new_tsd = nap.apply_bandpass_filter(tsd, (10, 20), fs=1000, mode='butter')

Currently, the filtering module provides two methods for frequency manipulation: butter for a recursive Butterworth filter and sinc for a Windowed-sinc convolution. This notebook provides a comparison of the two methods.

Basics#

We start by generating a signal with multiple frequencies (2, 10 and 50 Hz).

fs = 1000 # sampling frequency
t = np.linspace(0, 2, fs * 2)
f2 = np.cos(t*2*np.pi*2)
f10 = np.cos(t*2*np.pi*10)
f50 = np.cos(t*2*np.pi*50)

sig = nap.Tsd(t=t,d=f2+f10+f50 + np.random.normal(0, 0.5, len(t)))

Text(0.5, 0, 'Time (s)')

../_images/e258617a6cc7914aa032c5044271d8d6fda45f480989a5c798b0dad84658293e.png

We can compute the Fourier transform of sig to verify that all the frequencies are there.

psd = nap.compute_power_spectral_density(sig, fs, norm=True)

(0.0, 100.0)

../_images/66d7647b40b54230dfd5a21f7056ee8150eab403a3a0a66cb4d61dfe6d86406e.png

Let’s say we would like to see only the 10 Hz component. We can use the function apply_bandpass_filter with mode butter for Butterworth.

sig_butter = nap.apply_bandpass_filter(sig, (8, 12), fs, mode='butter')

Let’s compare it to the sinc mode for Windowed-sinc.

sig_sinc = nap.apply_bandpass_filter(sig, (8, 12), fs, mode='sinc', transition_bandwidth=0.003)

(0.0, 1.0)

../_images/2ebf2b10b97ae75dc2d4a6551844756907e276ab840bbe033f3cb9a9f53f5325.png

This gives similar results except at the edges.

Another use of filtering is to remove some frequencies. Here we can try to remove the 50 Hz component in the signal.

sig_butter = nap.apply_bandstop_filter(sig, cutoff=(45, 55), fs=fs, mode='butter')
sig_sinc = nap.apply_bandstop_filter(sig, cutoff=(45, 55), fs=fs, mode='sinc', transition_bandwidth=0.004)

(0.0, 1.0)

../_images/9d6a0e3a2efa813d48f7bd77b4d547be57479c54e37a60f9b183976ed33b5045.png

Let’s see what frequencies remain;

psd_butter = nap.compute_power_spectral_density(sig_butter, fs, norm=True)
psd_sinc = nap.compute_power_spectral_density(sig_sinc, fs, norm=True)

(0.0, 70.0)

../_images/bace54b0bcf5a0931c4a6eb93d2bc4ff814c41fd177cd9f43d116ee4938717f4.png

Inspecting frequency fesponses of a filter#

We can inspect the frequency response of a filter by plotting its power spectral density (PSD). To do this, we can use the get_filter_frequency_response function, which returns a pandas Series with the frequencies as the index and the PSD as values.

Let’s extract the frequency response of a Butterworth filter and a sinc low-pass filter.

# compute the frequency response of the filters
psd_butter = nap.get_filter_frequency_response(
    200, fs,"lowpass", "butter", order=8
)
psd_sinc = nap.get_filter_frequency_response(
    200, fs,"lowpass", "sinc", transition_bandwidth=0.1
)

Transition band butterworth filter: (164Hz, 290Hz)
Transition band sinc filter: (146Hz, 268Hz)

../_images/c29dafc651da1d8f612ed4b2674a360ab908bf7ffff76b0ee1d94b535ae62199.png

The frequency band with response close to one will be preserved by the filtering (pass band), the band with response close to zero will be discarded (stop band), and the band in between will be partially attenuated (transition band).

Note

Here, we define the transition band as the range where the amplitude attenuation is between 99% and 1%. The transition_bandwidth parameter of the sinc filter is approximately the width of the transition band normalized by the sampling frequency. In the example above, if you divide the transition band width of 122Hz by the sampling frequency of 1000Hz, you get 0.122, which is close to the 0.1 value set.

You can modulate the width of the transition band by setting the order parameter of the Butterworth filter or the transition_bandwidth parameter of the sinc filter. First, let’s get the frequency response for a Butterworth low pass filter with different order:

butter_freq = {
    order: nap.get_filter_frequency_response(250, fs, "lowpass", "butter", order=order)
    for order in [2, 4, 6]}

… and then the frequency response for the Windowed-sinc equivalent with different transition bandwidth.

sinc_freq = {
    tb: nap.get_filter_frequency_response(250, fs,"lowpass", "sinc", transition_bandwidth=tb)
    for tb in [0.002, 0.02, 0.2]}

Let’s plot the frequency response of both.

../_images/7b95d5a0d890a2b2bae4c1c4a66e4aa5c6893efcba7009ddedec7e8f673e4ced.png

Warning

In some cases, the transition bandwidth that is too high generates a kernel that is too short. The amplitude of the original signal will then be lower than expected. In this case, the solution is to decrease the transition bandwidth when using the windowed-sinc mode. Note that this increases the length of the kernel significantly. Let see it with the band pass filter.

sinc_freq = {
    tb:nap.get_filter_frequency_response((100, 200), fs, "bandpass", "sinc", transition_bandwidth=tb)
    for tb in [0.004, 0.2]}

../_images/9a9edc6d970bb653204b1e4085962aea1f6f2d4352c9dd3f0bfecfe80fa34580.png

Performances#

Let’s compare the performance of each when varying the number of time points and the number of dimensions.

Show code cell source Hide code cell source

from time import perf_counter

def get_mean_perf(tsd, mode, n=10):
    tmp = np.zeros(n)
    for i in range(n):
        t1 = perf_counter()
        _ = nap.apply_lowpass_filter(tsd, 0.25 * tsd.rate, mode=mode)
        t2 = perf_counter()
        tmp[i] = t2 - t1
    return [np.mean(tmp), np.std(tmp)]

def benchmark_time_points(mode):
    times = []
    for T in np.arange(1000, 100000, 20000):
        time_array = np.arange(T)/1000
        data_array = np.random.randn(len(time_array))
        startend = np.linspace(0, time_array[-1], T//100).reshape(T//200, 2)
        ep = nap.IntervalSet(start=startend[::2,0], end=startend[::2,1])
        tsd = nap.Tsd(t=time_array, d=data_array, time_support=ep)
        times.append([T]+get_mean_perf(tsd, mode))
    return np.array(times)

def benchmark_dimensions(mode):
    times = []
    for n in np.arange(1, 100, 10):
        time_array = np.arange(10000)/1000
        data_array = np.random.randn(len(time_array), n)
        startend = np.linspace(0, time_array[-1], 10000//100).reshape(10000//200, 2)
        ep = nap.IntervalSet(start=startend[::2,0], end=startend[::2,1])
        tsd = nap.TsdFrame(t=time_array, d=data_array, time_support=ep)
        times.append([n]+get_mean_perf(tsd, mode))
    return np.array(times)


times_sinc = benchmark_time_points(mode="sinc")
times_butter = benchmark_time_points(mode="butter")

dims_sinc = benchmark_dimensions(mode="sinc")
dims_butter = benchmark_dimensions(mode="butter")


plt.figure(figsize = (16, 5))
plt.subplot(121)
for arr, label in zip(
    [times_sinc, times_butter],
    ["Windowed-sinc", "Butter"],
    ):
    plt.plot(arr[:, 0], arr[:, 1], "o-", label=label)
    plt.fill_between(arr[:, 0], arr[:, 1] - arr[:, 2], arr[:, 1] + arr[:, 2], alpha=0.2)
plt.legend()
plt.xlabel("Number of time points")
plt.ylabel("Time (s)")
plt.title("Low pass filtering benchmark")
plt.subplot(122)
for arr, label in zip(
    [dims_sinc, dims_butter],
    ["Windowed-sinc", "Butter"],
    ):
    plt.plot(arr[:, 0], arr[:, 1], "o-", label=label)
    plt.fill_between(arr[:, 0], arr[:, 1] - arr[:, 2], arr[:, 1] + arr[:, 2], alpha=0.2)
plt.legend()
plt.xlabel("Number of dimensions")
plt.ylabel("Time (s)")
plt.title("Low pass filtering benchmark")

Text(0.5, 1.0, 'Low pass filtering benchmark')

../_images/07d65abb1d5971018995b849693402064efd95b013d132d490406002f55f42dd.png