Skip to content

Instantly share code, notes, and snippets.

@PeterNSteinmetz
Created August 17, 2023 18:15
Show Gist options
  • Save PeterNSteinmetz/47050c4796f9fd818d0a6b5d0a2f0d68 to your computer and use it in GitHub Desktop.
Save PeterNSteinmetz/47050c4796f9fd818d0a6b5d0a2f0d68 to your computer and use it in GitHub Desktop.
Filter to emulate background in human intracranial microwire recordings.
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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment