Python scipy.signal 模块,iirfilter() 实例源码

我们从Python开源项目中,提取了以下22个代码示例,用于说明如何使用scipy.signal.iirfilter()

项目:NeoAnalysis    作者:neoanalysis    | 项目源码 | 文件源码
def __band_filter(self,ite_data,fs,order,lowcut,highcut,zerophase,btype,ftype,rps=None):
        fe = fs/2.0
        low = lowcut/fe
        high = highcut/fe
        if low<0:
            low=0
        if high>1:
            high=1
        if ftype == "cheby1":
            rp = rps
            z,p,k = signal.iirfilter(order,[low,high],btype=btype,ftype=ftype,output="zpk",rp=rp)
        elif ftype == "cheby2":
            rs = rps
            z,p,k = signal.iirfilter(order,[low,high],btype=btype,ftype=ftype,output="zpk",rs=rs)
        elif ftype == "ellip":
            rp = rps[0]
            rs = rps[1]
            z,p,k = signal.iirfilter(order,[low,high],btype=btype,ftype=ftype,output="zpk",rp=rp,rs=rs)
        else:
            z,p,k = signal.iirfilter(order,[low,high],btype=btype,ftype=ftype,output="zpk")
        sos = signal.zpk2sos(z,p,k)
        ite_data = signal.sosfilt(sos,ite_data)
        if zerophase:
            ite_data = signal.sosfilt(sos,ite_data[::-1])[::-1]
        return ite_data
项目:NeoAnalysis    作者:neoanalysis    | 项目源码 | 文件源码
def __highpass_filter(self,ite_data,fs,order,lowcut,zerophase,btype,ftype,rps=None):
        fe = fs/2.0
        low = lowcut/fe
        if low<0:
            low=0
        if ftype == "cheby1":
            rp = rps
            z,p,k = signal.iirfilter(order,low,btype=btype,ftype=ftype,output="zpk",rp=rp)
        elif ftype == "cheby2":
            rs = rps
            z,p,k = signal.iirfilter(order,low,btype=btype,ftype=ftype,output="zpk",rs=rs)
        elif ftype == "ellip":
            rp = rps[0]
            rs = rps[1]
            z,p,k = signal.iirfilter(order,low,btype=btype,ftype=ftype,output="zpk",rp=rp,rs=rs)
        else:
            z,p,k = signal.iirfilter(order,low,btype=btype,ftype=ftype,output="zpk")
        sos = signal.zpk2sos(z,p,k)
        ite_data = signal.sosfilt(sos,ite_data)
        if zerophase:
            ite_data = signal.sosfilt(sos,ite_data[::-1])[::-1]
        return ite_data
项目:NeoAnalysis    作者:neoanalysis    | 项目源码 | 文件源码
def __lowpass_filter(self,ite_data,fs,order,highcut,zerophase,btype,ftype,rps=None):
        fe = fs/2.0
        high = highcut/fe
        if high>1:
            high=1
        if ftype == "cheby1":
            rp = rps
            z,p,k = signal.iirfilter(order,high,btype=btype,ftype=ftype,output="zpk",rp=rp)
        elif ftype == "cheby2":
            rs = rps
            z,p,k = signal.iirfilter(order,high,btype=btype,ftype=ftype,output="zpk",rs=rs)
        elif ftype == "ellip":
            rp = rps[0]
            rs = rps[1]
            z,p,k = signal.iirfilter(order,high,btype=btype,ftype=ftype,output="zpk",rp=rp,rs=rs)
        else:
            z,p,k = signal.iirfilter(order,high,btype=btype,ftype=ftype,output="zpk")
        sos = signal.zpk2sos(z,p,k)
        ite_data = signal.sosfilt(sos,ite_data)
        if zerophase:
            ite_data = signal.sosfilt(sos,ite_data[::-1])[::-1]
        return ite_data
项目:NeoAnalysis    作者:neoanalysis    | 项目源码 | 文件源码
def __lowpass(self,ite_data,highcut,fs,order,ftype,zerophase,**args):
        fe = fs/2.0
        high = highcut/fe
        if high>1:
            high=1

        if ftype == "cheby1":
            rp = args["rp"]
            z,p,k = signal.iirfilter(order,high,btype="lowpass",ftype=ftype,output="zpk",rp=rp)
        elif ftype == "cheby2":
            rs = args["rs"]
            z,p,k = signal.iirfilter(order,high,btype="lowpass",ftype=ftype,output="zpk",rs=rs)
        elif ftype == "ellip":
            rp = args["rp"]
            rs = args["rs"]
            z,p,k = signal.iirfilter(order,high,btype="lowpass",ftype=ftype,output="zpk",rp=rp,rs=rs)
        else:
            z,p,k = signal.iirfilter(order,high,btype="lowpass",ftype=ftype,output="zpk")
        sos = signal.zpk2sos(z,p,k)
        ite_data = signal.sosfilt(sos,ite_data)
        if zerophase:
            ite_data = signal.sosfilt(sos,ite_data[::-1])[::-1]
        return ite_data
项目:NeoAnalysis    作者:neoanalysis    | 项目源码 | 文件源码
def band_pass(ite_data,freqmin,freqmax,samp_freq,corners=32,zerophase=True):
    fe = samp_freq/2.0
    low = freqmin/fe
    high = freqmax/fe
    # raise error for illegal input
    if high - 1.0 > -1e-6:
        msg = ("Selected high corner frequency ({}) of bandpass is at or above Nyquist ({}). Applying a high-pass instead.").format(freqmax, fe)
        return False
    if low >1 :
        msg = "selected low corner requency is above Nyquist"
        return False

    z,p,k = signal.iirfilter(corners,[low,high],btype='band',ftype='butter',output='zpk')
    sos = zpk2sos(z,p,k)
    ite_data = sosfilt(sos,ite_data)
    if zerophase:
        ite_data = sosfilt(sos,ite_data[::-1])[::-1]
    return ite_data
项目:NeoAnalysis    作者:neoanalysis    | 项目源码 | 文件源码
def band_stop(ite_data,freqmin,freqmax,samp_freq,corners=4,zerophase=True):
    fe = samp_freq/2.0
    low = freqmin/fe
    high = freqmax/fe
    # raise error for illegal input
    if high - 1.0 > -1e-6:
        msg = ("Selected high corner frequency ({}) of bandpass is at or above Nyquist ({}). Applying a high-pass instead.").format(freqmax, fe)
        return False
    if low >1 :
        msg = "selected low corner requency is above Nyquist"
        return False

    z,p,k = signal.iirfilter(corners,[low,high],btype='bandstop',ftype='butter',output='zpk')
    sos = zpk2sos(z,p,k)
    ite_data = sosfilt(sos,ite_data)
    if zerophase:
        ite_data = sosfilt(sos,ite_data[::-1])[::-1]
    return ite_data
项目:NeoAnalysis    作者:neoanalysis    | 项目源码 | 文件源码
def __band_filter(self,ite_data,fs,order,lowcut,highcut,zerophase,btype,ftype,rps=None):
        fe = fs/2.0
        low = lowcut/fe
        high = highcut/fe
        if low<0:
            low=0
        if high>1:
            high=1
        if ftype == "cheby1":
            rp = rps
            z,p,k = signal.iirfilter(order,[low,high],btype=btype,ftype=ftype,output="zpk",rp=rp)
        elif ftype == "cheby2":
            rs = rps
            z,p,k = signal.iirfilter(order,[low,high],btype=btype,ftype=ftype,output="zpk",rs=rs)
        elif ftype == "ellip":
            rp = rps[0]
            rs = rps[1]
            z,p,k = signal.iirfilter(order,[low,high],btype=btype,ftype=ftype,output="zpk",rp=rp,rs=rs)
        else:
            z,p,k = signal.iirfilter(order,[low,high],btype=btype,ftype=ftype,output="zpk")
        sos = signal.zpk2sos(z,p,k)
        ite_data = signal.sosfilt(sos,ite_data)
        if zerophase:
            ite_data = signal.sosfilt(sos,ite_data[::-1])[::-1]
        return ite_data
项目:NeoAnalysis    作者:neoanalysis    | 项目源码 | 文件源码
def __highpass_filter(self,ite_data,fs,order,lowcut,zerophase,btype,ftype,rps=None):
        fe = fs/2.0
        low = lowcut/fe
        if low<0:
            low=0
        if ftype == "cheby1":
            rp = rps
            z,p,k = signal.iirfilter(order,low,btype=btype,ftype=ftype,output="zpk",rp=rp)
        elif ftype == "cheby2":
            rs = rps
            z,p,k = signal.iirfilter(order,low,btype=btype,ftype=ftype,output="zpk",rs=rs)
        elif ftype == "ellip":
            rp = rps[0]
            rs = rps[1]
            z,p,k = signal.iirfilter(order,low,btype=btype,ftype=ftype,output="zpk",rp=rp,rs=rs)
        else:
            z,p,k = signal.iirfilter(order,low,btype=btype,ftype=ftype,output="zpk")
        sos = signal.zpk2sos(z,p,k)
        ite_data = signal.sosfilt(sos,ite_data)
        if zerophase:
            ite_data = signal.sosfilt(sos,ite_data[::-1])[::-1]
        return ite_data
项目:NeoAnalysis    作者:neoanalysis    | 项目源码 | 文件源码
def __highpass(self,ite_data,lowcut,fs,order,ftype,zerophase,**args):
        fe = fs/2.0
        low = lowcut/fe
        if low<0:
            low=0

        if ftype == "cheby1":
            rp = args["rp"]
            z,p,k = signal.iirfilter(order,low,btype="highpass",ftype=ftype,output="zpk",rp=rp)
        elif ftype == "cheby2":
            rs = args["rs"]
            z,p,k = signal.iirfilter(order,low,btype="highpass",ftype=ftype,output="zpk",rs=rs)
        elif ftype == "ellip":
            rp = args["rp"]
            rs = args["rs"]
            z,p,k = signal.iirfilter(order,low,btype="highpass",ftype=ftype,output="zpk",rp=rp,rs=rs)
        else:
            z,p,k = signal.iirfilter(order,low,btype="highpass",ftype=ftype,output="zpk")
        sos = signal.zpk2sos(z,p,k)
        ite_data = signal.sosfilt(sos,ite_data)
        if zerophase:
            ite_data = signal.sosfilt(sos,ite_data[::-1])[::-1]
        return ite_data
项目:NeoAnalysis    作者:neoanalysis    | 项目源码 | 文件源码
def __lowpass(self,ite_data,highcut,fs,order,ftype,zerophase,**args):
        fe = fs/2.0
        high = highcut/fe
        if high>1:
            high=1

        if ftype == "cheby1":
            rp = args["rp"]
            z,p,k = signal.iirfilter(order,high,btype="lowpass",ftype=ftype,output="zpk",rp=rp)
        elif ftype == "cheby2":
            rs = args["rs"]
            z,p,k = signal.iirfilter(order,high,btype="lowpass",ftype=ftype,output="zpk",rs=rs)
        elif ftype == "ellip":
            rp = args["rp"]
            rs = args["rs"]
            z,p,k = signal.iirfilter(order,high,btype="lowpass",ftype=ftype,output="zpk",rp=rp,rs=rs)
        else:
            z,p,k = signal.iirfilter(order,high,btype="lowpass",ftype=ftype,output="zpk")
        sos = signal.zpk2sos(z,p,k)
        ite_data = signal.sosfilt(sos,ite_data)
        if zerophase:
            ite_data = signal.sosfilt(sos,ite_data[::-1])[::-1]
        return ite_data
项目:NeoAnalysis    作者:neoanalysis    | 项目源码 | 文件源码
def band_pass(ite_data,freqmin,freqmax,samp_freq,corners=32):
    fe = samp_freq/2.0
    low = freqmin/fe
    high = freqmax/fe
    # raise error for illegal input
    if high - 1.0 > -1e-6:
        msg = ("Selected high corner frequency ({}) of bandpass is at or above Nyquist ({}). Applying a high-pass instead.").format(freqmax, fe)
        return False
    if low >1 :
        msg = "selected low corner requency is above Nyquist"
        return False

    z,p,k = signal.iirfilter(corners,[low,high],btype='band',ftype='butter',output='zpk')
    sos = zpk2sos(z,p,k)
    ite_data = sosfilt(sos,ite_data)
    ite_data = sosfilt(sos,ite_data[::-1])[::-1]
    return ite_data
项目:NeoAnalysis    作者:neoanalysis    | 项目源码 | 文件源码
def __bandpass(self,ite_data,lowcut,highcut,fs,order,ftype,zerophase,**args):
        fe = fs/2.0
        low = lowcut/fe
        high = highcut/fe
        if low<0:
            low=0
        if high>1:
            high=1

        if ftype == "cheby1":
            rp = args["rp"]
            z,p,k = signal.iirfilter(order,[low,high],btype="band",ftype=ftype,output="zpk",rp=rp)
        elif ftype == "cheby2":
            rs = args["rs"]
            z,p,k = signal.iirfilter(order,[low,high],btype="band",ftype=ftype,output="zpk",rs=rs)
        elif ftype == "ellip":
            rp = args["rp"]
            rs = args["rs"]
            z,p,k = signal.iirfilter(order,[low,high],btype="band",ftype=ftype,output="zpk",rp=rp,rs=rs)
        else:
            z,p,k = signal.iirfilter(order,[low,high],btype="band",ftype=ftype,output="zpk")
        sos = signal.zpk2sos(z,p,k)
        ite_data = signal.sosfilt(sos,ite_data)
        if zerophase:
            ite_data = signal.sosfilt(sos,ite_data[::-1])[::-1]
        return ite_data
项目:NeoAnalysis    作者:neoanalysis    | 项目源码 | 文件源码
def __bandstop(self,ite_data,lowcut,highcut,fs,order,ftype,zerophase,**args):
        fe = fs/2.0
        low = lowcut/fe
        high = highcut/fe
        if low<0:
            low=0
        if high>1:
            high=1

        if ftype == "cheby1":
            rp = args["rp"]
            z,p,k = signal.iirfilter(order,[low,high],btype="bandstop",ftype=ftype,output="zpk",rp=rp)
        elif ftype == "cheby2":
            rs = args["rs"]
            z,p,k = signal.iirfilter(order,[low,high],btype="bandstop",ftype=ftype,output="zpk",rs=rs)
        elif ftype == "ellip":
            rp = args["rp"]
            rs = args["rs"]
            z,p,k = signal.iirfilter(order,[low,high],btype="bandstop",ftype=ftype,output="zpk",rp=rp,rs=rs)
        else:
            z,p,k = signal.iirfilter(order,[low,high],btype="bandstop",ftype=ftype,output="zpk")

        sos = signal.zpk2sos(z,p,k)
        ite_data = signal.sosfilt(sos,ite_data)
        if zerophase:
            ite_data = signal.sosfilt(sos,ite_data[::-1])[::-1]
        return ite_data
项目:NeoAnalysis    作者:neoanalysis    | 项目源码 | 文件源码
def __bandpass(self,ite_data,lowcut,highcut,fs,order,ftype,zerophase,**args):
        fe = fs/2.0
        low = lowcut/fe
        high = highcut/fe
        if low<0:
            low=0
        if high>1:
            high=1

        if ftype == "cheby1":
            rp = args["rp"]
            z,p,k = signal.iirfilter(order,[low,high],btype="band",ftype=ftype,output="zpk",rp=rp)
        elif ftype == "cheby2":
            rs = args["rs"]
            z,p,k = signal.iirfilter(order,[low,high],btype="band",ftype=ftype,output="zpk",rs=rs)
        elif ftype == "ellip":
            rp = args["rp"]
            rs = args["rs"]
            z,p,k = signal.iirfilter(order,[low,high],btype="band",ftype=ftype,output="zpk",rp=rp,rs=rs)
        else:
            z,p,k = signal.iirfilter(order,[low,high],btype="band",ftype=ftype,output="zpk")
        sos = signal.zpk2sos(z,p,k)
        ite_data = signal.sosfilt(sos,ite_data)
        if zerophase:
            ite_data = signal.sosfilt(sos,ite_data[::-1])[::-1]
        return ite_data
项目:NeoAnalysis    作者:neoanalysis    | 项目源码 | 文件源码
def __bandstop(self,ite_data,lowcut,highcut,fs,order,ftype,zerophase,**args):
        fe = fs/2.0
        low = lowcut/fe
        high = highcut/fe
        if low<0:
            low=0
        if high>1:
            high=1

        if ftype == "cheby1":
            rp = args["rp"]
            z,p,k = signal.iirfilter(order,[low,high],btype="bandstop",ftype=ftype,output="zpk",rp=rp)
        elif ftype == "cheby2":
            rs = args["rs"]
            z,p,k = signal.iirfilter(order,[low,high],btype="bandstop",ftype=ftype,output="zpk",rs=rs)
        elif ftype == "ellip":
            rp = args["rp"]
            rs = args["rs"]
            z,p,k = signal.iirfilter(order,[low,high],btype="bandstop",ftype=ftype,output="zpk",rp=rp,rs=rs)
        else:
            z,p,k = signal.iirfilter(order,[low,high],btype="bandstop",ftype=ftype,output="zpk")

        sos = signal.zpk2sos(z,p,k)
        ite_data = signal.sosfilt(sos,ite_data)
        if zerophase:
            ite_data = signal.sosfilt(sos,ite_data[::-1])[::-1]
        return ite_data
项目:kaggle-seizure-prediction    作者:sics-lm    | 项目源码 | 文件源码
def apply_ellip_filter(self, data, order, cutoff, btype):
        ripple = 3
        attenuation = 50
        cutoff /= self.nyq
        b, a = iirfilter(N=order, Wn=cutoff, rp=ripple, rs=attenuation,
                         btype=btype, ftype='ellip')
        return lfilter(b, a, data, axis=0)
项目:arlpy    作者:org-arl    | 项目源码 | 文件源码
def test_lfilter_gen(self):
        x = np.random.normal(0, 1, 1000)
        hb = np.array([0, 0, 1, 0], dtype=np.float)
        f = signal.lfilter_gen(hb, 1)
        y = [f.send(v) for v in x]
        self.assertArrayEqual(np.append([0, 0], x[:-2]), y)
        hb, ha = sp.iirfilter(4, 0.01, btype='lowpass')
        y1 = sp.lfilter(hb, ha, x)
        f = signal.lfilter_gen(hb, ha)
        y2 = [f.send(v) for v in x]
        self.assertArrayEqual(y1, y2, precision=6)
项目:cebl    作者:idfah    | 项目源码 | 文件源码
def __init__(self, lowFreq, highFreq, sampRate=1.0,
                 order=3, filtType='butter', zeroPhase=True,
                 dtype=None, **kwargs):
        """Construct a new IIR bandpass filter.
        """
        BandpassFilterBase.__init__(self, lowFreq, highFreq, sampRate, dtype=dtype)

        self.order     = order
        self.filtType  = filtType.lower()
        self.zeroPhase = zeroPhase

        if self.bandType not in ('allpass', 'allstop'):
            if self.bandType == 'lowpass':
                self.Wn = self.high
            elif self.bandType == 'highpass':
                self.Wn = self.low
            elif self.bandType == 'bandpass':
                self.Wn = (self.low, self.high)
            elif self.bandType == 'bandstop':
                self.Wn = (self.high, self.low)
            else:
                raise Exception('Invalid bandType: ' + str(self.bandType))

            self.numCoef, self.denomCoef = spsig.iirfilter(order, self.Wn,
                ftype=filtType, btype=self.bandType, **kwargs)

            self.numCoef = self.numCoef.astype(self.dtype, copy=False)
            self.denomCoef = self.denomCoef.astype(self.dtype, copy=False)

            self.initZi()
项目:yam    作者:trichter    | 项目源码 | 文件源码
def _filter_resp(freqmin, freqmax, corners=2, zerophase=False, sr=None,
                 N=None, whole=False):
    """
    Complex frequency response of Butterworth-Bandpass Filter.

    :param freqmin: Pass band low corner frequency.
    :param freqmax: Pass band high corner frequency.
    :param corners: Filter corners
    :param zerophase: If True, apply filter once forwards and once backwards.
        This results in twice the number of corners but zero phase shift in
        the resulting filtered trace.
    :param sr: Sampling rate in Hz.
    :param N,whole: passed to scipy.signal.freqz

    :return: frequencies and complex response
    """
    df = sr
    fe = 0.5 * df
    low = freqmin / fe
    high = freqmax / fe
    # raise for some bad scenarios
    if high > 1:
        high = 1.0
        msg = "Selected high corner frequency is above Nyquist. " + \
              "Setting Nyquist as high corner."
        log.warning(msg)
    if low > 1:
        msg = "Selected low corner frequency is above Nyquist."
        raise ValueError(msg)
    [b, a] = iirfilter(corners, [low, high], btype='band',
                       ftype='butter', output='ba')
    freqs, values = freqz(b, a, N, whole=whole)
    if zerophase:
        values *= np.conjugate(values)
    return freqs, values
项目:seis_tools    作者:romaguir    | 项目源码 | 文件源码
def filter_freqs(lowcut, highcut, fs, plot=False, corners=4):
   '''
   The frequency band information should technically be the amplitude spectrum
   of the filtered seismograms, but the amplitude spectrum of the bandpass
   filter is sufficient. Here we use a bandpass butterworth filter.

   The function freq band takes 3 arguments:
           lowcut = low end of bandpass (in Hz)
           highcut = high end of bandpass (in Hz)
           fs      = sample rate (1.0/dt)

   It returns two vectors:
           omega = frequency axis (in rad/s)
           amp   = frequency response of the filter
   '''
   #Define bandpass parameters
   nyquist = 0.5 * fs
   fmin    = lowcut/nyquist
   fmax    = highcut/nyquist

   #Make filter shape
   b, a = iirfilter(corners, [fmin,fmax], btype='band', ftype='butter')

   #Determine freqency response
   freq_range = np.linspace(0,0.15,200)
   w, h = freqz(b,a,worN=freq_range)

   omega    = fs * w                # in rad/s
   omega_hz = (fs * w) / (2*np.pi)  # in hz
   amp   = abs(h)

   if(plot == True):
      #Checks----------------
      #plt.semilogx(omega,amp)
      #plt.axvline(1/10.0)
      #plt.axvline(1/25.0)
      #plt.axhline(np.sqrt(0.5))

      plt.plot(omega,amp)
      plt.xlabel('frequency (rad/s)')
      plt.ylabel('amplitude')
      plt.show()

   return  omega, amp
项目:seis_tools    作者:romaguir    | 项目源码 | 文件源码
def get_filter_freqs(filter_type,freqmin,freqmax,sampling_rate,**kwargs):
   '''
   Returns frequency band information of the filter used in cross
   correlation delay time measurements.

   args-------------------------------------------------------------------------- 
   filter_type: type of filter used (e.g., bandpass, gaussian etc...)
   freqmin: minimum frequency of filter (in Hz)
   freqmax: maximum frequency of filter (in Hz)
   sampling_rate: sampling rate of seismograms (in Hz)
   kwargs-------------------------------------------------------------------------
   corners: number of corners used in filter (if bandpass). default = 2
   std_dev: standard deviation of gaussian filter (in Hz)
   plot: True or False

   returns-----------------------------------------------------------------------
   omega: frequency axis (rad/s)
   amp: frequency response of filter
   '''
   plot = kwargs.get('plot',False)
   corners = kwargs.get('corners',4)
   std_dev = kwargs.get('std_dev',0.1)
   mid_freq = kwargs.get('mid_freq',1/10.0)

   if filter_type == 'bandpass':
      nyquist = 0.5 * sampling_rate
      fmin = freqmin/nyquist
      fmax = freqmax/nyquist
      b, a = iirfilter(corners, [fmin,fmax], btype='band', ftype='butter')
      freq_range = np.linspace(0,0.15,200)
      w, h = freqz(b,a,worN=freq_range)
      omega = sampling_rate * w
      omega_hz = (sampling_rate * w) / (2*np.pi)
      amp = abs(h)

   elif filter_type == 'gaussian':
      fmin=freqmin
      fmax=freqmax
      omega_hz = np.linspace(0,0.5,200)
      omega = omega_hz*(2*np.pi)
      f_middle_hz = (fmin+fmax)/2.0
      f_middle = f_middle_hz*(2*np.pi)
      #f_middle_hz = mid_freq  #f_middle_hz = 10.0 was used in JGR manuscript
      #f_middle = f_middle_hz*(2*np.pi)
      print f_middle_hz,f_middle
      amp = np.exp(-1.0*((omega-f_middle)**2)/(2*(std_dev**2)))
      amp = amp/np.max(amp)

   if plot:
      fig,axes = plt.subplots(2)
      axes[0].plot(omega,amp)
      axes[0].set_xlabel('frequency (rad/s)')
      axes[0].set_ylabel('amplitude')
      axes[1].plot(omega_hz,amp)
      axes[1].set_xlabel('frequency (Hz)')
      axes[1].set_ylabel('amplitude')
      axes[1].axvline(fmin)
      axes[1].axvline(fmax)
      plt.show()

   return omega, amp
项目:seis_tools    作者:romaguir    | 项目源码 | 文件源码
def get_filter_freqs(filter_type,freqmin,freqmax,sampling_rate,**kwargs):
   '''
   Returns frequency band information of the filter used in cross
   correlation delay time measurements.

   args-------------------------------------------------------------------------- 
   filter_type: type of filter used (e.g., bandpass, gaussian etc...)
   freqmin: minimum frequency of filter (in Hz)
   freqmax: maximum frequency of filter (in Hz)
   sampling_rate: sampling rate of seismograms (in Hz)
   kwargs-------------------------------------------------------------------------
   corners: number of corners used in filter (if bandpass). default = 2
   std_dev: standard deviation of gaussian filter (in Hz)
   plot: True or False

   returns-----------------------------------------------------------------------
   omega: frequency axis (rad/s)
   amp: frequency response of filter
   '''
   plot = kwargs.get('plot',False)
   corners = kwargs.get('corners',4)
   std_dev = kwargs.get('std_dev',0.1)
   mid_freq = kwargs.get('mid_freq',1/10.0)

   if filter_type == 'bandpass':
      nyquist = 0.5 * sampling_rate
      fmin = freqmin/nyquist
      fmax = freqmax/nyquist
      b, a = iirfilter(corners, [fmin,fmax], btype='band', ftype='butter')
      freq_range = np.linspace(0,0.15,200)
      w, h = freqz(b,a,worN=freq_range)
      omega = sampling_rate * w
      omega_hz = (sampling_rate * w) / (2*np.pi)
      amp = abs(h)

   elif filter_type == 'gaussian':
      fmin=freqmin
      fmax=freqmax
      omega_hz = np.linspace(0,0.5,200)
      omega = omega_hz*(2*np.pi)
      f_middle_hz = (fmin+fmax)/2.0
      f_middle = f_middle_hz*(2*np.pi)
      #f_middle_hz = mid_freq  #f_middle_hz = 10.0 was used in JGR manuscript
      #f_middle = f_middle_hz*(2*np.pi)
      print f_middle_hz,f_middle
      amp = np.exp(-1.0*((omega-f_middle)**2)/(2*(std_dev**2)))
      amp = amp/np.max(amp)

   if plot:
      fig,axes = plt.subplots(2)
      axes[0].plot(omega,amp)
      axes[0].set_xlabel('frequency (rad/s)')
      axes[0].set_ylabel('amplitude')
      axes[1].plot(omega_hz,amp)
      axes[1].set_xlabel('frequency (Hz)')
      axes[1].set_ylabel('amplitude')
      axes[1].axvline(fmin)
      axes[1].axvline(fmax)
      plt.show()

   return omega, amp