# generate linear time-domain filter coefficients, common to both H1 and L1. # First, define some functions: # This function will generate digital filter coefficients for bandstops (notches). # Understanding it requires some signal processing expertise, which we won't get into here. def iir_bandstops ( fstops , fs , order = 4 ): """ellip notch filter fstops is a list of entries of the form [frequency (Hz), df, df2] where df is the pass width and df2 is the stop width (narrower than the pass width). Use caution if passing more than one freq at a time, because the filter response might behave in ways you don't expect. """ nyq = 0.5 * fs # Zeros zd, poles pd, and gain kd for the digital filter zd = np . array ([]) pd = np . array ([]) kd = 1 # Notches for fstopData in fstops : fstop = fstopData [ 0 ] df = fstopData [ 1 ] df2 = fstopData [ 2 ] low = ( fstop - df ) / nyq high = ( fstop + df ) / nyq low2 = ( fstop - df2 ) / nyq high2 = ( fstop + df2 ) / nyq z , p , k = iirdesign ([ low , high ], [ low2 , high2 ], gpass = 1 , gstop = 6 , ftype = 'ellip' , output = 'zpk' ) zd = np . append ( zd , z ) pd = np . append ( pd , p ) # Set gain to one at 100 Hz...better not notch there bPrelim , aPrelim = zpk2tf ( zd , pd , 1 ) outFreq , outg0 = freqz ( bPrelim , aPrelim , 100 / nyq ) # Return the numerator and denominator of the digital filter b , a = zpk2tf ( zd , pd , k ) return b , a def get_filter_coefs ( fs ): # assemble the filter b,a coefficients: coefs = [] # bandpass filter parameters lowcut = 43 highcut = 260 order = 4 # bandpass filter coefficients nyq = 0.5 * fs low = lowcut / nyq high = highcut / nyq bb , ab = butter ( order , [ low , high ], btype = 'band' ) coefs . append (( bb , ab )) # Frequencies of notches at known instrumental spectral line frequencies. # You can see these lines in the ASD above, so it is straightforward to make this list. notchesAbsolute = np . array ( [ 14.0 , 34.70 , 35.30 , 35.90 , 36.70 , 37.30 , 40.95 , 60.00 , 120.00 , 179.99 , 304.99 , 331.49 , 510.02 , 1009.99 ]) # notch filter coefficients: for notchf in notchesAbsolute : bn , an = iir_bandstops ( np . array ([[ notchf , 1 , 0.1 ]]), fs , order = 4 ) coefs . append (( bn , an )) # Manually do a wider notch filter around 510 Hz etc. bn , an = iir_bandstops ( np . array ([[ 510 , 200 , 20 ]]), fs , order = 4 ) coefs . append (( bn , an )) # also notch out the forest of lines around 331.5 Hz bn , an = iir_bandstops ( np . array ([[ 331.5 , 10 , 1 ]]), fs , order = 4 ) coefs . append (( bn , an )) return coefs # and then define the filter function: def filter_data ( data_in , coefs ): data = data_in . copy () for coef in coefs : b , a = coef # filtfilt applies a linear filter twice, once forward and once backwards. # The combined filter has linear phase. data = filtfilt ( b , a , data ) return data