Subversion Repositories pentevo

Rev

Blame | Last modification | View Log | Download | RSS feed

  1. #!/usr/bin/env python
  2.  
  3. import sys
  4. from sympy import *
  5. import math
  6. import cmath
  7. import random
  8.  
  9. #e24 = [1.0, 1.1, 1.2, 1.3, 1.5, 1.6, 1.8, 2.0, 2.2, 2.4, 2.7, 3.0, 3.3, 3.6, 3.9, 4.3, 4.7, 5.1, 5.6, 6.2, 6.8, 7.5, 8.2, 9.1];
  10. #e6  = [1.0,                1.5,                2.2,                3.3,                4.7,                6.8               ];
  11.  
  12.  
  13.  
  14. def play_mfb3_lowpass():
  15.  
  16.     # as in mfb3.png
  17.  
  18.     r1,r2,r3,r4,c1,c2,c3 = symbols('r1 r2 r3 r4 c1 c2 c3', real=True, positive=True)
  19.  
  20.     # angular frequency (omega)
  21.     w = symbols('w', real=True, positive=True)
  22.  
  23.     ho,hu = symbols('ho hu')
  24.     hon,hun = symbols('hon hun')
  25.  
  26.     ho = I*r3
  27.     hu = (c1*c2*c3*r1*r2*r3*r4*w**3 - I*c1*c3*r1*r2*r3*w**2 - I*c1*c3*r1*r2*r4*w**2 - I*c1*c3*r1*r3*r4*w**2 - c1*r1*r2*w - I*c2*c3*r1*r3*r4*w**2 - I*c2*c3*r2*r3*r4*w**2 - c3*r1*r3*w - c3*r1*r4*w - c3*r2*r3*w - c3*r2*r4*w - c3*r3*r4*w + I*r1 + I*r2)
  28.  
  29.     hu = expand(hu/(I*(r1+r2)))
  30.  
  31.     huc = hu.as_coefficients_dict(w)
  32.  
  33.     print('hu={}'.format(huc))
  34.  
  35.     ps0 = simplify(huc[w**0])
  36.     ps1 = simplify(huc[w**1])
  37.     ps2 = simplify(huc[w**2])
  38.     ps3 = simplify(huc[w**3])
  39.  
  40.     print('ps0={}\nps1={}\nps2={}\nps3={}'.format(ps0,ps1,ps2,ps3))
  41.  
  42.     #w1,w2,q = symbols('w1 w2 q')
  43.  
  44.     #weq1 = Eq( ps1, I*(1/w1 + 1/(w2*q)) )
  45.     #weq2 = Eq( ps2, (-1/(q*w1*w2)-1/(w2*w2)) )
  46.     #weq3 = Eq( ps3, -I/(w1*w2*w2) )
  47.  
  48.     #w_solve = solve( [weq1,weq2,weq3], [w1,w2,q] )
  49.  
  50.     #print(w_solve)
  51.  
  52. def play_mfb3_param():
  53.    
  54.     r1,r2,r3,r4,c1,c2,c3 = symbols('r1 r2 r3 r4 c1 c2 c3')
  55.     A,B,C = symbols('A B C')
  56.     w1,w2,q = symbols('w1 w2 q')
  57.     beta,gamma,kappa = symbols('beta gamma kappa')
  58.     #beta=1/w1 gamma=1/w2 kappa=1/q
  59.  
  60.     #eq1 = Eq( (c1*r1*r2 + c3*r1*r3 + c3*r1*r4 + c3*r2*r3 + c3*r2*r4 + c3*r3*r4)/(r1 + r2), (1/w1 + 1/(w2*q)) )
  61.     #eq2 = Eq( c3*(c1*r1*r2*r3 + c1*r1*r2*r4 + c1*r1*r3*r4 + c2*r1*r3*r4 + c2*r2*r3*r4)/(r1 + r2), (1/(q*w1*w2) + 1/(w2*w2)) )
  62.     #eq3 = Eq( c1*c2*c3*r1*r2*r3*r4/(r1 + r2), 1/(w1*w2*w2) )
  63.  
  64.     #eq1 = Eq( A, (1/w1 + 1/(w2*q)) )
  65.     #eq2 = Eq( B, (1/(q*w1*w2) + 1/(w2*w2)) )
  66.     #eq3 = Eq( C, 1/(w1*w2*w2) )
  67.  
  68.     #p_solve = solve( [eq1,eq2,eq3], [w1,w2,q], dict=True, manual=True, simplify=False )
  69.     #print(p_solve)
  70.  
  71.     eq1 = Eq( A, (beta + gamma*kappa) )
  72.     eq2 = Eq( B, (kappa*beta*gamma + gamma*gamma) )
  73.     eq3 = Eq( C, (beta*gamma*gamma) )
  74.     bgk = solve( [eq1,eq2,eq3], [beta,gamma,kappa], dict=True, manual=True, simplify=False, domain=S.Complexes )
  75.     print(bgk)
  76.  
  77.     for b in bgk:
  78.         beta_o  = b[beta]
  79.         gamma_o = b[gamma]
  80.         kappa_o = b[kappa]
  81.  
  82.         print('\n=====')
  83.         print('orig_beta={}'.format(beta_o))
  84.         print('orig_gamma={}'.format(gamma_o))
  85.         print('orig_kappa={}'.format(kappa_o))
  86.        
  87.         beta_s = simplify(beta_o)
  88.         print('simp_beta={}'.format(beta_s))
  89.         gamma_s = simplify(gamma_o)
  90.         print('simp_gamma={}'.format(gamma_s))
  91.         kappa_s = simplify(kappa_o)
  92.         print('simp_kappa={}'.format(kappa_s))
  93.  
  94.  
  95.  
  96. def main():
  97.  
  98.     #play_mfb3_lowpass()
  99.     play_mfb3_param()
  100.  
  101.  
  102.  
  103. if __name__=="__main__":
  104.     main()
  105.  
  106.  
  107.  
  108.