我们从Python开源项目中,提取了以下22个代码示例,用于说明如何使用scipy.signal.iirfilter()。
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
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
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
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
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
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
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
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
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
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
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)
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)
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()
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
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
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