Python signal.iirfilter函数代码示例

def preprocess_morlet(data, metadata, notch_filter=True, downsampling_factor=32, sfreq=250, freqs=[10], n_cycles=50):

    b1, a1 = signal.iirfilter(1, [59.0 / 125.0, 61.0 / 125.0], btype="bandstop")
    b2, a2 = signal.iirfilter(1, 3.0 / 125.0, btype="highpass")

    def do_morlet(arr):

        if notch_filter:
            arr = signal.lfilter(b1, a1, arr)
            arr = signal.lfilter(b2, a2, arr)

        cwt = wavelet_transform(
            arr[:, np.newaxis], sfreq=sfreq, freqs=freqs, n_cycles=n_cycles, include_phase=False, log_mag=True
        )

        arr = cwt.mean(axis=1)

        if downsampling_factor > 1:
            decimated = signal.decimate(arr, downsampling_factor)
        else:
            decimated = arr

        return decimated

    return preprocess_general(process=do_morlet, data=data, metadata=metadata)

def preprocess_stft(data, metadata, notch_filter=True,
                    box_width=128, downsampling_factor=32,
                    sfreq=250, low_f=8, high_f=12, kaiser_beta=14):

    b1, a1 = signal.iirfilter(1, [59.0/125.0, 61.0/125.0], btype='bandstop')
    b2, a2 = signal.iirfilter(1, 3.0/125.0, btype='highpass')

    def do_stft(arr):

        if notch_filter:
            arr = signal.lfilter(b1, a1, arr)
            arr = signal.lfilter(b2, a2, arr)
            if len(arr) < box_width:
                raise ValueError("The buffer_size used by the connector should "
                                 "be higher than box_width. buffer_size = "
                                 "%s | box_width = %s" % (len(arr), box_width))

        out = stft(arr[:, np.newaxis],
                   box_width=box_width, step=downsampling_factor, pad_width=0,
                   kaiser_beta=kaiser_beta, include_phase=False, log_mag=True)

        fftfreq = np.fft.rfftfreq(box_width, d=1/float(sfreq))
        good = np.logical_and(fftfreq >= low_f, fftfreq <= high_f)
        out = out[:, good].mean(axis=1)

        return out



    return preprocess_general(process=do_stft,
                              data=data, metadata=metadata)

def compute_theory(stala, stalo, evla, evlo, gamma_r, gamma_s, gamma_d, iiM, tt, t1, t2, channel):
    G = 6.674e-11
    R = 6370000.

    # station coordinates
    stala *= np.pi / 180
    stalo *= np.pi / 180

    # event coordinates
    evla *= np.pi / 180
    evlo *= np.pi /180

    rake   = np.array([[np.cos(gamma_r), -np.sin(gamma_r), 0], [np.sin(gamma_r), np.cos(gamma_r), 0], [0, 0, 1]])
    strike = np.array([[np.cos(gamma_s), np.sin(gamma_s), 0], [-np.sin(gamma_s), np.cos(gamma_s), 0], [0, 0, 1]])
    dip    = np.array([[np.cos(gamma_d), 0, -np.sin(gamma_d)], [0, 1, 0], [np.sin(gamma_d), 0, np.cos(gamma_d)]])

    ex = strike.dot(dip).dot(np.array([0, 0, 1]))
    ez = strike.dot(dip).dot(rake).dot(np.array([0, 1, 0]))

    rst = R * np.array([np.cos(stala)*np.sin(stalo), np.cos(stala)*np.cos(stalo), np.sin(stala)])
    rev = R * np.array([np.cos(evla)*np.sin(evlo), np.cos(evla)*np.cos(evlo), np.sin(evla)])
    dr  = rst - rev

    er  = rst / np.linalg.norm(rst)
    er0 =  dr / np.linalg.norm(dr)

    eNS = np.array([np.sin(stala)*np.sin(stalo), np.sin(stala)*np.cos(stalo), -np.cos(stala)])
    eNS = eNS / np.linalg.norm(eNS)

    eWE = np.array([np.cos(stalo), -np.sin(stalo), 0.])
    eWE = eWE / np.linalg.norm(eWE)

    a = ex.T.dot(er0)*ez
    b = ez.T.dot(er0)*ex
    c = ex.T.dot(er0)*ez.T.dot(er0)*er0

    vec = - 6*G / np.linalg.norm(dr)**4 * ( a + b - 5*c )
    th_UD = er.T.dot(vec) * iiM
    th_WE = eWE.T.dot(vec) * iiM
    th_NS = eNS.T.dot(vec) * iiM

    if channel == 'LHZ':
        th = th_UD
    elif channel == 'LHN':
        th = th_NS
    else:
        th = th_WE

    # filter band applied to the seismic data during preprocessing
    fcut1 = 0.001 / (0.5 * 1.0)
    [b, a] = signal.iirfilter(4, fcut1, btype='high', ftype='bessel', output='ba')
    th = signal.lfilter(b, a, th)

    fcut2 = 2*0.059904 / (0.5 * 1.0)
    [b, a] = signal.iirfilter(8, fcut2, btype='low', ftype='bessel', output='ba')
    th = signal.lfilter(b, a, th)

    mask = np.nonzero( (tt >= t1*np.ones(len(tt))) * (tt <= t2*np.ones(len(tt))) )
    return np.mean(th[mask]*1.0E8)

    def set_sampling_ferquency(self, fs, channels, bpfilter, notchfilter):
        """ Set the sampling frequency and filters for individual channels.

        Parameters:
        fs -- sampling frequency
        channels -- list of booleans: channels[0] == True: enable filter for channel 0
        bpfilter -- tuple: parameters for the band pass filter (hp, lp, fs, order) or None
        notchfilter -- tuple: parameters for the band stop filter (hp, lp, fs, order) or None

        """
        # we have: hp, lp, fs, order, typ
        # signal.iirfilter(order/2, [hp/(fs/2), lp/(fs/2)], ftype='butter', btype='band')
        # we get 18 coeffs and put them in as '<d' in the buffer
        # struct.pack('<'+'d'*18, *coeffs)

        # special filter: means no filter
        null_filter = "\x00\x00\x00\x00\x00\x00\xf0\x3f"+"\x00\x00\x00\x00\x00\x00\x00\x00"*17

        if bpfilter:
            bp_hp, bp_lp, bp_fs, bp_order = bpfilter
            bp_b, bp_a = iirfilter(bp_order/2, [bp_hp/(bp_fs/2), bp_lp/(bp_fs/2)], ftype='butter', btype='band')
            bp_filter = list(bp_b)
            bp_filter.extend(list(bp_a))
            bp_filter = struct.pack("<"+"d"*18, *bp_filter)
        else:
            bp_filter = null_filter

        if notchfilter:
            bs_hp, bs_lp, bs_fs, bs_order = notchfilter
            bs_b, bs_a = iirfilter(bs_order/2, [bs_hp/(bs_fs/2), bs_lp/(bs_fs/2)], ftype='butter', btype='bandstop')
            bs_filter = list(bs_b)
            # the notch filter has (always?) an order of 4 so fill the gaps with
            # zeros
            if len(bs_filter) < 9:
                diff = 9 - len(bs_filter)
                bs_filter.extend([0.0 for i in range(diff)])
            bs_filter.extend(list(bs_a))
            if len(bs_filter) < 18:
                diff = 18 - len(bs_filter)
                bs_filter.extend([0.0 for i in range(diff)])
            bs_filter = struct.pack("<"+"d"*18, *bs_filter)
        else:
            bs_filter = null_filter

        # set the filters for all channels
        if bpfilter == notchfilter == None:
            self.devh.controlMsg(CX_OUT, 0xc6, value=0x01, buffer=bp_filter)
            self.devh.controlMsg(CX_OUT, 0xc7, value=0x01, buffer=bs_filter)
        else:
            idx = 1
            for i in channels:
                if i:
                    self.devh.controlMsg(CX_OUT, 0xc6, value=idx, buffer=bp_filter)
                    self.devh.controlMsg(CX_OUT, 0xc7, value=idx, buffer=bs_filter)
                idx += 1

        # set the sampling frequency
        self.devh.controlMsg(CX_OUT, 0xb6, value=fs, buffer=0)

	def __init__(self, order=3, wn=0.01, coeff_scale=1, rp=1, rs=60, ftype='butter', dtype=np.float):
		if order%2 == 0:
			print 'No even order filters allowed.'
			return
		# Filter settings
		self.order = order
		self.wn = wn
		self.rp = rp
		self.rs = rs
		self.dtype = dtype
		self.ftype = ftype
		self.coeff_scale = coeff_scale
		gammas = []
		psi = np.tan(np.pi*self.wn/2.0)
		psi2 = psi*psi
		if self.ftype == 'butter':
			(z,p,k) = signal.iirfilter(self.order, psi, btype='lowpass', analog=1, ftype='butter', output='zpk')
			filter_roots = np.sort(p)
		elif self.ftype == 'bessel':
			print 'Please avoid using Bessel filters as they don\'t translate well to LWDFs.'
			(z,p,k) = signal.iirfilter(self.order, psi, btype='lowpass', analog=1, ftype='bessel', output='zpk')
			filter_roots = np.sort(p)
		elif self.ftype == 'cheby1':
			(z,p,k) = signal.iirfilter(self.order, psi, rp=1, btype='lowpass', analog=1, ftype='cheby1', output='zpk')
			filter_roots = np.sort(p)
		elif self.ftype == 'cheby2':
			(z,p,k) = signal.iirfilter(self.order, psi, rs=self.rs, btype='lowpass', analog=1, ftype='cheby2', output='zpk')
			filter_roots = np.sort(p)
		else:
			print 'Invalid filter type.'
			return
		# Separate the real pole from the complex poles
		real_index = 0;
		for i in range(self.order):
			if abs(filter_roots[i].imag) <= 1e-16:
				real_index = i
				break
		complex_roots = np.concatenate((filter_roots[0:real_index],filter_roots[real_index+1:]))
		# Put in the real pole's gamma value
		h_B = -1.0*filter_roots[real_index].real
		gammas.append((1.0 - h_B) / (1.0 + h_B))
		# Calculate coefficients of the individual Hurwitz polynomials
		for i in (range((order-1)/2)):
			h_A = -2.0*(complex_roots[2*i].real)
			h_B = abs(complex_roots[2*i])**2
			gammas.append((h_A - h_B - 1.0)/(h_A + h_B + 1.0))
			gammas.append((1.0 - h_B) / (1.0 + h_B))
		# Construct filter
		for i in range(self.order):
			self.adaptors.append(Adaptor(gammas[i],self.coeff_scale,self.dtype))
			self.registers.append(0)

    def __engine(self):
        """
            Engine of the function    
        """
        
        # ---------------------------------------------------------------------
        # Generate the base signal
        self.nSmp = round(self.fR * self.tS)  # The number of samples in the output signal
        self.mSig = np.random.randn(self.nSigs, self.nSmp)   # Generate the noise

        # ---------------------------------------------------------------------
        # Filter the signal with a low pass filter, if it is needed
        if self.wasParamGiven('fMax') and self.wasParamGiven('fMin'):

            # Design a iir low pass filter
            iCFP_l = self.fMin/(0.5*self.fR)   # Compute the filter parameter for the low cutoff frequency 
            iCFP_h = self.fMax/(0.5*self.fR)   # Compute the filter parameter for the high cutoff frequency 

            (vN, vD) = scsig.iirfilter(self.nFiltOrd, [iCFP_l, iCFP_h], btype='bandpass', ftype=self.strFilt,
                                       rs=self.iRs, rp=self.iRp)

            # Apply the filter
            self.mSig = scsig.lfilter(vN, vD, self.mSig)
            
        elif self.wasParamGiven('fMax'):

            # Design a iir low pass filter
            iCFP = self.fMax/(0.5*self.fR)   # Compute the filter parameter for the cutoff frequency 
            (vN, vD) = scsig.iirfilter(self.nFiltOrd, iCFP, btype='lowpass', ftype=self.strFilt,
                                       rs=self.iRs, rp=self.iRp)

            # Apply the filter
            self.mSig = scsig.lfilter(vN, vD, self.mSig)

        elif self.wasParamGiven('fMin'):

            # Design a iir low pass filter
            iCFP = self.fMin/(0.5*self.fR)   # Compute the filter parameter for the cutoff frequency 
            (vN, vD) = scsig.iirfilter(self.nFiltOrd, iCFP, btype='highpass', ftype=self.strFilt,
                                       rs=self.iRs, rp=self.iRp)

            # Apply the filter
            self.mSig = scsig.lfilter(vN, vD, self.mSig)

        # ---------------------------------------------------------------------

        # Adjust the signal power
        (self.mSig, self.vP) = self._adjPower(self.mSig, self.iP)
        
        return

    def _get_a(self):
        # Combine value/units from GUI to form actual cutoff frequency.
        val = self.filter_cutoff
        units = self.filter_cutoff_units
        if (units == "Hz"):
            fc = val
        elif (units == "kHz"):
            fc = val * 1e3
        elif (units == "MHz"):
            fc = val * 1e6
        else:
            fc = val * 1e9

        # Generate the filter coefficients.
        if (self.filter_type == "FIR"):
            w = fc/(self.sample_rate/2)
            b = firwin(self.Ntaps, w)
            a = [1]
        elif (self.filter_type == "IIR"):
            (b, a) = iirfilter(self.Ntaps - 1, fc/(self.sample_rate/2), btype='lowpass')
        else:
            a = self.usr_a[0]
            b = self.usr_b[0]
        if (self.filter_type != "custom"):
            self.a_str = reduce(lambda string, item: string + "%+06.3f  "%item, a, "")
            self.b_str = reduce(lambda string, item: string + "%+06.3f  "%item, b, "")
        self.b = b
        return a

    def __init__(self, fp, fs, gpass, gstop, ftype, btype): 
        #Variables init.
        self.fp = fp
        self.fs = fs
        self.gpass =  gpass
        self.gstop = gstop
        self.ftype = ftype
        self.btype = btype

        #Filter type for plot's title.
        types_dict = {"butter":"Butterworth", "cheby1":"Chebyshev I", "cheby2":"Chebyshev II", "ellip": "Cauer"}
        self.ftype_plot = types_dict[ftype]

        self.__wsk()
        self.__filter_order()
            
        #Designing filter.
        (self.b, self.a) = signal.iirfilter(self.ord,
                                            self.wn,
                                            rp=self.gpass, 
                                            rs=self.gstop, 
                                            btype=self.btype, 
                                            analog=True, 
                                            output='ba', 
                                            ftype=ftype)

        #Frequency response of analog filter.
        (self.w, self.h) = signal.freqs(self.b, self.a, worN=1000)

        #Denormalizing variabels for ploting. Pulsation to frequency.
        self.w = (self.w * (self.sampling_w / 2)) / (2 * pi)
        self.wn = (self.wn * (self.sampling_w / 2)) / (2 * pi)

def noise(t, amplitude, center_frequency, bandwidth):
    noise = np.random.uniform(low=-1, high=1, size=len(t))
    fl = center_frequency-bandwidth
    fh = center_frequency+bandwidth
    Wn = np.divide([fl, fh], fs/2.)
    b, a = iirfilter(8, Wn)
    return filtfilt(b, a, amplitude*noise)

  def test_slowsphere(self):
    '''
    When the input is stationary, the result should be similar to the
    symmetrical whitening transform but no low-frequency changes should
    exist.
    '''
    # init
    S = np.random.randn(1e4, 4)
    A = np.where(np.eye(4), 1, np.random.rand(4, 4))
    xs = np.dot(S, A)
    xs = xs-np.mean(xs, axis=0)

    # perform slow sphering
    xs2 = slow_sphere(xs, signal.iirfilter(2, .1, btype='low'), 10)

    # test slowness property
    wstep = 20
    sigs = np.asarray([cov0(xs2[i:i+wstep]) for i in 
      range(0, xs2.shape[0], wstep)])

    spec = np.log(np.abs(np.apply_along_axis(
      lambda x: np.mean(np.abs(stft(x, 256, 256)), axis=0), 0, sigs)))

    self.assert_(np.mean(spec[1:10]) < np.mean(spec[10:]))

    # test whiteness property 
    np.testing.assert_almost_equal(
      np.cov(xs2, rowvar=False), np.eye(4), decimal=1)
  
    # test that whitening preserves channel mapping
    np.testing.assert_equal(
      np.argmax(np.corrcoef(S.T, xs2.T)[:4, 4:], axis=1),
      np.arange(4))

	def compute_filter(self):
		freqRatio = self.frequency / self.nyquist
		print self.frequency, self.nyquist, freqRatio
		print type(self.frequency), type(self.nyquist), type(freqRatio)

		wn = freqRatio
		r = 0.1
		B, A = np.zeros(3), np.zeros(3)
		A[0],A[1],A[2] = 1.0, -2.0*r*np.cos(2*np.pi*wn), r*r
		B[0],B[1],B[2] = 1.0, -2.0*np.cos(2*np.pi*wn), 1.0


		self._filter_b = B
		self._filter_a = A

		#%Compute zeros
		#zeros = np.array([np.exp(0+1j *np.pi*freqRatio ), np.exp( 0-1j*np.pi*freqRatio )])

		#%Compute poles
		#poles = (1-self.notchWidth) * zeros

		#self._filter_b = np.poly( zeros ) # Get moving average filter coefficients
		#self._filter_a = np.poly( poles ) # Get autoregressive filter coefficients

		self._filter_b, self._filter_a = iirfilter(4, (45. / self.nyquist, 55. / self.nyquist), 'bandstop')

def lowpass(freq, df, corners=4):
    """
    Butterworth-Lowpass Filter.

    Filter data removing data over certain frequency ``freq`` using ``corners``
    corners.
    The filter uses :func:`scipy.signal.iirfilter` (for design)
    and :func:`scipy.signal.sosfilt` (for applying the filter).

    :type data: numpy.ndarray
    :param data: Data to filter.
    :param freq: Filter corner frequency.
    :param df: Sampling rate in Hz.
    :param corners: Filter corners / order.
    :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.
    :return: Filtered data.
    """
    fe = 0.5 * df
    f = freq / fe
    # raise for some bad scenarios
    if f > 1:
        f = 1.0
        msg = "Selected corner frequency is above Nyquist. " + \
              "Setting Nyquist as high corner."
        warnings.warn(msg)
    z, p, k = iirfilter(corners, f, btype='lowpass', ftype='butter',
                        output='zpk')
    sos = zpk2sos(z, p, k)
    return sos

    def update_plots(self, fs, data):
        self.current_update += 1
        data = signal.detrend(data.ravel())

        # Plot RMS
        if self._coefs is None:
            self._coefs = signal.iirfilter(2, (400.0/(fs/2), 40e3/(fs/2)))
            b, a = self._coefs
            self._zf = signal.lfiltic(b, a, data[:len(a)-1], data[:len(b)-1])
        b, a = self._coefs

        data, self._zf = signal.lfilter(b, a, data, zi=self._zf)
        rms = np.mean(data**2)**0.5
        db_rms = db(rms)-self.paradigm.mic_sens_dbv-db(20e-6)
        self.append_data(time=self.current_time, rms=db_rms)
        self.current_time += len(data)/fs

        self.current_spl = db_rms
        self.current_spl_average = self.rms_data.get_data('rms')[-60:].mean()
        self.overall_spl_average = self.rms_data.get_data('rms').mean()

        w_frequency = psd_freq(data, fs)
        w_psd = psd(data, fs, 'hamming')
        w_psd_db = db(w_psd)-self.paradigm.mic_sens_dbv-db(20e-6)
        self.rms_data.update_data(frequency=w_frequency, psd=w_psd_db)

 def __init__(self):
     self.filter_parameters = filter_parameters()
     # self.num, self.denom = iirdesign(wp=[self.filter_parameters.lowcut, self.filter_parameters.highcut], ws=[0.005,0.35], gpass=self.filter_parameters.mloss, gstop=self.filter_parameters.matten, analog=False, ftype='butter', output='ba')
     # self.num, self.denom = iirdesign(wp=[self.filter_parameters.lowcut, self.filter_parameters.highcut], ws=[0.01,0.35], gpass=0.05, gstop=5, analog=False, ftype='cheby1', output='ba')
     # self.num, self.denom = iirdesign(wp=[self.filter_parameters.lowcut, self.filter_parameters.highcut], ws=[0.005,0.35], gpass=1.7, gstop=20, analog=False, ftype='cheby2', output='ba')
     # self.num, self.denom = iirdesign(wp=[self.filter_parameters.lowcut, self.filter_parameters.highcut], ws=[0.005,0.35], gpass=0.4, gstop=28, analog=False, ftype='ellip', output='ba')
     self.num, self.denom = iirfilter(5, [self.filter_parameters.lowcut, self.filter_parameters.highcut], btype='bandpass', analog=False, ftype='bessel', output='ba')

    def run(self):
        out = 0
        z1 = 0
        z2 = 0
        outputlist = []
        for line in self.signal:
            z1 = line + z1 - z2
            out = self.quantization(z1)
            z2 = out
            outputlist.append(out)
        spectrum = numpy.fft.fft(outputlist[0:0+self.fftsample])
        spplot = (numpy.sqrt(numpy.power(spectrum.real,2)+numpy.power(spectrum.imag,2)))
        freqlist = numpy.fft.fftfreq(self.fftsample,d=1/self.samplerate)
        b, a = signal.iirfilter(4, 1000 / (self.samplerate / 2), btype = 'lowpass', analog = False, ftype = 'butter', output = 'ba')
        w,h = signal.freqz(b,a,self.fftsample)
        bpltlist = spplot
        #bpltlist *= abs(h)
        print len(spplot)
        print spplot
        plt.plot(freqlist[0:self.fftsample/2],spplot[0:self.fftsample/2])
        plt.axis([0, self.samplerate/2.0, 0, 500])
        plt.xlim()

        plt.xlabel("frequency [Hz]")
        plt.ylabel("amplitude spectrum")
        plt.show()