Coverage for C: \ Users \ peaco \ OneDrive \ Documents \ GitHub \ mth5 \ mth5 \ timeseries \ ts_filters.py: 57%
265 statements
« prev ^ index » next coverage.py v7.13.1, created at 2026-01-10 00:01 -0800
1#!/usr/bin/env python
3"""
4time series filters
6"""
8# =================================================================
9import numpy as np
10from loguru import logger
11from matplotlib import pyplot as plt
12from matplotlib.lines import Line2D
13from scipy import signal
16# =================================================================
19def butter_bandpass(lowcut, highcut, sample_rate, order=5):
20 """
21 Butterworth bandpass filter using scipy.signal
23 Transforms band corners to angular frequencies
25 :param lowcut: low cut frequency in Hz (3dB point)
26 :type lowcut: float
27 :param highcut: high cut frequency in Hz (3dB point)
28 :type highcut: float
29 :param sample_rate: Sample rate
30 :type sample_rate: float
31 :param order: Butterworth order, defaults to 5
32 :type order: int, optional
33 :return: SOS scipy.signal format
34 :rtype: scipy.signal.SOS?
36 """
37 nyq = 0.5 * sample_rate
38 if (highcut is None) and (lowcut is None):
39 msg = f"Butterworth bandpass undefined with edges ({lowcut}, {highcut})\n"
40 raise ValueError(msg)
42 # Transforms band corners to angular frequencies
43 if lowcut is not None:
44 low = lowcut / nyq
45 if highcut is not None:
46 high = highcut / nyq
47 if lowcut and highcut:
48 sos = signal.butter(
49 order, [low, high], analog=False, btype="bandpass", output="sos"
50 )
51 return sos
53 msg = f"Butterworth bandpass requested with edges ({lowcut}, {highcut})\n"
54 if highcut is None:
55 msg += "Upper band edge not defined, will treat as a High Pass Filter\n"
56 logger.info(msg)
57 return signal.butter(order, low, analog=False, btype="highpass", output="sos")
58 elif lowcut is None:
59 msg += "Lower band edge not defined, will treat as a Low Pass Filter\n"
60 logger.info(msg)
61 return signal.butter(order, high, analog=False, btype="lowpass", output="sos")
64def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
65 """
67 :param data: 1D time series data
68 :type data: np.ndarray
69 :param lowcut: low cut frequency in Hz
70 :type lowcut: float
71 :param highcut: high cut frequency in Hz
72 :type highcut: float
73 :param fs: Sample rate
74 :type fs: float
75 :param order: Butterworth order, defaults to 5
76 :type order: int, optional
77 :return: filtered data
78 :rtype: np.ndarray
80 """
81 if (highcut is None) and (lowcut is None):
82 msg = f"Butterworth bandpass undefined with edges ({lowcut}, {highcut})\n"
83 msg = f"{msg} Returning original data"
84 logger.warning(msg)
85 return np.copy(data)
87 sos = butter_bandpass(lowcut, highcut, fs, order=order)
88 y = signal.sosfiltfilt(sos, data)
90 return y
93def low_pass(data, low_pass_freq, cutoff_freq, sample_rate):
94 """
96 :param data: 1D time series data
97 :type data: np.ndarray
98 :param low_pass_freq: low pass frequency in Hz
99 :type low_pass_freq: float
100 :param cutoff_freq: cut off frequency in Hz
101 :type cutoff_freq: float
102 :param sampling_rate: Sample rate in samples per second
103 :type sampling_rate: float
104 :return: lowpass filtered data
105 :rtype: np.ndarray
107 """
108 nyq = 0.5 * sample_rate
110 filt_order, wn = signal.buttord(low_pass_freq / nyq, cutoff_freq / nyq, 3, 40)
112 b, a = signal.butter(filt_order, wn, btype="low")
113 data_filtered = signal.filtfilt(b, a, data)
115 return data_filtered
118def zero_pad(input_array, power=2, pad_fill=0):
119 """
121 :param input_array: 1D array
122 :type input_array: np.ndarray
123 :param power: base power to used to pad to, defaults to 2 which is optimal
124 for the FFT
125 :type power: int, optional
126 :param pad_fill: fill value for padded values, defaults to 0
127 :type pad_fill: float, optional
128 :return: zero padded array
129 :rtype: np.ndarray
131 """
133 len_array = input_array.shape[0]
134 if power == 2:
135 npow = int(np.ceil(np.log2(len_array)))
136 if power == 10:
137 npow = int(np.ceil(np.log10(len_array)))
138 if npow > 32:
139 logger.warning(
140 "Exceeding memory allocation inherent in your computer 2**32. "
141 "Limiting the zero pad to 2**32"
142 )
143 pad_array = np.zeros(power**npow)
144 if pad_fill != 0:
145 pad_array[:] = pad_fill
146 pad_array[0:len_array] = input_array
148 return pad_array
151class RemoveInstrumentResponse:
152 """
153 Remove instrument response from the given channel response filter
155 The order of operations is important (if applied):
157 1) detrend
158 2) zero mean
159 3) zero pad
160 4) time window
161 5) frequency window
162 6) remove response
163 7) undo time window
164 8) bandpass
166 :param ts: time series data to remove response from
167 :type ts: np.ndarray((N,) , dtype=float)
168 :param time_array: time index that corresponds to the time series
169 :type time_array: np.ndarray((N,) , dtype=np.datetime[ns])
170 :param sample_interval: seconds per sample (time interval between samples)
171 :type sample_interval: float
172 :param channel_response: Channel response filter with all filters
173 included to convert from counts to physical units
174 :type channel_response: `class`:mt_metadata.timeseries.filters.ChannelResponse`
176 **kwargs**
178 :param plot: to plot the calibration process [ False | True ]
179 :type plot: boolean, default True
180 :param detrend: Remove linar trend of the time series
181 :type detrend: boolean, default True
182 :param zero_mean: Remove the mean of the time series
183 :type zero_mean: boolean, default True
184 :param zero_pad: pad the time series to the next power of 2 for efficiency
185 :type zero_pad: boolean, default True
186 :param t_window: Time domain window name see `scipy.signal.windows` for options
187 :type t_window: string, default None
188 :param t_window_params: Time domain window parameters, parameters can be
189 found in `scipy.signal.windows`
190 :type t_window_params: dictionary
191 :param f_window: Frequency domain window name see `scipy.signal.windows` for options
192 :type f_window: string, default None
193 :param f_window_params: Frequency window parameters, parameters can be
194 found in `scipy.signal.windows`
195 :type f_window_params: dictionary
196 :param bandpass: bandpass frequency and order {"low":, "high":, "order":,}
197 :type bandpass: dictionary
200 """
202 def __init__(
203 self,
204 ts,
205 time_array,
206 sample_interval,
207 channel_response,
208 **kwargs,
209 ):
210 """
212 :param ts:
213 :param time_array:
214 :param sample_interval:
215 :param channel_response:
216 :param filters_to_remove: optional list of specific filters to remove. If not provided, filters will be
217 taken from channel_response.
218 """
219 self.logger = logger
220 self.ts = ts
221 self.time_array = time_array
222 self.sample_interval = sample_interval
223 self.channel_response = channel_response
224 self.plot = False
225 self.detrend = True
226 self.zero_mean = True
227 self.zero_pad = True
228 self.t_window = None
229 self.t_window_params = {}
230 self.f_window = None
231 self.f_window_params = {}
232 self.bandpass = {}
233 self.fig = None
234 self.nrows = None
235 self.subplot_dict = {}
236 self.include_decimation = False
237 self.include_delay = False
239 for key, value in kwargs.items():
240 setattr(self, key, value)
242 def _subplots(self, x, y, color, num, label):
243 """
244 helper function to make subplots for if plotting is desired
246 :param x: x array
247 :type x: np.ndarray
248 :param y: y array
249 :type y: np.ndarray
250 :param color: color of the line
251 :type color: tuple
252 :param num: subplot number
253 :type num: integer
254 :param label: legend label
255 :type label: string
257 """
258 ax_t = self.fig.get_axes()[0]
259 ax_f = self.fig.get_axes()[1]
260 ax = self.fig.add_subplot(self.nrows, 2, num, sharex=ax_t)
261 ax2 = self.fig.add_subplot(self.nrows, 2, num + 1, sharex=ax_f)
263 # plot time series
264 line = Line2D([0], [0], color=color)
265 ax.plot(x, y, color=color)
266 f = np.fft.rfftfreq(x.size, d=self.sample_interval)
268 # plot spectra
269 data = np.fft.rfft(y)
270 ax2.loglog(f, abs(data), color=color)
271 ax.legend(
272 [line],
273 [label],
274 loc="upper left",
275 borderaxespad=0.01,
276 borderpad=0.1,
277 markerscale=0.02,
278 ).set_zorder(1000)
280 @staticmethod
281 def get_window(window, window_params, size):
282 """
283 Get window from scipy.signal
285 :param window: name of the window
286 :type window: string
287 :param window_params: dictionary of window parameters
288 :type window_params: dictionary
289 :param size: number of points in the window
290 :type size: integer
291 :return: window function
292 :rtype: class:`scipy.signal`
294 """
295 return getattr(signal.windows, window)(size, **window_params)
297 def apply_detrend(self, ts):
298 """
299 Detrend time series using scipy.detrend('linear')
301 :param ts: input time series
302 :type ts: np.ndarray
303 :return: detrended time series
304 :rtype: np.ndarray
306 """
307 ts = signal.detrend(np.nan_to_num(ts), type="linear")
309 if self.plot:
310 self._subplots(
311 self.time_array,
312 ts,
313 (0.45, 0.1, 0.5),
314 self.subplot_dict["detrend"],
315 "Detrended",
316 )
317 return ts
319 def apply_zero_mean(self, ts):
320 """
321 Remove the mean from the time series
323 :param ts: input time series
324 :type ts: np.ndarray
325 :return: zero mean time series
326 :rtype: np.ndarray
328 """
329 ts = ts - ts.mean()
331 if self.plot:
332 self._subplots(
333 self.time_array,
334 ts,
335 (0.55, 0.1, 0.4),
336 self.subplot_dict["zero_mean"],
337 "Zero Mean",
338 )
339 return ts
341 def apply_zero_pad(self, ts):
342 """
343 zero pad to power of 2, at the end of the time series to make the
344 FFT more efficient
346 :param ts: input time series
347 :type ts: np.ndarray
348 :return: zero padded time series
349 :rtype: np.ndarray
351 """
353 pad_ts = zero_pad(ts)
355 if self.plot:
356 self._subplots(
357 self.time_array,
358 pad_ts[0 : self.time_array.size],
359 (0.7, 0.1, 0.25),
360 self.subplot_dict["pad"],
361 "Zero Pad",
362 )
363 return pad_ts
365 def apply_t_window(self, ts):
366 """
367 Apply a window in the time domain. Get the available windows from
368 `scipy.signal.windows`
370 :param ts: input time series
371 :type ts: np.ndarray
372 :return: windowed time series
373 :rtype: np.ndarray
375 """
376 w = self.get_window(self.t_window, self.t_window_params, self.ts.size)
377 ts = ts * w
379 if self.plot:
380 self._subplots(
381 self.time_array,
382 ts,
383 (0.7, 0.1, 0.25),
384 self.subplot_dict["t_window"],
385 f"t_window {self.t_window.capitalize()}",
386 )
387 return ts
389 def apply_f_window(self, data):
390 """
391 Apply a frequency domain window. Get the available windows from
392 `scipy.signal.windows`
394 Need to create a window twice the size of the input because we are
395 only taking the rfft which gives just half the spectra
396 and then take only half the window
398 :param data: input spectra
399 :type data: np.ndarray
400 :return: windowed spectra
401 :rtype: np.ndarray
403 """
405 w = self.get_window(self.f_window, self.f_window_params, 2 * data.size)[
406 data.size :
407 ]
408 data = data * w
409 if self.plot:
410 f = np.fft.rfftfreq(2 * data.size, d=self.sample_interval)[1:]
411 ax_t = self.fig.get_axes()[0]
412 ax_f = self.fig.get_axes()[1]
413 ax = self.fig.add_subplot(
414 self.nrows, 2, self.subplot_dict["f_window"], sharex=ax_t
415 )
416 ax2 = self.fig.add_subplot(
417 self.nrows, 2, self.subplot_dict["f_window"] + 1, sharex=ax_f
418 )
420 ax.plot(
421 self.time_array,
422 abs(np.fft.irfft(data)[0 : self.time_array.size]),
423 color=(0.85, 0.1, 0.15),
424 )
425 ax2.loglog(f, abs(data), color=(0.85, 0.1, 0.15))
426 axt2 = ax2.twinx()
427 axt2.semilogx(f, w, color=(0.5, 0.5, 0.5), zorder=0)
428 line = Line2D([0], [0], color=(0.85, 0.1, 0.15))
429 ax.legend(
430 [line],
431 [f"Freq Window {self.f_window.capitalize()}"],
432 loc="upper left",
433 borderaxespad=0.01,
434 borderpad=0.1,
435 )
436 return data
438 def apply_bandpass(self, ts):
439 """
440 apply a bandpass filter to the calibrated data
443 :param ts: calibrated time series
444 :type ts: np.ndarray
445 :return: bandpassed time series
446 :rtype: np.ndarray
448 """
449 try:
450 filter_order = self.bandpass["order"]
451 except KeyError:
452 filter_order = 5
453 ts = butter_bandpass_filter(
454 ts,
455 self.bandpass["low"],
456 self.bandpass["high"],
457 1.0 / self.sample_interval,
458 order=filter_order,
459 )
461 if self.plot:
462 self._subplots(
463 self.time_array,
464 ts,
465 (0.875, 0.1, 0.1),
466 self.subplot_dict["bandpass"],
467 "Band Pass",
468 )
469 return ts
471 def _get_subplot_count(self):
472 """
473 helper function to get subplot information
475 :return: dictionary of subplot information
476 :rtype: dictionary
478 """
479 order = [
480 "detrend",
481 "zero_mean",
482 "t_window",
483 "pad",
484 "f_window",
485 "bandpass",
486 ]
487 pdict = {
488 "pad": self.zero_pad,
489 "zero_mean": self.zero_mean,
490 "detrend": self.detrend,
491 "t_window": self.t_window,
492 "f_window": self.f_window,
493 "bandpass": self.bandpass,
494 }
495 subplot_dict = {}
496 count = 3
497 nrows = 2
498 for key in order:
499 value = pdict[key]
500 if value:
501 subplot_dict[key] = count
502 count += 2
503 nrows += 1
504 self.nrows = nrows
506 return subplot_dict
508 def remove_instrument_response(
509 self,
510 operation="divide",
511 include_decimation=None,
512 include_delay=None,
513 filters_to_remove=[],
514 ):
515 """
516 Remove instrument response following the recipe provided
518 :return: calibrated time series
519 :rtype: np.ndarray
521 """
522 # if filters to include not specified, get from self
523 if not filters_to_remove:
524 self.logger.debug("No explicit list of filters was passed to remove")
525 self.logger.debug("Will determine filters to remove ... ")
526 if include_decimation is None:
527 include_decimation = self.include_decimation
528 if include_delay is None:
529 include_delay = self.include_delay
530 filters_to_remove = self.channel_response.get_list_of_filters_to_remove(
531 include_decimation=include_decimation, include_delay=include_delay
532 )
533 if filters_to_remove is []:
534 raise ValueError("There are no filters in channel_response to remove")
536 ts = np.copy(self.ts)
537 f = np.fft.rfftfreq(ts.size, d=self.sample_interval)
538 step = 1
539 if self.plot:
540 self.subplot_dict = self._get_subplot_count()
541 self.fig = plt.figure(figsize=[10, 12])
542 self.fig.clf()
543 ax = self.fig.add_subplot(self.nrows, 2, 1)
544 ax2 = self.fig.add_subplot(self.nrows, 2, 2)
545 (l1,) = ax.plot(self.time_array, ts, color="k", lw=2)
546 ax.set_xlim((self.time_array[0], self.time_array[-1]))
547 (l2,) = ax2.loglog(f, abs(np.fft.rfft(ts)), "k", lw=2)
548 ax2.set_xlim((f[0], f[-1]))
549 ax.legend(
550 [l1],
551 ["Original"],
552 loc="upper left",
553 borderaxespad=0.01,
554 borderpad=0.1,
555 )
556 # detrend
557 if self.detrend:
558 ts = self.apply_detrend(ts)
559 self.logger.debug(f"Step {step}: Applying Linear Detrend")
560 step += 1
561 # zero mean
562 if self.zero_mean:
563 ts = self.apply_zero_mean(ts)
564 self.logger.debug(f"Step {step}: Removing Mean")
565 step += 1
566 # filter in time domain
567 if self.t_window is not None:
568 ts = self.apply_t_window(ts)
569 self.logger.debug(f"Step {step}: Applying {self.t_window} Time Window")
570 step += 1
571 if self.plot:
572 wax = self.fig.get_axes()[self.subplot_dict["t_window"] - 1].twinx()
573 (tw,) = wax.plot(
574 self.time_array,
575 self.get_window(self.t_window, self.t_window_params, ts.size),
576 color=(0.75, 0.75, 0.75),
577 zorder=0,
578 )
579 if self.zero_pad:
580 # pad the time series to a power of 2, this may be overkill, especially for long time series
581 ts = self.apply_zero_pad(ts)
582 self.logger.debug(f"Step {step}: Applying Zero Padding")
583 step += 1
584 # get the real frequencies of the FFT -- zero pad may have changed ts.size
585 f = np.fft.rfftfreq(ts.size, d=self.sample_interval)
587 # compute the complex response given the frequency range of the FFT
588 # the complex response assumes frequencies are in reverse order and flip them on input
589 # so we need to flip the complex reponse so it aligns with the fft.
590 cr = self.channel_response.complex_response(
591 f,
592 filters_list=filters_to_remove,
593 )[::-1]
594 # remove the DC term at frequency == 0
595 cr[-1] = abs(cr[-2]) + 0.0j
597 data = np.fft.rfft(ts)
598 # remove DC term
599 data[-1] = abs(data[-1]) + 0.0j
601 # if a window is requested then create it here and mulitply by the data
602 # the windows are designed to be symmetrical about frequency = 0
603 # here we are taking only the real part of the FFT so we cut the window in half
604 if self.f_window is not None:
605 data = self.apply_f_window(data)
606 self.logger.debug(f"Step {step}: Applying {self.f_window} Frequency Window")
607 step += 1
609 # calibrate the time series, compute real part of fft, divide out
610 # channel response, inverse fft
611 if operation == "divide":
612 calibrated_ts = np.fft.irfft(data / cr)[0 : self.ts.size]
613 self.logger.debug(
614 f"Step {step}: Removing Calibration via divide channel response"
615 )
616 step += 1
617 elif operation == "multiply":
618 calibrated_ts = np.fft.irfft(data * cr)[0 : self.ts.size]
619 self.logger.warning(
620 f"Instrument response being applied rather that expected "
621 f"operation of removing the response"
622 )
623 step += 1
624 else:
625 msg = f"Operation {operation} not recognized method of instrument response correction"
626 logger.error(msg)
627 raise Exception
629 # If a time window was applied, need to un-apply it to reconstruct the signal.
630 if self.t_window is not None:
631 w = self.get_window(self.t_window, self.t_window_params, calibrated_ts.size)
632 with np.errstate(divide="ignore", invalid="ignore"):
633 calibrated_ts = calibrated_ts / w
634 self.logger.debug(f"Step {step}: Un-applying Time Window")
635 step += 1
636 if self.bandpass:
637 calibrated_ts = self.apply_bandpass(calibrated_ts)
638 self.logger.debug(f"Step {step}: Applying Bandpass Filter")
639 step += 1
640 if self.plot:
641 self._subplots(
642 self.time_array[0 : calibrated_ts.size],
643 calibrated_ts,
644 (1, 0.1, 0.1),
645 self.nrows * 2 - 1,
646 "Calibrated",
647 )
648 self.fig.get_axes()[-2].set_ylabel(self.channel_response.units_in)
649 if self.t_window is not None:
650 wax = self.fig.get_axes()[-2].twinx()
651 (tw,) = wax.plot(
652 self.time_array, 1.0 / w, color=(0.5, 0.5, 0.5), zorder=0
653 )
654 self.fig.tight_layout()
655 plt.show()
656 return calibrated_ts
659def adaptive_notch_filter(
660 bx,
661 df=100,
662 notches=[50, 100],
663 notchradius=0.5,
664 freqrad=0.9,
665 rp=0.1,
666 dbstop_limit=5.0,
667):
668 """
670 :param bx: time series to filter
671 :type bx: np.ndarray
672 :param df: sample rate in samples per second, defaults to 100
673 :type df: float, optional
674 :param notches: list of frequencies to locate notches at in Hz,
675 defaults to [50, 100]
676 :type notches: list, optional
677 :param notchradius: notch radius, defaults to 0.5
678 :type notchradius: float, optional
679 :param freqrad: radius to search for a peak at the notch frequency,
680 defaults to 0.9
681 :type freqrad: float, optional
682 :param rp: ripple of Chebyshev type 1 filter, lower numbers means less
683 ripples, defaults to 0.1
684 :type rp: float, optional
685 :param dbstop_limit: limits the difference between the peak at the
686 notch and surrounding spectra. Any difference
687 above dbstop_limit will be filtered, anything
688 less will not, defaults to 5.0
689 :type dbstop_limit: float, optional
690 :return: notch filtered data
691 :rtype: np.ndarray
692 :return: list of notch frequencies
693 :rtype: list
696 Example
697 ---------
699 >>> import RemovePeriodicNoise_Kate as rmp
700 >>> # make a variable for the file to load in
701 >>> fn = r"/home/MT/mt01_20130101_000000.BX"
702 >>> # load in file, if the time series is not an ascii file
703 >>> # might need to add keywords to np.loadtxt or use another
704 >>> # method to read in the file
705 >>> bx = np.loadtxt(fn)
706 >>> # create a list of frequencies to filter out
707 >>> freq_notches = [50, 150, 200]
708 >>> # filter data
709 >>> bx_filt, filt_lst = rmp.adaptiveNotchFilter(bx, df=100.
710 >>> ... notches=freq_notches)
711 >>> #save the filtered data into a file
712 >>> np.savetxt(r"/home/MT/Filtered/mt01_20130101_000000.BX", bx_filt)
714 Notes:
716 Most of the time the default parameters work well, the only thing
717 you need to change is the notches and perhaps the radius. I would
718 test it out with a few time series to find the optimum parameters.
719 Then make a loop over all you time series data. Something like
721 >>> import os
722 >>> dirpath = r"/home/MT"
723 >>> #make a director to save filtered time series
724 >>> save_path = r"/home/MT/Filtered"
725 >>> if not os.path.exists(save_path):
726 >>> os.mkdir(save_path)
727 >>> for fn in os.listdir(dirpath):
728 >>> bx = np.loadtxt(os.path.join(dirpath, fn)
729 >>> bx_filt, filt_lst = rmp.adaptiveNotchFilter(bx, df=100.
730 >>> ... notches=freq_notches)
731 >>> np.savetxt(os.path.join(save_path, fn), bx_filt)
733 """
735 bx = np.array(bx)
737 if type(notches) is list:
738 notches = np.array(notches)
739 elif type(notches) in [float, int]:
740 notches = np.array([notches], dtype=np.float)
741 df = float(df) # make sure df is a float
742 dt = 1.0 / df # sampling rate
744 # transform data into frequency domain to find notches
745 BX = np.fft.fft(zero_pad(bx))
746 n = len(BX) # length of array
747 dfn = df / n # frequency step
748 dfnn = int(freqrad / dfn) # radius of frequency search
749 fn = notchradius # filter radius
750 freq = np.fft.fftfreq(n, dt)
752 filtlst = []
753 for notch in notches:
754 if notch > freq.max():
755 break
756 else:
757 fspot = int(round(notch / dfn))
758 nspot = np.where(
759 abs(BX) == max(abs(BX[max([fspot - dfnn, 0]) : min([fspot + dfnn, n])]))
760 )[0][0]
762 med_bx = np.median(
763 abs(BX[max([nspot - dfnn * 10, 0]) : min([nspot + dfnn * 10, n])]) ** 2
764 )
766 # calculate difference between peak and surrounding spectra in dB
767 dbstop = 10 * np.log10(abs(BX[nspot]) ** 2 / med_bx)
768 if np.nan_to_num(dbstop) == 0.0 or dbstop < dbstop_limit:
769 filtlst.append("No need to filter \n")
770 else:
771 filtlst.append([freq[nspot], dbstop])
772 ws = 2 * np.array([freq[nspot] - fn, freq[nspot] + fn]) / df
773 wp = 2 * np.array([freq[nspot] - 2 * fn, freq[nspot] + 2 * fn]) / df
774 ford, wn = signal.cheb1ord(wp, ws, 1, dbstop)
775 b, a = signal.cheby1(1, 0.5, wn, btype="bandstop")
776 bx = signal.filtfilt(b, a, bx)
777 return bx, filtlst