Changeset 13330 for fact/tools/pyscripts/pyfact/fir_filter.py

Timestamp:

04/11/12 08:55:46 (13 years ago)

Author:

neise

Message:

debugged filter artifacts in the first few slices, by hand.... ugly but works

File:

• fact/tools/pyscripts/pyfact/fir_filter.py (modified) (2 diffs)

fact/tools/pyscripts/pyfact/fir_filter.py

@@ -30 +30 @@
             
             *data* 1D or 2D numpy array
-        """
-        length = max(len(self.a),len(self.b))-1
-        if length > 0:
-            if ( data.ndim == 1):
-                initial = np.ones(length)
-                initial *= data[0]
-            elif ( data.ndim == 2):
-                initial = np.ones( (data.shape[0], length) )
-                for i in range(data.shape[0]):
-                    initial[i,:] *= data[i,0]
-            else:
-                print 'HELP.'
-                pass
             
-            filtered, zf = signal.lfilter(self.b, self.a, data, zi=initial)
+            remark:
+            I did not understand how to use the initial filter conditions of lfilter()
+            to produce the output, I expected.
+            So I apply the filters as follows. 
+            the filter *delay* is equal to its length-1
+            Then I extend the input data by this delay-length, adding copies of the 
+            first value.
+            Then the filter runs ovter this extended data.
+            The output will have a filter artifact in the first samples, which 
+            will be cut off anyway, because they were artificially added before.
+        """
+        delay = max(len(self.a),len(self.b))-1
+        
+        if ( data.ndim == 1):
+            initial = np.ones(delay)
+            initial *= data[0]
+        elif ( data.ndim == 2):
+            initial = np.ones( (data.shape[0], delay) )
+            for i in range(data.shape[0]):
+                initial[i,:] *= data[i,0]
         else:
-            filtered= signal.lfilter(self.b, self.a, data)            
-        filtered = filtered.reshape(data.shape)
+            print 'HELP.'
+            pass
+        data = np.hstack( (initial,data) )
+
+        filtered= signal.lfilter(self.b, self.a, data)
+        if ( data.ndim == 1):
+            filtered = filtered[delay:]
+        elif ( data.ndim == 2):
+            filtered = filtered[:,delay:]
+
         return filtered
 
@@ -123 +137 @@
     
 
+def _test_SlidingAverage2():
+    """ test the sliding average function
+    use a step function as input
+    """
+    from plotters import Plotter
+    from generator import SignalGenerator
+    generate = SignalGenerator()
+    plot = Plotter('_test_SlidingAverage')
+
+    safilter = SlidingAverage(8)
+    print safilter
+
+    signal = np.ones( (6,20) ) * 3.0
+    filtered = safilter(signal)
+    #plot( [signal[0], filtered[0]] , ['original', 'filtered'] )
+
+    raw_input('press any key to go on')
+    plt.close(plot.figure)
+
+
+
 def _test_CFD():
     """ test the remove signal function