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 [81]:

    

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)**4, z ) 
Ao = sy.poly((z-a)**3, z )
Acl = Ac*Ao

# Assume anti-aliasing filter corresponding to a delay of two sample periods
A = sy.poly( (z**2*(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**4 - 2.22140275816017*z**3 + 1.22140275816017*z**2, z, domain='RR')
Poly(0.0214027581601699*z + 0.022877793471864, z, domain='RR')
Poly(z**3 - 3*a*z**2 + 3*a**2*z - a**3, z, domain='ZZ[a]')
Poly(1.0*z**4 - 1.79731585646889*z**3 + 1.21137910796793*z**2 - 0.36287181315765*z + 0.0407622039783662, z, domain='RR')

In [82]:

    

R = sy.poly(z**3 + r1*z**2+r2*z+r3, z)
S = sy.poly(s0*z**3 + s1*z**2 + s2*z + s3, z)
print A*R+B*S
print Acl

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


    

    




Poly(1.0*z**7 + (1.0*r1 - 2.22140275816017)*z**6 + (-2.22140275816017*r1 + 1.0*r2 + 1.22140275816017)*z**5 + (1.22140275816017*r1 - 2.22140275816017*r2 + 1.0*r3 + 0.0214027581601699*s0)*z**4 + (1.22140275816017*r2 - 2.22140275816017*r3 + 0.022877793471864*s0 + 0.0214027581601699*s1)*z**3 + (1.22140275816017*r3 + 0.022877793471864*s1 + 0.0214027581601699*s2)*z**2 + (0.022877793471864*s2 + 0.0214027581601699*s3)*z + 0.022877793471864*s3, z, domain='RR[r1,r2,r3,s0,s1,s2,s3]')
Poly(1.0*z**7 + (-3.0*a - 1.79731585646889)*z**6 + (3.0*a**2 + 5.39194756940666*a + 1.21137910796793)*z**5 + (-1.0*a**3 - 5.39194756940666*a**2 - 3.6341373239038*a - 0.36287181315765)*z**4 + (1.79731585646889*a**3 + 3.6341373239038*a**2 + 1.08861543947295*a + 0.0407622039783662)*z**3 + (-1.21137910796793*a**3 - 1.08861543947295*a**2 - 0.122286611935099*a)*z**2 + (0.36287181315765*a**3 + 0.122286611935099*a**2)*z - 0.0407622039783662*a**3, z, domain='RR[a]')

In [83]:

    

sol=sy.solve(dioph, (r1,r2, r3, s0,s1,s2, s3))
print sol


    

    




{s3: -1.78173668839598*a**3, s0: -0.869546949487873*a**3 + 56.6553501855205*a**2 - 70.0836348720752*a + 30.6100629635936, r3: -0.981389296931198*a**3 + 0.0596799465734188*a**2 - 1.29614940062047*a + 0.534452974786953, s2: 3.82121721355068e-29*a**2*(4.58706499600301e+29*a + 1.39882392611259e+29), s1: -16.9535027365645*a**3 - 55.7707162945703*a**2 + 63.853790915937*a - 28.5334483115476, r1: -3.0*a + 0.424086901691283, r2: 3.0*a**2 - 1.27226070507385*a + 0.932044162924379}

In [84]:

    

len(dioph)


    

    
Out[84]:





7

In [85]:

    

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


    

In [86]:

    

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


    

In [87]:

    

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


    

    




function [out1, out2, out3, out4, out5, out6, out7, out8, out9] = PS4_ship_heading_RST_example_aliasing(a)
  %PS4_SHIP_HEADING_RST_EXAMPLE_ALIASING  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 = -3.0*a + 0.424086901691283;
  out2 = 3.0*a.^2 - 1.27226070507385*a + 0.932044162924379;
  out3 = -0.981389296931198*a.^3 + 0.0596799465734188*a.^2 - 1.29614940062047*a + 0.534452974786953;
  out4 = -0.869546949487873*a.^3 + 56.6553501855205*a.^2 - 70.0836348720752*a + 30.6100629635936;
  out5 = -16.9535027365645*a.^3 - 55.7707162945703*a.^2 + 63.853790915937*a - 28.5334483115476;
  out6 = 3.82121721355068e-29*a.^2.*(4.58706499600301e+29*a + 1.39882392611259e+29);
  out7 = -1.78173668839598*a.^3;
  out8 = 2.07661465204603;
  out9 = a;

end

In [91]:

    

Ac.subs(z, 1)


    

    
Out[91]:





0.0919536423197624

In [92]:

    

(1-0.45)**4


    

    
Out[92]:





0.09150625000000003

In [94]:

    

B.subs(z,1)


    

    
Out[94]:





0.0442805516320339

In [ ]: