def freqz_(sys, w, dt=8e-9):
"""
This function computes the frequency response of a zpk system at an
array of frequencies.
It loosely mimicks 'scipy.signal.frequresp'.
Parameters
----------
system: (zeros, poles, k)
zeros and poles both in rad/s, k is the actual coefficient, not DC gain
w: np.array
frequencies in rad/s
dt: sampling time
Returns
-------
np.array(..., dtype=np.complex) with the response
"""
z, p, k = sys
b, a = sig.zpk2tf(z, p, k)
_, h = sig.freqz(b, a, worN=w*dt)
return h
python类freqz()的实例源码
def frequencyResponse(self, freqs=None):
if self.bandType == 'allpass':
return spsig.freqz(1, worN=freqs)
if self.bandType == 'allstop':
return spsig.freqz(0, worN=freqs)
numCoef = self.numCoef
denomCoef = self.denomCoef
if self.zeroPhase:
# http://www.mathworks.com/matlabcentral/newsreader/view_thread/245017
numCoef = np.convolve(numCoef,numCoef[::-1])
denomCoef = np.convolve(denomCoef,denomCoef[::-1])
# freqz does not preserve dtype of arguments, report bug XXX - idfah
return spsig.freqz(numCoef, denomCoef, worN=freqs)
def another():
# plot the figure of signals
fig = plt.figure(1)
ax = fig.add_subplot(311)
ax.set_title('input signal')
ax.plot(t[1:Npts],x[1:Npts]) # just plot part of the signal
ax = fig.add_subplot(312)
ax.set_title('expected signal')
ax.plot(t[1:Npts],xo[1:Npts])
ax = fig.add_subplot(313)
ax.set_title('filtered signal')
ax.plot(t[1:Npts],y[1:Npts])
fig.savefig('pic1.png')
# plot the mag & phase response of the LPF
w,h = signal.freqz(b,1)
hdb = 20 * np.log10(np.abs(h))
hphs = np.unwrap(np.arctan2(np.imag(h),np.real(h)))
fig2=plt.figure(2)
ax2 = fig2.add_subplot(211)
ax2.set_title('frequency response')
ax2.plot(w/max(w),hdb)
ax2 = fig2.add_subplot(212)
ax2.set_title('phase response')
ax2.plot(w/max(w),hphs)
fig2.savefig('pic2.png')
def test_absorption_filter(self):
b = uwa.absorption_filter(200000)
w, h = sp.freqz(b, 1, 4)
h = 20*np.log10(np.abs(h))
self.assertEqual(list(np.round(h)), [0.0, -2.0, -8.0, -17.0])
def freqz(b, a=1, fs=2.0, worN=None, whole=False):
"""Plot frequency response of a filter.
This is a convenience function to plot frequency response, and internally uses
:func:`scipy.signal.freqz` to estimate the response. For further details, see the
documentation for :func:`scipy.signal.freqz`.
:param b: numerator of a linear filter
:param a: denominator of a linear filter
:param fs: sampling rate in Hz (optional, normalized frequency if not specified)
:param worN: see :func:`scipy.signal.freqz`
:param whole: see :func:`scipy.signal.freqz`
:returns: (frequency vector, frequency response vector)
>>> import arlpy
>>> arlpy.signal.freqz([1,1,1,1,1], fs=120000);
>>> w, h = arlpy.signal.freqz([1,1,1,1,1], fs=120000)
"""
import matplotlib.pyplot as plt
w, h = _sig.freqz(b, a, worN, whole)
f = w*fs/(2*_np.pi)
fig = plt.figure()
ax1 = fig.add_subplot(111)
plt.plot(f, 20*_np.log10(abs(h)), 'b')
plt.ylabel('Amplitude [dB]', color='b')
plt.xlabel('Frequency [Hz]')
plt.grid()
ax1.twinx()
angles = _np.unwrap(_np.angle(h))
plt.plot(f, angles, 'g')
plt.ylabel('Angle (radians)', color='g')
plt.axis('tight')
plt.show()
return w, h
def plotfilt(b, fs:int, ofn=None):
if fs is None:
fs=1 #normalized freq
L = b.size
fg,axs = subplots(2,1,sharex=False)
freq, response = signal.freqz(b)
response_dB = 20*np.log10(abs(response))
if response_dB.max()>0:
logging.error('filter may be unstable')
axs[0].plot(freq*fs/(2*np.pi), response_dB)
axs[0].set_title(f'filter response {L} taps')
axs[0].set_ylim((-100,None))
axs[0].set_ylabel('|H| [db]')
axs[0].set_xlabel('frequency [Hz]')
t = np.arange(0, L/fs, 1/fs)
axs[1].plot(t,b)
axs[1].set_xlabel ('time [sec]')
axs[1].set_title('impulse response')
axs[1].set_ylabel('amplitude')
axs[1].autoscale(True,tight=True)
fg.tight_layout()
if ofn:
ofn = Path(ofn).expanduser()
ofn = ofn.with_suffix('.png')
print('writing',ofn)
fg.savefig(str(ofn),dpi=100,bbox_inches='tight')
def frequencyResponse(self, freqs=None):
return spsig.freqz(self.impulseResponse, worN=freqs)
def sosfreqz(sos, ws = None):
if ws is None:
H = [1] * 512
else:
H = [1] * len(ws)
for i in range(len(sos)):
w, h = freqz(sos[i,:3], sos[i, 3:], worN = ws)
H *= h
return w, H
def MOEar(self, correctionType = 'ELC'):
""" Method to approximate middle-outer ear transfer function for linearly scaled
frequency representations, using an FIR approximation of order 600 taps.
As appears in :
- A. Härmä, and K. Palomäki, ''HUTear – a free Matlab toolbox for modeling of human hearing'',
in Proceedings of the Matlab DSP Conference, pp 96-99, Espoo, Finland 1999.
Arguments :
correctionType : (string) String which specifies the type of correction :
'ELC' - Equal Loudness Curves at 60 dB (default)
'MAP' - Minimum Audible Pressure at ear canal
'MAF' - Minimum Audible Field
Returns :
LTq : (ndarray) 1D Array containing the transfer function, without the DC sub-band.
"""
# Parameters
firOrd = self.nfft
Cr, fr, Crdb = self.OutMidCorrection(correctionType, firOrd, self.fs)
Cr[self.nfft - 1] = 0.
# FIR Design
A = firwin2(firOrd, fr, Cr, nyq = self.fs/2)
B = 1
_, LTq = freqz(A, B, firOrd, self.fs)
LTq = 20. * np.log10(np.abs(LTq))
LTq -= max(LTq)
return LTq[:self.nfft/2 + 1]
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_attenuation(h, freq, gain):
"""Compute minimum attenuation at stop frequency"""
from scipy.signal import freqz
_, filt_resp = freqz(h.ravel(), worN=np.pi * freq)
filt_resp = np.abs(filt_resp) # use amplitude response
filt_resp /= np.max(filt_resp)
filt_resp[np.where(gain == 1)] = 0
idx = np.argmax(filt_resp)
att_db = -20 * np.log10(filt_resp[idx])
att_freq = freq[idx]
return att_db, att_freq
def tf_coefficients(self, frequencies=None, coefficients=None,
delay=False):
"""
computes implemented transfer function - assuming no delay and
infinite precision (actually floating-point precision)
Returns the discrete transfer function realized by coefficients at
frequencies.
Parameters
----------
coefficients: np.array
coefficients as returned from iir module
frequencies: np.array
frequencies to compute the transfer function for
dt: float
discrete sampling time (seconds)
zoh: bool
If true, zero-order hold implementation is assumed. Otherwise,
the delay is expected to depend on the index of biquad.
Returns
-------
np.array(..., dtype=np.complex)
"""
if frequencies is None:
frequencies = self.frequencies
frequencies = np.asarray(frequencies, dtype=np.float64)
if coefficients is None:
fcoefficients = self.coefficients
else:
fcoefficients = coefficients
# discrete frequency
w = frequencies * 2 * np.pi * self.dt * self.loops
# the higher stages have progressively more delay to the output
if delay:
delay_per_cycle = np.exp(-1j * self.dt * frequencies * 2 * np.pi)
h = np.zeros(len(w), dtype=np.complex128)
for i in range(len(fcoefficients)):
sos = np.asarray(fcoefficients[i], dtype=np.float64)
# later we can use sig.sosfreqz (very recent function, dont want
# to update scipy now)
ww, hh = sig.freqz(sos[:3], sos[3:], worN=np.asarray(w,
dtype=np.float64))
if delay:
hh *= delay_per_cycle ** i
h += hh
return h
def test_temporal_learning(flen=[5]):
'''
Function that computes CSP filters then uses those with the temporal filter MTL
idea, and confirms that the output has a spectral profile that is similar to expected.
Generate y values from trace of filtered covariance to ensure that is not an issue
'''
def covmat(x,y):
return (1/(x.shape[1]-1)*x.dot(y.T))
D = scio.loadmat('/is/ei/vjayaram/code/python/pyMTL_MunichMI.mat')
data = D['T'].ravel()
labels = D['Y'].ravel()
fig = plt.figure()
fig.suptitle('Recovered frequency filters for various filter lengths')
model = TemporalBRC(max_prior_iter=100)
# compute cross-covariance matrices and generate X
for find,freq in enumerate(flen):
X = []
for tind,d in enumerate(data):
d = d.transpose((2,0,1))
x = np.zeros((d.shape[0], freq))
nsamples = d.shape[2]
for ind, t in enumerate(d):
for j in range(freq):
C = covmat(t[:,0:(nsamples-j)],t[:,j:])
x[ind,j] = np.trace(C + C.T)
X.append(x)
labels[tind] = labels[tind].ravel()
model.fit_multi_task(X,labels)
b = numerics.solve_fir_coef(model.prior.mu.flatten())[0]
# look at filter properties
w,h = freqz(b[1:],worN=100)
w = w*500/2/np.pi # convert to Hz
ax1 = fig.add_subplot(3,3,find+1)
plt.plot(w, 20 * np.log10(abs(h)), 'b')
plt.ylabel('Amplitude [dB]', color='b')
plt.xlabel('Frequency [Hz]')
ax2 = ax1.twinx()
angles = np.unwrap(np.angle(h))
plt.plot(w, angles, 'g')
plt.ylabel('Angle (radians)', color='g')
plt.grid()
plt.axis('tight')
plt.show()
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
