Changeset 864 for SHX/trunk


Ignore:
Timestamp:
11/15/12 17:50:53 (8 years ago)
Author:
marcus
Message:
  • updating filter command, not yet finished
Location:
SHX/trunk/SeismicHandler
Files:
4 edited

Legend:

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

    r841 r864  
    66 
    77import numpy as np 
     8import ctypes as C 
    89from SeismicHandler.basics.command import BaseCommand 
    910from SeismicHandler.modules.traces import traces_from_list, add_traces, \ 
     
    1112from SeismicHandler.basics.messages import log_message, get_runtime 
    1213from SeismicHandler.basics.error import ShxError 
    13  
     14from SeismicHandler.lib.headers import COMPLEX, FFT_RATFCT, shc 
     15from obspy.signal.util import nextpow2 
    1416 
    1517provides = {"filter": "_filter"} 
     
    113115                } 
    114116 
     117                y = _filterFLF(f, tf.data, tf.stats.delta) 
    115118                tf.simulate( 
    116119                    paz_simulate=simfilter, 
     
    119122                    shsim=True 
    120123                ) 
     124                tf.data = y 
    121125 
    122126                msg = "%s filtered by %s" % (tf.id, f.filename) 
     
    133137 
    134138 
     139def _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 
    135191if __name__ == "__main__": 
    136192    import doctest 
  • SHX/trunk/SeismicHandler/lib/headers.py

    r789 r864  
    163163    ] 
    164164 
     165 
     166class COMPLEX(C.Structure): 
     167    """ 
     168    Complex number structure. 
     169 
     170    typedef struct sh_complex { 
     171        REAL    re;    /* real part */ 
     172        REAL    im;    /* imaginary part */ 
     173    } COMPLEX; 
     174    """ 
     175 
     176    _fields_ = [ 
     177        ("re", C.c_float), 
     178        ("im", C.c_float), 
     179    ] 
     180 
     181 
     182class FFT_RATFCT(C.Structure): 
     183    """ 
     184    Mapper for transfer functions (Seismic Handler, ffusrdef.h) 
     185 
     186    #define FFC_MAXDEGREE 25 
     187    typedef struct { 
     188        REAL     norm;                 /* normalisation */ 
     189        int      no_of_zeroes;         /* number of zeroes */ 
     190        COMPLEX  zero[FFC_MAXDEGREE];  /* zeroes of transfer function */ 
     191        int      no_of_poles;          /* number of poles */ 
     192        COMPLEX  pole[FFC_MAXDEGREE];  /* poles of transfer function */ 
     193    } FFT_RATFCT; 
     194    """ 
     195 
     196    _fields_ = [ 
     197        ("norm", C.c_float), 
     198        ("no_of_zeros", C.c_int), 
     199        ("zero", COMPLEX * 25), 
     200        ("no_of_poles", C.c_int), 
     201        ("pole", COMPLEX * 25), 
     202    ] 
     203 
     204 
    165205# prototypes 
    166206 
     
    321361    np.ctypeslib.ndpointer(np.float32, flags="C_CONTIGUOUS"), 
    322362] 
     363 
     364# filter operations 
     365shc.ff_frqmul.argtypes = [ 
     366    C.c_long, C.POINTER(COMPLEX), C.c_float, C.c_int, C.POINTER(FFT_RATFCT) 
     367] 
     368shc.ff_costaper.argtypes = [ 
     369    C.c_float, C.c_long, C.POINTER(COMPLEX) 
     370] 
     371shc.nr_realft.argtypes = [ 
     372    np.ctypeslib.ndpointer(np.float32, flags="C_CONTIGUOUS"), C.c_long, 
     373    C.c_int 
     374] 
  • SHX/trunk/SeismicHandler/modules/filters.py

    r798 r864  
    1 #! /usr/bin/env python 
    21# -*- coding: utf-8 -*- 
    32 
Note: See TracChangeset for help on using the changeset viewer.