[med-svn] [python-mne] 02/353: zero phase, not yet working correctly

[Yaroslav Halchenko]

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 e3a2b8a16293c902ed83d5f26765fb8bdddddac3
Author: Martin Luessi <mluessi at nmr.mgh.harvard.edu>
Date:   Fri Oct 14 09:13:00 2011 -0400

    zero phase, not yet working correctly
---
 mne/filter.py | 91 ++++++++++++++++++++++++++++++++---------------------------
 1 file changed, 50 insertions(+), 41 deletions(-)

diff --git a/mne/filter.py b/mne/filter.py
index b492e3c..90a88e0 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):
+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.
@@ -25,9 +25,11 @@ def overlap_add_filter(x, h, N_fft=None):
     
     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_x)))
-            cost = np.ceil(N_x / (N - N_h + 1)) * N  * (np.log2(N) + 1)                
+                             np.floor(np.log2(N_tot)))
+            cost = np.ceil(N_tot / (N - N_h + 1)) * N  * (np.log2(N) + 1)                
             N_fft = N[np.argmin(cost)]
         else:
             # Use only a single block
@@ -41,32 +43,53 @@ def overlap_add_filter(x, h, N_fft=None):
         raise ValueError('N_fft is too short, has to be at least len(h)')
 
     # Filter in frequency domain    
-    # FIXME: abs?
-    h_fft = np.abs(fft(np.concatenate((h, np.zeros(N_fft - N_h)))))   
-
+    h_fft = fft(np.concatenate((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)
 
-    # Go through segments
+    # Number of segements (including fractional segments)
     num_segments = int(np.ceil(N_x / float(N_seg)))
+    
+    filter_input = x
 
-    for seg_ind in range(num_segments):
-        seg = x[seg_ind*N_seg:(seg_ind+1)*N_seg]
-        seg_fft = fft(np.concatenate((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] \
-                += np.real(ifft(h_fft * seg_fft))            
-        else:
-            # Last segment
-            x_filtered[seg_ind*N_seg:] \
-                += np.real(ifft(h_fft * seg_fft))[:N_x-seg_ind*N_seg]
- 
+    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)
+            
+        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)))))
+    
+    
+            if seg_ind*N_seg+N_fft < N_x:                               
+                x_filtered[seg_ind*N_seg:seg_ind*N_seg+N_fft] \
+                    += np.real(ifft(h_fft * seg_fft))            
+            else:
+                # Last segment
+                x_filtered[seg_ind*N_seg:] \
+                    += np.real(ifft(h_fft * seg_fft))[:N_x-seg_ind*N_seg]
+    if zero_phase:
+        # flip signal back
+        x_filtered = x_filtered[::-1]
+        
     return x_filtered           
         
     
 
-def band_pass_filter(x, Fs, Fp1, Fp2, ov_add):
+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
@@ -127,28 +150,14 @@ def band_pass_filter(x, Fs, Fp1, Fp2, ov_add):
     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))
-#        
-    
-    if ov_add:
-        N = int(2 * Fs)
-        B = signal.firwin2(N, [0, Ns1, Np1, Np2, Ns2, 1], [0, 0, 1, 1, 0, 0])
+    N = len(x)
 
-        xf = overlap_add_filter(x, B)    
-    else:
-        N = int(2 * Fs)
-        B = signal.firwin2(N, [0, Ns1, Np1, Np2, Ns2, 1], [0, 0, 1, 1, 0, 0])
+    B = signal.firwin2(N, [0, Ns1, Np1, Np2, Ns2, 1], [0, 0, 1, 1, 0, 0])
 
-        B = ifft(np.abs(fft(B)))
-        # Make zero-phase filter function
-        
-        xf = (signal.convolve(x, B, 'full'))
-        
+    # Make zero-phase filter function
+    H = np.abs(fft(B))
+
+    xf = np.real(ifft(fft(x) * H))
     xf = xf[:Norig]
     x = x[:Norig]
 

-- 
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/python-mne.git