Changeset 867 for SHX/trunk


Ignore:
Timestamp:
11/22/12 11:26:49 (8 years ago)
Author:
marcus
Message:
  • more on filters
Location:
SHX/trunk/SeismicHandler
Files:
3 added
2 edited

Legend:

Unmodified
Added
Removed
  • SHX/trunk/SeismicHandler/commands/filter_.py

    r864 r867  
    66 
    77import numpy as np 
    8 import ctypes as C 
    98from SeismicHandler.basics.command import BaseCommand 
    109from SeismicHandler.modules.traces import traces_from_list, add_traces, \ 
     
    1211from SeismicHandler.basics.messages import log_message, get_runtime 
    1312from SeismicHandler.basics.error import ShxError 
    14 from SeismicHandler.lib.headers import COMPLEX, FFT_RATFCT, shc 
    15 from obspy.signal.util import nextpow2 
     13 
    1614 
    1715provides = {"filter": "_filter"} 
     
    115113                } 
    116114 
    117                 y = _filterFLF(f, tf.data, tf.stats.delta) 
    118115                tf.simulate( 
    119116                    paz_simulate=simfilter, 
     
    122119                    shsim=True 
    123120                ) 
    124                 tf.data = y 
    125121 
    126122                msg = "%s filtered by %s" % (tf.id, f.filename) 
     
    137133 
    138134 
    139 def _filterFLF(filter, data, sampling_rate, taper=1.0): 
    140     """ 
    141     Apply a Infinite Impulse response (IIR) filter. (FLF file) 
    142     """ 
    143     dt = 1.0 / float(sampling_rate) 
    144     serieslth = nextpow2(len(data)) 
    145     serieslth_2 = serieslth / 2 
    146     # add extra entry as realft starts with index 1! Numerical Recipies 
    147     y = np.zeros(serieslth + 1, dtype='float32') 
    148     y[1:len(data) + 1] = data # copy the data, data are casted to dtype of y 
    149     # perform FFT 
    150     shc.nr_realft(y, serieslth_2, 1) 
    151     # remove dummy, Y is a reference to y 
    152     Y = y[1:] 
    153     # set every second entry to it negative 
    154     # start with 3th entry, zeroth is offset, first is nyquist 
    155     # second is first real value, third is first imag value 
    156     Y[3::2] = -Y[3::2] 
    157     # create structure for transfer function 
    158     # poles and zeros cannot be casted to that structure 
    159     # they are not casted in C either 
    160     fct_ffv = FFT_RATFCT() 
    161     fct_ffv.norm = filter.normalization 
    162     fct_ffv.zero = (COMPLEX * 25) () 
    163     fct_ffv.pole = (COMPLEX * 25) () 
    164     fct_ffv.no_of_zeros = len(filter.input_coeff) 
    165     fct_ffv.no_of_poles = len(filter.output_coeff) 
    166     # set poles & zeros 
    167     for pos, coeff in enumerate(filter.input_coeff): 
    168         fct_ffv.zero[pos] = COMPLEX(coeff.real, coeff.imag) 
    169     for pos, coeff in enumerate(filter.output_coeff): 
    170         fct_ffv.pole[pos] = COMPLEX(coeff.real, coeff.imag) 
    171     # apply filter 
    172     shc.ff_frqmul(serieslth_2, C.cast(Y.ctypes.data, C.POINTER(COMPLEX)), 
    173                    2.0 * np.pi / (float(serieslth) * dt), 1, fct_ffv) 
    174     # taper the trace 
    175     if  taper >= 0.0 and taper < 1.0: 
    176         shc.ff_costaper(taper, serieslth_2, 
    177                          C.cast(Y.ctypes.data, C.POINTER(COMPLEX))) 
    178     # back transformation 
    179     Y[3::2] = -Y[3::2] 
    180     # now we need the dummy index again, take the y reference 
    181     shc.nr_realft(y, serieslth_2, -1); 
    182     # cut to the correct length 
    183     y = y[1:len(data)+1] 
    184     # save some memory, Y still has serieslth size 
    185     del Y 
    186     # normalize due to Numerical Recipies fft 
    187     y *= 2.0 / float(serieslth) 
    188     return y 
    189  
    190  
    191135if __name__ == "__main__": 
    192136    import doctest 
  • SHX/trunk/SeismicHandler/lib/headers.py

    r864 r867  
    363363 
    364364# filter operations 
    365 shc.ff_frqmul.argtypes = [ 
    366     C.c_long, C.POINTER(COMPLEX), C.c_float, C.c_int, C.POINTER(FFT_RATFCT) 
    367 ] 
    368 shc.ff_costaper.argtypes = [ 
    369     C.c_float, C.c_long, C.POINTER(COMPLEX) 
    370 ] 
    371 shc.nr_realft.argtypes = [ 
    372     np.ctypeslib.ndpointer(np.float32, flags="C_CONTIGUOUS"), C.c_long, 
    373     C.c_int 
    374 ] 
     365#shc.ff_frqmul.argtypes = [ 
     366#    C.c_long, C.POINTER(COMPLEX), C.c_float, C.c_int, C.POINTER(FFT_RATFCT) 
     367#] 
     368#shc.ff_costaper.argtypes = [ 
     369#    C.c_float, C.c_long, C.POINTER(COMPLEX) 
     370#] 
     371#shc.nr_realft.argtypes = [ 
     372#    np.ctypeslib.ndpointer(np.float32, flags="C_CONTIGUOUS"), C.c_long, 
     373#    C.c_int 
     374#] 
Note: See TracChangeset for help on using the changeset viewer.