[med-svn] [python-mne] 03/353: zero-phase using overlap-add working, code needs cleanup
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 5e2766fd27935ad4a10a5eef9d241245cfdbdbc9
Author: Martin Luessi <mluessi at nmr.mgh.harvard.edu>
Date: Fri Oct 14 17:06:46 2011 -0400
zero-phase using overlap-add working, code needs cleanup
---
mne/filter.py | 149 ++++++++++++++++++++++++++++++++++++++++++++++++----------
1 file changed, 123 insertions(+), 26 deletions(-)
diff --git a/mne/filter.py b/mne/filter.py
index 90a88e0..7f6f171 100644
--- a/mne/filter.py
+++ b/mne/filter.py
@@ -3,7 +3,7 @@ from scipy import signal
from scipy.fftpack import fft, ifft
-def overlap_add_filter(x, h, N_fft=None, zero_phase=True):
+def _overlap_add_filter(x, h, N_fft=None, zero_phase=True):
""" Calculate linear convolution using overlap-add FFTs
Implements the linear convolution between x and h using overlap-add FFTs.
@@ -17,16 +17,24 @@ def overlap_add_filter(x, h, N_fft=None, zero_phase=True):
N_fft : int
Length of the FFT. If None, the best size is determined automatically.
"""
-
- # https://ccrma.stanford.edu/~jos/fp/Forward_Backward_Filtering.html
- # http://www.mathworks.com/matlabcentral/fileexchange/17061-filtfilthd/content/filtfilthd.m
+
N_h = len(h)
- N_x = len(x)
+
+ # Extend the signal at the edges (see scipy.signal.filtfitlt)
+ #if N_h < 3 * len(x):
+ 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)
+
if N_fft == None:
if N_x > N_h:
N_tot = 2*N_x if zero_phase else N_x
-
+
N = 2**np.arange(np.ceil(np.log2(N_h)),
np.floor(np.log2(N_tot)))
cost = np.ceil(N_tot / (N - N_h + 1)) * N * (np.log2(N) + 1)
@@ -34,45 +42,43 @@ def overlap_add_filter(x, h, N_fft=None, zero_phase=True):
else:
# 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)')
# Filter in frequency domain
- h_fft = fft(np.concatenate((h, np.zeros(N_fft - N_h))))
-
+ h_fft = fft(np.r_[h, np.zeros(N_fft - N_h)])
+
if zero_phase:
# We will apply the filter in forward and backward direction: Scale
# frequency response of the filter so that the shape of the amplitude
# response stays the same when it is applied twice
# be careful not to divide by too small numbers
- #idx = np.where(np.abs(h_fft) > 1e-6)
- #h_fft[idx] = h_fft[idx] / np.sqrt(np.abs(h_fft[idx]))
- pass
-
- x_filtered = np.zeros_like(x)
+ idx = np.where(np.abs(h_fft) > 1e-6)
+ h_fft[idx] = h_fft[idx] / np.sqrt(np.abs(h_fft[idx]))
+
+ x_filtered = np.zeros_like(x_ext)
# Number of segements (including fractional segments)
num_segments = int(np.ceil(N_x / float(N_seg)))
-
- filter_input = x
+
+ filter_input = x_ext
for pass_no in range(2) if zero_phase else range(1):
-
+
if pass_no == 1:
# second pass: flip signal
- filter_input = x_filtered[::-1]
- x_filtered = np.zeros_like(x)
-
+ filter_input = np.flipud(x_filtered)
+ x_filtered = np.zeros_like(x_ext)
+
for seg_ind in range(num_segments):
seg = filter_input[seg_ind*N_seg:(seg_ind+1)*N_seg]
- seg_fft = fft(np.concatenate((seg, np.zeros(N_fft - len(seg)))))
-
+ seg_fft = fft(np.r_[seg, np.zeros(N_fft - len(seg))])
if seg_ind*N_seg+N_fft < N_x:
x_filtered[seg_ind*N_seg:seg_ind*N_seg+N_fft] \
@@ -81,13 +87,104 @@ def overlap_add_filter(x, h, N_fft=None, zero_phase=True):
# Last segment
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]
+
if zero_phase:
# flip signal back
- x_filtered = x_filtered[::-1]
-
+ x_filtered = np.flipud(x_filtered)
+
return x_filtered
-
+
+def _filter(x, Fs, freq, gain):
+
+ 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
+
+ # Make x EVEN
+ Norig = len(x)
+ if Norig % 2 == 1:
+ x = np.r_[x, 1]
+ N = len(x)
+
+ H = signal.firwin2(N, freq, gain)
+
+ # Make zero-phase filter function
+ B = np.abs(fft(H))
+
+ xf = np.real(ifft(fft(x) * B))
+ xf = xf[:Norig]
+ x = x[:Norig]
+ else:
+ # Use overlap-add filter
+ N = int(10 * Fs)
+ H = signal.firwin2(N, freq, gain)
+ xf = _overlap_add_filter(x, H, zero_phase=True)
+
+ return xf
+
+
+def new_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
+ """
+ Fs = float(Fs)
+ Fp1 = float(Fp1)
+ Fp2 = float(Fp2)
+
+ # Default values in Hz
+ Fs1 = Fp1 - 0.5
+ Fs2 = Fp2 + 0.5
+
+ xf = _filter(x, Fs, [0, Fs1, Fp1, Fp2, Fs2, Fs/2], [0, 0, 1, 1, 0, 0])
+
+ return xf
+
def band_pass_filter(x, Fs, Fp1, Fp2):
"""Bandpass filter for the signal x.
--
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