mne_c_filter.py

import numpy as np
import matplotlib.pyplot as plt


def _create_filter_response(sfreq, lowpass, highpass, lowpass_width=5.,
                            highpass_width=0.5):
    """Create a frequency response"""
    filter_size = 4096
    filter_taper_size = 2048
    resp_size = (filter_size + 2 * filter_taper_size) // 2 + 1

    freq_resp = np.ones(resp_size)
    # Start simple first
    highpasss = int(((resp_size - 1) * highpass) / (0.5 * sfreq))
    lowpasss = int(((resp_size - 1) * lowpass) / (0.5 * sfreq))

    lowpass_widths = int(((resp_size - 1) * lowpass_width) / (0.5 * sfreq))
    lowpass_widths = (lowpass_widths + 1) // 2        # What user specified
    highpass_widths = int(((resp_size-1) * highpass_width) / (0.5 * sfreq))
    highpass_widths = (highpass_widths + 1) // 2
    if highpasss > highpass_widths + 1:
        w = highpass_widths
        mult = 1.0 / w
        add = 3.0
        s = highpasss - w + 1
        freq_resp[:s] = 0.
        for k in range(-w + 1, w):
            if s >= 0 and s < resp_size:
                c = np.cos(np.pi / 4. * (k * mult + add))
                freq_resp[s] *= c * c
            s += 1
    else:
        raise NotImplementedError

    w = lowpass_widths
    mult = 1.0/w
    add = 1.0
    s = lowpasss-w+1
    for k in range(-w + 1, w):
        if s >= 0 and s < resp_size:
            c = np.cos(np.pi / 4. * (k * mult + add))
            freq_resp[s] *= c*c
        s += 1
    freq_resp[s:] = 0.0
    return freq_resp


freq_resp = _create_filter_response(1000., 40., 1)
plt.plot(freq_resp)