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
评论列表
文章目录
