[med-svn] [python-mne] 04/353: All filters using overlap-add
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:24:23 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 02baee5e69fa838e5617b7982768e5322a310c1b
Author: Martin Luessi <mluessi at nmr.mgh.harvard.edu>
Date: Tue Oct 25 17:23:08 2011 -0400
All filters using overlap-add
---
mne/filter.py | 222 ++++++++++++++++++----------------------------------------
1 file changed, 67 insertions(+), 155 deletions(-)
diff --git a/mne/filter.py b/mne/filter.py
index 7f6f171..c37e18d 100644
--- a/mne/filter.py
+++ b/mne/filter.py
@@ -4,9 +4,12 @@ from scipy.fftpack import fft, ifft
def _overlap_add_filter(x, h, N_fft=None, zero_phase=True):
- """ Calculate linear convolution using overlap-add FFTs
+ """ Filter using overlap-add FFTs.
- Implements the linear convolution between x and h using overlap-add FFTs.
+ Filters the signal x using a filter with the impulse response h.
+ If zero_phase==True, the amplitude response is scaled and the filter is
+ applied in forward and backward direction, resulting in a zero-phase
+ filter.
Parameters
----------
@@ -16,21 +19,27 @@ def _overlap_add_filter(x, h, N_fft=None, zero_phase=True):
Filter impule response (FIR filter coefficients)
N_fft : int
Length of the FFT. If None, the best size is determined automatically.
+ zero_phase : boolean
+ If True: the filter is applied in forward and backward direction,
+ resulting in a zero-phase filter
+
+ Returns
+ -------
+ xf : 1d array
+ x filtered
"""
N_h = len(h)
+
+ # Extend the signal by mirroring the edges to reduce transient filter
+ # response
+ N_edge = min(N_h, len(x))
- # Extend the signal at the edges (see scipy.signal.filtfitlt)
- #if N_h < 3 * len(x):
+ x_ext = np.r_[2*x[0]-x[N_edge-1:0:-1], x, 2*x[-1]-x[-2:-N_edge-1:-1]]
- edge = 3 * N_h
-
- x_ext = np.r_[2*x[0]-x[edge-1:0:-1], x, 2*x[-1]-x[-2:-edge-1:-1]]
-
- print 'length x_ext: %d edge: %d' % (len(x_ext), edge)
-
N_x = len(x_ext)
+ # Determine FFT length to use
if N_fft == None:
if N_x > N_h:
N_tot = 2*N_x if zero_phase else N_x
@@ -43,10 +52,6 @@ def _overlap_add_filter(x, h, N_fft=None, zero_phase=True):
# Use only a single block
N_fft = 2**np.ceil(np.log2(N_x + N_h - 1))
- print 'FFT length: %d' % (N_fft)
-
- N_seg = N_fft - N_h + 1
-
if N_fft <= 0:
raise ValueError('N_fft is too short, has to be at least len(h)')
@@ -64,8 +69,11 @@ def _overlap_add_filter(x, h, N_fft=None, zero_phase=True):
x_filtered = np.zeros_like(x_ext)
+ # Segment length for signal x
+ N_seg = N_fft - N_h + 1
+
# Number of segements (including fractional segments)
- num_segments = int(np.ceil(N_x / float(N_seg)))
+ num_segments = np.ceil(N_x / float(N_seg))
filter_input = x_ext
@@ -88,8 +96,8 @@ def _overlap_add_filter(x, h, N_fft=None, zero_phase=True):
x_filtered[seg_ind*N_seg:] \
+= np.real(ifft(h_fft * seg_fft))[:N_x-seg_ind*N_seg]
- # Remove edges that we added
- x_filtered = x_filtered[edge-1:-edge+1]
+ # Remove mirrored edges that we added
+ x_filtered = x_filtered[N_edge-1:-N_edge+1]
if zero_phase:
# flip signal back
@@ -99,16 +107,39 @@ def _overlap_add_filter(x, h, N_fft=None, zero_phase=True):
def _filter(x, Fs, freq, gain):
+ """ Filter signal using gain control points in the frequency domain.
+
+ The filter impulse response is constructed from a Hamming window (window
+ used in "firwin2" function) to avoid ripples in the frequency reponse
+ (windowing is a smoothing in frequency domain).
+ If the signal is shorter than 10 seconds, it is filtered directly using
+ FFTs (zero-phase). If the signal is longer, zero-phase overlap-add
+ filtering is used.
+
+ Parameters
+ ----------
+ x : 1d array
+ Signal to filter
+ Fs : float
+ Sampling rate in Hz
+ freq : 1d array
+ Frequency sampling points in Hz
+ gain : 1d array
+ Filter gain at frequency sampling points
+
+ Returns
+ -------
+ xf : 1d array
+ x filtered
+ """
assert x.ndim == 1
# normalize frequencies
freq = [f / (Fs / 2) for f in freq]
if len(x) < 10*Fs:
- # TODO: how to decide which filter to use
-
- # Short signal: use direct FFT filter
+ # Use direct FFT filtering for short signals
# Make x EVEN
Norig = len(x)
@@ -133,14 +164,11 @@ def _filter(x, Fs, freq, gain):
return xf
-def new_band_pass_filter(x, Fs, Fp1, Fp2):
+def band_pass_filter(x, Fs, Fp1, Fp2):
"""Bandpass filter for the signal x.
- An acausal fft algorithm is applied (i.e. no phase shift). The filter
- functions is constructed from a Hamming window (window used in "firwin2"
- function) to avoid ripples in the frequency reponse (windowing is a
- smoothing in frequency domain)
-
+ Applies a zero-phase bandpass filter to the signal x.
+
Parameters
----------
x : 1d array
@@ -186,88 +214,10 @@ def new_band_pass_filter(x, Fs, Fp1, Fp2):
return xf
-def band_pass_filter(x, Fs, Fp1, Fp2):
- """Bandpass filter for the signal x.
-
- An acausal fft algorithm is applied (i.e. no phase shift). The filter
- functions is constructed from a Hamming window (window used in "firwin2"
- function) to avoid ripples in the frequency reponse (windowing is a
- smoothing in frequency domain)
-
- Parameters
- ----------
- x : 1d array
- Signal to filter
- Fs : float
- sampling rate
- Fp1 : float
- low cut-off frequency
- Fp2 : float
- high cut-off frequency
-
- Returns
- -------
- xf : array
- x filtered
-
- Notes
- -----
- The passbands (Fp1 Fp2) frequencies are defined in Hz as
- ----------
- /| | \
- / | | \
- / | | \
- / | | \
- ---------- | | -----------------
- | |
- Fs1 Fp1 Fp2 Fs2
-
- DEFAULTS values
- Fs1 = Fp1 - 0.5 in Hz
- Fs2 = Fp2 + 0.5 in Hz
- """
- Fp1 = float(Fp1)
- Fp2 = float(Fp2)
-
- # Default values in Hz
- Fs1 = Fp1 - 0.5
- Fs2 = Fp2 + 0.5
-
- assert x.ndim == 1
-
- # Make x EVEN
- Norig = len(x)
- if Norig % 2 == 1:
- x = np.r_[x, 1]
-
- # Normalize frequencies
- Ns1 = Fs1 / (Fs / 2)
- Ns2 = Fs2 / (Fs / 2)
- Np1 = Fp1 / (Fs / 2)
- Np2 = Fp2 / (Fs / 2)
-
- # Construct the filter function H(f)
- N = len(x)
-
- B = signal.firwin2(N, [0, Ns1, Np1, Np2, Ns2, 1], [0, 0, 1, 1, 0, 0])
-
- # Make zero-phase filter function
- H = np.abs(fft(B))
-
- xf = np.real(ifft(fft(x) * H))
- xf = xf[:Norig]
- x = x[:Norig]
-
- return xf
-
-
def low_pass_filter(x, Fs, Fp):
"""Lowpass filter for the signal x.
- An acausal fft algorithm is applied (i.e. no phase shift). The filter
- functions is constructed from a Hamming window (window used in "firwin2"
- function) to avoid ripples in the frequency reponse (windowing is a
- smoothing in frequency domain)
+ Applies a zero-phase lowpass filter to the signal x.
Parameters
----------
@@ -296,41 +246,20 @@ def low_pass_filter(x, Fs, Fp):
Fp Fp+0.5
"""
+ Fs = float(Fs)
Fp = float(Fp)
-
- assert x.ndim == 1
-
- # Make x EVEN
- Norig = len(x)
- if Norig % 2 == 1:
- x = np.r_[x, 1]
-
- # Normalize frequencies
- Ns = (Fp + 0.5) / (Fs / 2)
- Np = Fp / (Fs / 2)
-
- # Construct the filter function H(f)
- N = len(x)
-
- B = signal.firwin2(N, [0, Np, Ns, 1], [1, 1, 0, 0])
-
- # Make zero-phase filter function
- H = np.abs(fft(B))
-
- xf = np.real(ifft(fft(x) * H))
- xf = xf[:Norig]
- x = x[:Norig]
-
+
+ Fstop = Fp + 0.5
+
+ xf = _filter(x, Fs, [0, Fp, Fstop, Fs/2], [1, 1, 0, 0])
+
return xf
def high_pass_filter(x, Fs, Fp):
"""Highpass filter for the signal x.
- An acausal fft algorithm is applied (i.e. no phase shift). The filter
- functions is constructed from a Hamming window (window used in "firwin2"
- function) to avoid ripples in the frequency reponse (windowing is a
- smoothing in frequency domain)
+ Applies a zero-phase highpass filter to the signal x.
Parameters
----------
@@ -359,29 +288,12 @@ def high_pass_filter(x, Fs, Fp):
Fp-0.5 Fp
"""
+
+ Fs = float(Fs)
Fp = float(Fp)
+
+ Fstop = Fp - 0.5
- assert x.ndim == 1
-
- # Make x ODD
- Norig = len(x)
- if Norig % 2 == 0:
- x = np.r_[x, 1]
-
- # Normalize frequencies
- Ns = (Fp - 0.5) / (Fs / 2)
- Np = Fp / (Fs / 2)
-
- # Construct the filter function H(f)
- N = len(x)
-
- B = signal.firwin2(N, [0, Ns, Np, 1], [0, 0, 1, 1])
-
- # Make zero-phase filter function
- H = np.abs(fft(B))
-
- xf = np.real(ifft(fft(x) * H))
- xf = xf[:Norig]
- x = x[:Norig]
+ xf = _filter(x, Fs, [0, Fstop, Fp, Fs/2], [0, 0, 1, 1])
return xf
--
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