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: [email protected] <[email protected]> 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 _______________________________________________ debian-med-commit mailing list [email protected] http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/debian-med-commit
