In [1]:

    

import numpy as np
import sympy as sy
from sympy import I, exp
from sympy.utilities.codegen import codegen


    

In [2]:

    

z = sy.symbols("z", real=False)
r1,r2,r3, s0,s1, s2,s3 = sy.symbols("r1,r2,r3, s0,s1,s2,s3", real=True)
a, h = sy.symbols("a, h", real=True)


    

In [11]:

    

hh = 0.2
pcont = -4
pd1 = np.exp(hh*pcont)
#Ac = sy.poly(z**2*(z - pd1)*(z-pd2), z ) 
Ac = sy.poly((z - pd1)**2, z ) 
Ao = sy.poly((z-a), z )
Acl = Ac*Ao

# Assume anti-aliasing filter corresponding to a delay of two sample periods
A = sy.poly( ((z-1)*(z-exp(h))).subs(h,hh), z )
B = sy.poly( ((-1+exp(h)-h)*z + 1-(1-h)*exp(h)).subs(h,hh), z)
print pd1
print A
print B
print Ao
print Ac


    

    




0.449328964117
Poly(1.0*z**2 - 2.22140275816017*z + 1.22140275816017, z, domain='RR')
Poly(0.0214027581601699*z + 0.022877793471864, z, domain='RR')
Poly(z - a, z, domain='ZZ[a]')
Poly(1.0*z**2 - 0.898657928234443*z + 0.201896517994655, z, domain='RR')

In [15]:

    

R = sy.poly(z + r1, z)
S = sy.poly(s0*z + s1, z)
print A*R+B*S
print Acl

dioph=(A*R+B*S-Acl).all_coeffs()


    

    




Poly(1.0*z**3 + (1.0*r1 + 0.0214027581601699*s0 - 2.22140275816017)*z**2 + (-2.22140275816017*r1 + 0.022877793471864*s0 + 0.0214027581601699*s1 + 1.22140275816017)*z + 1.22140275816017*r1 + 0.022877793471864*s1, z, domain='RR[r1,s0,s1]')
Poly(1.0*z**3 + (-1.0*a - 0.898657928234443)*z**2 + (0.898657928234443*a + 0.201896517994655)*z - 0.201896517994655*a, z, domain='RR[a]')

In [16]:

    

sol=sy.solve(dioph, (r1, s0,s1))
print sol


    

    




{s1: 17.1460808200645*a - 29.3065846755898, s0: -23.9942021415893*a + 36.1547059971148, r1: -0.486457894317333*a + 0.548934401117437}

In [14]:

    

len(dioph)


    

    
Out[14]:





4

In [18]:

    

# Reorganize solution expression for matlab code generation
sol_expr = ('PS4_ship_heading_RST_example_noaliasing', [sol[r1], sol[s0], sol[s1], Ac.subs(z, 1)/B.subs(z,1), a])


    

In [19]:

    

# Export to matlab code
[(m_name, m_code)] = codegen(sol_expr, 'octave')


    

In [20]:

    

print m_code
with open("/home/kjartan/Dropbox/undervisning/tec/MR2007/matlab/PS4_ship_heading_RST_example_noaliasing.m", "w") as text_file:
    text_file.write(m_code)


    

    




function [out1, out2, out3, out4, out5] = PS4_ship_heading_RST_example_noaliasing(a)
  %PS4_SHIP_HEADING_RST_EXAMPLE_NOALIASING  Autogenerated by sympy
  %   Code generated with sympy 0.7.6
  %
  %   See http://www.sympy.org/ for more information.
  %
  %   This file is part of 'project'

  out1 = -0.486457894317333*a + 0.548934401117437;
  out2 = -23.9942021415893*a + 36.1547059971148;
  out3 = 17.1460808200645*a - 29.3065846755898;
  out4 = 6.84812132152482;
  out5 = a;

end

In [21]:

    

Ac.subs(z, 1)


    

    
Out[21]:





0.303238589760212

In [92]:

    

(1-0.45)**4


    

    
Out[92]:





0.09150625000000003

In [94]:

    

B.subs(z,1)


    

    
Out[94]:





0.0442805516320339

In [23]:

    

print dioph


    

    




[1.0*a + 1.0*r1 + 0.0214027581601699*s0 - 1.32274482992573, -0.898657928234443*a - 2.22140275816017*r1 + 0.022877793471864*s0 + 0.0214027581601699*s1 + 1.01950624016551, 0.201896517994655*a + 1.22140275816017*r1 + 0.022877793471864*s1]

In [24]:

    

0.55**2/(0.021+0.023)


    

    
Out[24]:





6.875000000000002

In [ ]: