Source Code for Module qubx.fast.filter

  1  """Compiled routines for filtering. 
  2   
  3  Copyright 2008-2011 Research Foundation State University of New York  
  4  This file is part of QUB Express.                                           
  5   
  6  QUB Express is free software; you can redistribute it and/or modify           
  7  it under the terms of the GNU General Public License as published by  
  8  the Free Software Foundation, either version 3 of the License, or     
  9  (at your option) any later version.                                   
 10   
 11  QUB Express is distributed in the hope that it will be useful,                
 12  but WITHOUT ANY WARRANTY; without even the implied warranty of        
 13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         
 14  GNU General Public License for more details.                          
 15   
 16  You should have received a copy of the GNU General Public License,    
 17  named LICENSE.txt, in the QUB Express program directory.  If not, see         
 18  <http://www.gnu.org/licenses/>.                                       
 19   
 20  """ 
 21   
 22  from qubx.fast.fast_utils import * 
 23   
 24  qubfast.qub_filterdatad.argtypes = (c_double_p, c_double_p, c_int, c_double, c_double) 
 25  qubfast.qub_filterdatad.restype = c_int 
 26  qubfast.qub_filterdataf.argtypes = (c_float_p, c_float_p, c_int, c_float, c_float) 
 27  qubfast.qub_filterdataf.restype = c_int 
 28   
 29  # slow!: 
 30  #import numpy 
 31  #import scipy.signal 
 32  #from math import * 
 33  #from qubx.types import memoize 
 34   
 35  #def make_kernel(fc): 
 36  #    sigma = 0.1325 / fc 
 37  #    if sigma < 0.6: 
 38  #        kernel = [0.0] * 3 
 39  #        kernel[0] = kernel[2] = sigma * sigma / 2.0 
 40  #        kernel[1] = 1.0 - kernel[0] 
 41  #        n = 1 
 42  #    else: 
 43  #        n = int(4 * sigma) 
 44  #        nfir = 2*n + 1 
 45  #        kernel = [0.0] * nfir 
 46  #        kernel[n] = 1.0 / (sqrt(2*pi) * sigma) 
 47   #       for k in xrange(1, n+1): 
 48  #            kernel[n+k] = kernel[n-k] = kernel[n] * exp(-((k/sigma)**2 / 2.0)) 
 49      #wrap: 
 50      #kernel[:n+1] = kernel[n:] 
 51      #kernel[n+1:] = reversed(kernel[1:n]) 
 52       
 53  #    return numpy.array(kernel) 
 54   
 55  #get_kernel = memoize(make_kernel) 
 56   
 57  #def ovrlap(samples, kernel): 
 58  #    ndata = samples.shape[0] 
 59  #    m = kernel.shape[0] 
 60  #    mpow2 = 4 
 61  #    while mpow2 < m: 
 62  #        mpow2 *= 2 
 63  #    nfft = 32 
 64  #    while nfft < m: 
 65  #        nfft *= 2 
 66  #    nd = nfft - (m - 1) 
 67  #    output = numpy.zeros(dtype=samples.dtype, shape=samples.shape) 
 68  #    piece = numpy.zeros(dtype=samples.dtype, shape=(nfft,)) 
 69  #    kpiece = numpy.zeros(dtype=samples.dtype, shape=(mpow2,)) 
 70  #    kpiece[:m] = kernel[:] 
 71  #    for np in xrange(0, ndata, nd): 
 72  #        piece[:] = 0.0 
 73  #        pstart = (m-1)/2 
 74  #        actual = min(nd, ndata-np) 
 75  #        piece[pstart:pstart+actual] = samples[np:np+actual] 
 76  #        piece = scipy.signal.fftconvolve(piece, kpiece, 'same') 
 77  #        nn = min(pstart-1, np-1) 
 78   #       output[np-nn-1:np] += piece[pstart-nn-1:pstart] 
 79  #        output[np:np+actual] += piece[pstart:pstart+actual] 
 80  #    return output 
 81   
 82 -def filter(samples, sampling_sec, freq_Hz): 
 83      """Returns a numpy.array of samples, low-pass filtered at freq_Hz.""" 
 84  #    kernel = get_kernel(sampling_sec * freq_Hz) 
 85  #    if kernel.shape[0] > 4095: 
 86  #        return ovrlap(samples, kernel) 
 87  #    output = numpy.convolve(samples, kernel, 'same') 
 88  #    return output.astype('float32') 
 89   
 90      try: 
 91          dtype = samples.dtype 
 92      except: 
 93          dtype = None 
 94      if dtype == numpy.float32: 
 95          output = numpy.zeros(shape=samples.shape, dtype=numpy.float32) 
 96          qubfast.qub_filterdataf(samples.ctypes.data_as(c_float_p), output.ctypes.data_as(c_float_p), 
 97                                  len(samples), sampling_sec, freq_Hz) 
 98          return output 
 99      elif dtype == numpy.float64: 
100          output = numpy.zeros(shape=samples.shape, dtype=numpy.float64) 
101          qubfast.qub_filterdatad(samples.ctypes.data_as(c_double_p), output.ctypes.data_as(c_double_p), 
102                                  len(samples), sampling_sec, freq_Hz) 
103          return output 
104      else: 
105          input = numpy.array(samples, dtype=numpy.float32) 
106          output = numpy.zeros(shape=samples.shape, dtype=numpy.float32) 
107          qubfast.qub_filterdataf(input.ctypes.data_as(c_float_p), output.ctypes.data_as(c_float_p), 
108                                  len(samples), sampling_sec, freq_Hz) 
109          return output 
110   
111   
112  qubfast.DFII32_New.argtypes = () 
113  qubfast.DFII32_New.restype = c_void_p 
114  qubfast.DFII32_Free.argtypes = (c_void_p,) 
115  qubfast.DFII32_Free.restype = None 
116  qubfast.DFII32_Init.argtypes = (c_char_p, c_void_p, c_float) 
117  qubfast.DFII32_Init.restype = c_int 
118  qubfast.DFII32.argtypes = (c_float, c_void_p) 
119  qubfast.DFII32.restype = c_float 
120   
121 -class DFII32(object): 
122 -    def __init__(self): 
123          self.obj = qubfast.DFII32_New() 
124 -    def __del__(self): 
125          qubfast.DFII32_Free(self.obj) 
126 -    def init(self, filename, initial_val=0.0): 
127          return bool(qubfast.DFII32_Init(filename, self.obj, initial_val)) 
128 -    def next(self, x): 
129          return qubfast.DFII32(x, self.obj) 
130 -    def nextT(self, x): 
131          return qubfast.DFIIT32(x, self.obj) 
132   
133   
134  qubfast.DFII64_New.argtypes = () 
135  qubfast.DFII64_New.restype = c_void_p 
136  qubfast.DFII64_Free.argtypes = (c_void_p,) 
137  qubfast.DFII64_Free.restype = None 
138  qubfast.DFII64_Init.argtypes = (c_char_p, c_void_p, c_double) 
139  qubfast.DFII64_Init.restype = c_int 
140  qubfast.DFII64.argtypes = (c_double, c_void_p) 
141  qubfast.DFII64.restype = c_double 
142   
143 -class DFII64(object): 
144 -    def __init__(self): 
145          self.obj = qubfast.DFII64_New() 
146 -    def __del__(self): 
147          qubfast.DFII64_Free(self.obj) 
148 -    def init(self, filename, initial_val=0.0): 
149          return bool(qubfast.DFII64_Init(filename, self.obj, initial_val)) 
150 -    def next(self, x): 
151          return qubfast.DFII64(x, self.obj) 
152 -    def nextT(self, x): 
153          return qubfast.DFIIT64(x, self.obj) 
154