[med-svn] [python-mne] 220/353: ENH: filter trans bw, warning
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:25:01 UTC 2015
This is an automated email from the git hooks/post-receive script.
yoh pushed a commit to tag 0.4
in repository python-mne.
commit c656aa82b1a95e3de04341dc20835d1d15270e4b
Author: mluessi at nmr.mgh.harvard.edu <mluessi at nmr.mgh.harvard.edu>
Date: Sat Jun 23 03:46:57 2012 -0400
ENH: filter trans bw, warning
---
mne/fiff/raw.py | 20 +++++++++------
mne/filter.py | 76 ++++++++++++++++++++++++++++++++++++++++++++++-----------
2 files changed, 74 insertions(+), 22 deletions(-)
diff --git a/mne/fiff/raw.py b/mne/fiff/raw.py
index e0f4e9a..b55068f 100644
--- a/mne/fiff/raw.py
+++ b/mne/fiff/raw.py
@@ -363,7 +363,7 @@ class Raw(object):
self.apply_function(hilbert, picks, np.complex64, n_jobs, verbose)
def filter(self, l_freq, h_freq, picks=None, filter_length=None,
- n_jobs=1, verbose=5):
+ trans_bw=0.5, n_jobs=1, verbose=5):
"""Filter a subset of channels.
Applies a zero-phase band-pass filter to the channels selected by
@@ -392,7 +392,8 @@ class Raw(object):
(n_times: number of timepoints in Raw object) the filter length
used is n_times. Otherwise, overlap-add filtering with a
filter of the specified length is used (faster for long signals).
-
+ trans_bw : float
+ Width of the transition band in Hz.
n_jobs: int (default: 1)
Number of jobs to run in parallel.
@@ -407,14 +408,17 @@ class Raw(object):
if picks is None:
picks = pick_types(self.info, meg=True, eeg=True)
if l_freq is None and h_freq is not None:
- self.apply_function(low_pass_filter, picks, None, n_jobs, verbose, fs,
- h_freq, filter_length=filter_length)
+ self.apply_function(low_pass_filter, picks, None, n_jobs, verbose,
+ fs, h_freq, filter_length=filter_length,
+ trans_bw=trans_bw)
if l_freq is not None and h_freq is None:
- self.apply_function(high_pass_filter, picks, None, n_jobs, verbose, fs,
- l_freq, filter_length=filter_length)
+ self.apply_function(high_pass_filter, picks, None, n_jobs, verbose,
+ fs, l_freq, filter_length=filter_length,
+ trans_bw=trans_bw)
if l_freq is not None and h_freq is not None:
- self.apply_function(band_pass_filter, picks, None, n_jobs, verbose, fs,
- l_freq, h_freq, filter_length=filter_length)
+ self.apply_function(band_pass_filter, picks, None, n_jobs, verbose,
+ fs, l_freq, h_freq,
+ filter_length=filter_length, trans_bw=trans_bw)
@deprecated('band_pass_filter is deprecated please use raw.filter instead')
def band_pass_filter(self, picks, l_freq, h_freq, filter_length=None,
diff --git a/mne/filter.py b/mne/filter.py
index 99eb8cc..0246fab 100644
--- a/mne/filter.py
+++ b/mne/filter.py
@@ -1,6 +1,7 @@
import warnings
import numpy as np
from scipy.fftpack import fft, ifft
+from scipy.signal import freqz
from .utils import firwin2 # back port for old scipy
@@ -134,6 +135,19 @@ def _overlap_add_filter(x, h, n_fft=None, zero_phase=True):
return x_filtered
+def _filter_attenuation(h, freq, gain):
+ """Compute minimum attenuation at stop frequency"""
+
+ _, filt_resp = freqz(h, worN=np.pi * freq)
+ filt_resp = np.abs(filt_resp) # use amplitude response
+ filt_resp /= np.max(filt_resp)
+ filt_resp[np.where(gain == 1)] = 0
+ idx = np.argmax(filt_resp)
+ att_db = -20 * np.log10(filt_resp[idx])
+ att_freq = freq[idx]
+
+ return att_db, att_freq
+
def _filter(x, Fs, freq, gain, filter_length=None):
"""Filter signal using gain control points in the frequency domain.
@@ -161,10 +175,15 @@ def _filter(x, Fs, freq, gain, filter_length=None):
xf : 1d array
x filtered
"""
+
+ # issue a warning if attenuation is less than this
+ min_att_db = 20
+
assert x.ndim == 1
# normalize frequencies
- freq = [f / (Fs / 2) for f in freq]
+ freq = np.array([f / (Fs / 2) for f in freq])
+ gain = np.array(gain)
if filter_length is None or len(x) <= filter_length:
# Use direct FFT filtering for short signals
@@ -180,6 +199,12 @@ def _filter(x, Fs, freq, gain, filter_length=None):
H = firwin2(N, freq, gain)
+ att_db, att_freq = _filter_attenuation(H, freq, gain)
+ if att_db < min_att_db:
+ att_freq *= Fs / 2
+ warnings.warn('Attenuation at stop frequency %0.1fHz is only '
+ '%0.1fdB.' % (att_freq, att_db))
+
# Make zero-phase filter function
B = np.abs(fft(H))
@@ -196,12 +221,21 @@ def _filter(x, Fs, freq, gain, filter_length=None):
N += 1
H = firwin2(N, freq, gain)
+
+ att_db, att_freq = _filter_attenuation(H, freq, gain)
+ att_db += 6 # the filter is applied twice (zero phase)
+ if att_db < min_att_db:
+ att_freq *= Fs / 2
+ warnings.warn('Attenuation at stop frequency %0.1fHz is only '
+ '%0.1fdB. Increase filter_length for higher '
+ 'attenuation.' % (att_freq, att_db))
+
xf = _overlap_add_filter(x, H, zero_phase=True)
return xf
-def band_pass_filter(x, Fs, Fp1, Fp2, filter_length=None):
+def band_pass_filter(x, Fs, Fp1, Fp2, filter_length=None, trans_bw=0.5):
"""Bandpass filter for the signal x.
Applies a zero-phase bandpass filter to the signal x.
@@ -220,6 +254,8 @@ def band_pass_filter(x, Fs, Fp1, Fp2, filter_length=None):
Length of the filter to use. If None or "len(x) < filter_length", the
filter length used is len(x). Otherwise, overlap-add filtering with a
filter of the specified length is used (faster for long signals).
+ trans_bw : float
+ Width of the transition band in Hz.
Returns
-------
@@ -238,17 +274,21 @@ def band_pass_filter(x, Fs, Fp1, Fp2, filter_length=None):
| |
Fs1 Fp1 Fp2 Fs2
- DEFAULTS values
- Fs1 = Fp1 - 0.5 in Hz
- Fs2 = Fp2 + 0.5 in Hz
+ Where
+ Fs1 = Fp1 - trans_bw in Hz
+ Fs2 = Fp2 + trans_bw in Hz
"""
Fs = float(Fs)
Fp1 = float(Fp1)
Fp2 = float(Fp2)
- # Default values in Hz
- Fs1 = max(Fp1 - 0.5, 0.)
- Fs2 = Fp2 + 0.5
+ Fs1 = Fp1 - trans_bw
+ Fs2 = Fp2 + trans_bw
+
+ if Fs1 <= 0:
+ raise ValueError('Filter specification invalid: Lower stop frequency '
+ 'too low (%0.1fHz). Increase Fp1 or reduce '
+ 'transition bandwidth (trans_bw)' % Fs1)
xf = _filter(x, Fs, [0, Fs1, Fp1, Fp2, Fs2, Fs / 2], [0, 0, 1, 1, 0, 0],
filter_length)
@@ -256,7 +296,7 @@ def band_pass_filter(x, Fs, Fp1, Fp2, filter_length=None):
return xf
-def low_pass_filter(x, Fs, Fp, filter_length=None):
+def low_pass_filter(x, Fs, Fp, filter_length=None, trans_bw=0.5):
"""Lowpass filter for the signal x.
Applies a zero-phase lowpass filter to the signal x.
@@ -273,6 +313,8 @@ def low_pass_filter(x, Fs, Fp, filter_length=None):
Length of the filter to use. If None or "len(x) < filter_length", the
filter length used is len(x). Otherwise, overlap-add filtering with a
filter of the specified length is used (faster for long signals).
+ trans_bw : float
+ Width of the transition band in Hz.
Returns
-------
@@ -289,20 +331,20 @@ def low_pass_filter(x, Fs, Fp, filter_length=None):
| \
| -----------------
|
- Fp Fp+0.5
+ Fp Fp+trans_bw
"""
Fs = float(Fs)
Fp = float(Fp)
- Fstop = Fp + 0.5
+ Fstop = Fp + trans_bw
xf = _filter(x, Fs, [0, Fp, Fstop, Fs / 2], [1, 1, 0, 0], filter_length)
return xf
-def high_pass_filter(x, Fs, Fp, filter_length=None):
+def high_pass_filter(x, Fs, Fp, filter_length=None, trans_bw=0.5):
"""Highpass filter for the signal x.
Applies a zero-phase highpass filter to the signal x.
@@ -319,6 +361,8 @@ def high_pass_filter(x, Fs, Fp, filter_length=None):
Length of the filter to use. If None or "len(x) < filter_length", the
filter length used is len(x). Otherwise, overlap-add filtering with a
filter of the specified length is used (faster for long signals).
+ trans_bw : float
+ Width of the transition band in Hz.
Returns
-------
@@ -335,14 +379,18 @@ def high_pass_filter(x, Fs, Fp, filter_length=None):
/ |
---------- |
|
- Fp-0.5 Fp
+ Fp-trans_bw Fp
"""
Fs = float(Fs)
Fp = float(Fp)
- Fstop = Fp - 0.5
+ Fstop = Fp - trans_bw
+ if Fstop <= 0:
+ raise ValueError('Filter specification invalid: Stop frequency too low'
+ '(%0.1fHz). Increase Fp or reduce transition '
+ 'bandwidth (trans_bw)' % Fstop)
xf = _filter(x, Fs, [0, Fstop, Fp, Fs / 2], [0, 0, 1, 1], filter_length)
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/python-mne.git
More information about the debian-med-commit
mailing list