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
python类iirfilter()的实例源码
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 __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 __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 __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):
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 __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
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