Filter to emulate background in human intracranial microwire recordings.

StdNoiseFilter.py

from math import pi as PI
import warnings

import numpy as np
import pandas as pd
from scipy import signal

class StdNoiseFilter:
    """
    IIR filter to provide standard noise implementation given sampling rate in Hz.

    Emulates observed sample channel characteristics when filtering white noise.
    """
    def __init__(self, sampRate = 24000, gainFactor = 1):
        """
        Construct filter from sampling rate and gain factor
        :param sampRate:
        :param gainFactor:
        """
        self._sr = sampRate

        pars = np.array([0.494 * gainFactor, sampRate, 2.1990, 2.2003, 2.2010, 9.999997, 100, 100, 533.5, 2822, 5021])

        # sos form with arbitrary gain empirically chosen to produce 1e-6 std in
        # passband for 1e-6 std input with 24,000 sampling rate
        self.sos = StdNoiseFilter.buildFromPars(pars)

    def buildFromPars(pars):
        """
        Build an sos filter from specified gain and corners.

        Used primarily for optimzation against known spectra. Parameters are provided as
        array for this use.

        :param pars: parameters for construction. np.array with 11 elements:
            0: gain, filter design gain
            1: sr, sampling rate (samples/s)
            2: f0, lowest frequency of high pass filter (Hz)
            3: p1, frequency to start decrease of gain near top N=-1
            4: p2, frequency to decrease gain further, N=-2
            5: z3, frequency to start decrease of slope, N=-1
            6: z4, frequency to decrease slope even further, N=0
            7: p5, frequency to start turn down of slope, N=-1
            8: z6, frequency to decrease slope again, roughly bottom of passband, N=0
            9: p7, frequency to turn down slope, roughly top of passband, N=-1
            10: f1, highest frequency for corner of low pass filter
        :return: sos filter in standard form
        """
        nyqFreq = pars[1] / 2

        # high pass filter to produce bump at low frequencies
        z, p, k = signal.butter(1, pars[2] / nyqFreq, 'high', output='zpk')  # frequency in rad / sample
        # output pole lies at 1 - f0 / F_nyq * pi rad / sample
        # single zero at 1 lies to right of all frequencies and determines slope at
        # lower frequencies

        # pole to start decrease in gain, N = -1
        p1 = 1 - pars[3] / nyqFreq * PI

        # pole to decrease gain a bit more, N = -2
        p2 = 1 - pars[4] / nyqFreq * PI

        # zero to decrease slope, N = -1
        z3 = 1 - pars[5] / nyqFreq * PI

        # zero to decrease slope, N = 0
        z4 = 1 - pars[6] / nyqFreq * PI

        # pole to decrease gain, N = -1
        p5 = 1 - pars[7] / nyqFreq * PI

        # zero to decrease slope, N = 0
        z6 = 1 - pars[8] / nyqFreq * PI

        # pole to decrease gain, N = -1
        p7 = 1 - pars[9] / nyqFreq * PI

        # butterworth low pass for upper anti - aliasing equivalent at 4000 Hz
        z10, p10, k10 = signal.butter(1, 4000 / nyqFreq, output='zpk')

        # combine poles and zeros
        zf = (z[0], z3, z4, z6, z10[0])
        pf = (np.real(p[0]), p1, p2, p5, p7, np.real(p10[0]))

        # sos form with arbitrary gain empirically chosen to produce 1e-6 std in
        # passband for 1e-6 std input with 24,000 sampling rate
        sos = signal.zpk2sos(zf, pf, pars[0])

        # hack on dc coefficient at 0,0 and coefficient at 0,2
        # sos[0,0] = 0
        # sos[0,2] = 1

        return sos

    def impResp(self, numSamps=5000):
        """
        Calculate impulse response of filter.

        :param numSamps number of samples to compute for impulse response
        :return: impulse response as pandas. Series with time in seconds as index
        """
        y = np.zeros(numSamps)
        y[0] = 1.0
        y = signal.sosfilt(self.sos, y)

        res = pd.Series(y)
        res.index = res.index / self._sr

        return res

    def filtLength(self, endRatio=0.001):
        """
        Compute effective length of filter by determining sample at which impulse
        response decays to given ratio from maximum absolute value (default = 0.001).
        Looks at final decay, ignoring any transient zero crossings.

        :return: effective length of filter
        """
        checkLen = 5000
        found = False
        while (not found and checkLen < 1e6):
            impr = self.impResp(checkLen)
            absv = abs(impr.values)
            maxAV = max(absv)
            lastGTI = round(self._sr * max(impr.index[absv >= maxAV * endRatio]))
            if lastGTI < impr.values.size-1:
                found = True
            else:
                checkLen *= 2

        if not found:
            warnings.warn("Could not find end of impulse response in {} samples.".format(checkLen))
            return np.nan
        else:
            return lastGTI