Commit 12153408 authored by hang yu's avatar hang yu
Browse files

more utilities supporting vectfit

parent 75185ca1
......@@ -45,7 +45,7 @@ def cc(z):
def model(s, poles, residues, d, h):
return sum(r/(s-p) for p, r in zip(poles, residues)) + d + s*h
def to_zpk(poles, residues, d, h):
def to_zpk(poles, residues, d, h, r_th=1e-8):
"""
HY: note that this vectfit returns in the partial fraction expansion form,
which may not be the most intuitive form.
......@@ -53,7 +53,43 @@ def to_zpk(poles, residues, d, h):
"""
bb, aa = sig.invres(residues, poles, d)
zz, pp, kk = sig.tf2zpk(bb, aa)
return zz, pp, kk
# sort and pair zeros and poles
z1=np.where(np.imag(zz)>r_th)
z1=zz[z1]
z2=np.where(np.abs(np.imag(zz)) <= r_th)
z2=np.sort(np.real(zz[z2]))
p1=np.where(np.imag(pp)>r_th)
p1=pp[p1]
p2=np.where(np.abs(np.imag(pp)) <= r_th)
p2=np.sort(np.real(pp[p2]))
np1, np2=len(p1), len(p2)
nz1, nz2=len(z1), len(z2)
idx_z=np.argsort(np.abs(z1))
z1 = z1[idx_z]
idx_p=np.argsort(np.abs(p1))
p1 = p1[idx_p]
zz = np.zeros(2 * len(z1) + len(z2), dtype=np.complex)
pp = np.zeros(2 * len(p1) + len(p2), dtype=np.complex)
for i in range (len(z1)):
zz[2*i] = z1[i]
zz[2*i+1] = np.conj(z1[i])
for i in range(len(p1)):
pp[2*i] = p1[i]
pp[2*i+1] = np.conj(p1[i])
zz[2*len(z1):]=z2
pp[2*len(p1):]=p2
return zz, pp, np.real(kk)
def vectfit_step(f, s, poles, ww=None):
"""
......@@ -257,24 +293,24 @@ def vectfit_auto_rescale(f, s, printparams=False, **kwargs):
# discard similar zpk pairs
def discard_similar_pz(zpk_s, cut=0.05, f_match=0.):
def discard_similar_pz(zpk_s, cut=0.05, f_match=0., r_th=1.e-8):
"""
If a pair of zeros is too close to a pair of poles, discard both of them.
"""
zz, pp, kk=zpk_s
z1=np.where(np.imag(zz)>1e-6)
z1=np.where(np.imag(zz)>r_th)
z1=zz[z1]
z2=np.where(np.abs(np.imag(zz)) <= 1e-6)
z2=zz[z2]
z2=np.where(np.abs(np.imag(zz)) <= r_th)
z2=np.sort(np.real(zz[z2]))
p1=np.where(np.imag(pp)>1e-6)
p1=np.where(np.imag(pp)>r_th)
p1=pp[p1]
p2=np.where(np.abs(np.imag(pp)) <= 1e-6)
p2=pp[p2]
p2=np.where(np.abs(np.imag(pp)) <= r_th)
p2=np.sort(np.real(pp[p2]))
np1, np2=len(p1), len(p2)
nz1, nz2=len(z1), len(z2)
......@@ -316,7 +352,7 @@ def discard_similar_pz(zpk_s, cut=0.05, f_match=0.):
return (zz, pp, np.real(kk))
def discard_features_low_coh(zpk_s, freq, coh,
coh_cut=0.5, f_match=0.):
coh_cut=0.5, f_match=0., r_th=1.e-8):
"""
for zeros/poles at frequencies where the coherence is below coh_cut
discard those zeros and poles
......@@ -331,17 +367,17 @@ def discard_features_low_coh(zpk_s, freq, coh,
zz, pp, kk=zpk_s
z1=np.where(np.imag(zz)>1e-6)
z1=np.where(np.imag(zz)>r_th)
z1=zz[z1]
z2=np.where(np.abs(np.imag(zz)) <= 1e-6)
z2=zz[z2]
z2=np.where(np.abs(np.imag(zz)) <= r_th)
z2=np.sort(np.real(zz[z2]))
p1=np.where(np.imag(pp)>1e-6)
p1=np.where(np.imag(pp)>r_th)
p1=pp[p1]
p2=np.where(np.abs(np.imag(pp)) <= 1e-6)
p2=pp[p2]
p2=np.where(np.abs(np.imag(pp)) <= r_th)
p2=np.sort(np.real(pp[p2]))
np1, np2=len(p1), len(p2)
nz1, nz2=len(z1), len(z2)
......@@ -405,6 +441,102 @@ def discard_features_low_coh(zpk_s, freq, coh,
return (zz, pp, np.real(kk))
def discard_hf_zeros(zpk_s, freq, f_cut, f_match=0., r_th=1.e-8):
"""
discard zeros above a certain freq f_cut
"""
w_cut = 2.*np.pi*f_cut
zz, pp, kk=zpk_s
z1=np.where(np.imag(zz)>r_th)
z1=zz[z1]
z2=np.where(np.abs(np.imag(zz)) <= r_th)
z2=np.sort(np.real(zz[z2]))
nz1, nz2=len(z1), len(z2)
z_del=[]
for i in range(nz1):
if np.abs(z1[i])>w_cut:
z_del.append(i)
kk *= (f_match-z1[i])*(f_match-np.conj(z1[i]))
z1=np.delete(z1, z_del)
idx_z=np.argsort(np.abs(z1))
z1 = z1[idx_z]
z_del=[]
for i in range(nz2):
if np.abs(z2[i])>w_cut:
z_del.append(i)
kk *= (f_match-z2[i])
z2=np.delete(z2, z_del)
idx_z=np.argsort(np.abs(z2))
z2 = z2[idx_z]
zz = np.zeros(2 * len(z1) + len(z2), dtype=np.complex)
for i in range (len(z1)):
zz[2*i] = z1[i]
zz[2*i+1] = np.conj(z1[i])
zz[2*len(z1):]=z2
return (zz, pp, np.real(kk))
def zpk_to_par_dict(zpk_s, r_th=1.e-8):
"""
convert a zpk in s domain to a par_dict used for sysID
"""
zz, pp, kk=zpk_s
z1=np.where(np.imag(zz)>r_th)
z1=zz[z1]
z2=np.where(np.abs(np.imag(zz)) <= r_th)
z2=np.sort(np.real(zz[z2]))
p1=np.where(np.imag(pp)>r_th)
p1=pp[p1]
p2=np.where(np.abs(np.imag(pp)) <= r_th)
p2=np.sort(np.real(pp[p2]))
np1, np2=len(p1), len(p2)
nz1, nz2=len(z1), len(z2)
zs_f, zs_Q = np.zeros(nz1), np.zeros(nz1)
ps_f, ps_Q = np.zeros(np1), np.zeros(np1)
for i in range(nz1):
zs_f[i], zs_Q[i] = get_res_f0_Q(np.real(z1[i]), np.imag(z1[i]))
for i in range(np1):
ps_f[i], ps_Q[i] = get_res_f0_Q(np.real(p1[i]), np.imag(p1[i]))
par_dict = {'zs_real':z2,
'ps_real':p2,
'zs_f':zs_f,
'zs_Q':zs_Q,
'ps_f':ps_f,
'ps_Q':ps_Q,
'ks':np.real(kk)}
return par_dict
def get_res_f0_Q(a, b):
"""
a, b are the real and imag parts of the complex pole
"""
w0 = np.sqrt(a**2 + b**2)
f0 = w0/2/np.pi
Q = -w0/2/a
return f0, Q
if __name__ == '__main__':
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment