In [1]:

    

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


    

In [12]:

    

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


    

In [25]:

    

pd1 = 0.4
pd2 = pd1
hh = 0.2
aa = 0
# Will need characteristic polynomial of order 9, since in the left hand side the order of 
# the polynomial will be n_A + n_R = 5+4 = 9. R and S will then have 4+5=9 unknowns 
Ac = sy.poly(z**3*(z - pd1)*(z-pd2), z ) 
#Ao = sy.poly(((z-a)**4).subs(a, aa), z )
Ao = sy.poly(((z-a)**4), z )
Acl = Ac*Ao

# Assume anti-aliasing filter corresponding to a delay of two sample periods + extra delay of 1 sample period
A = sy.poly( (z**3*(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 A
print B
print Ao
print Ac


    

    




Poly(1.0*z**5 - 2.22140275816017*z**4 + 1.22140275816017*z**3, z, domain='RR')
Poly(0.0214027581601699*z + 0.022877793471864, z, domain='RR')
Poly(z**4 - 4*a*z**3 + 6*a**2*z**2 - 4*a**3*z + a**4, z, domain='ZZ[a]')
Poly(1.0*z**5 - 0.8*z**4 + 0.16*z**3, z, domain='RR')

In [26]:

    

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

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


    

    




Poly(1.0*z**9 + (1.0*r1 - 2.22140275816017)*z**8 + (-2.22140275816017*r1 + 1.0*r2 + 1.22140275816017)*z**7 + (1.22140275816017*r1 - 2.22140275816017*r2 + 1.0*r3)*z**6 + (1.22140275816017*r2 - 2.22140275816017*r3 + 1.0*r4 + 0.0214027581601699*s0)*z**5 + (1.22140275816017*r3 - 2.22140275816017*r4 + 0.022877793471864*s0 + 0.0214027581601699*s1)*z**4 + (1.22140275816017*r4 + 0.022877793471864*s1 + 0.0214027581601699*s2)*z**3 + (0.022877793471864*s2 + 0.0214027581601699*s3)*z**2 + (0.022877793471864*s3 + 0.0214027581601699*s4)*z + 0.022877793471864*s4, z, domain='RR[r1,r2,r3,r4,s0,s1,s2,s3,s4]')
Poly(1.0*z**9 + (-4.0*a - 0.8)*z**8 + (6.0*a**2 + 3.2*a + 0.16)*z**7 + (-4.0*a**3 - 4.8*a**2 - 0.64*a)*z**6 + (1.0*a**4 + 3.2*a**3 + 0.96*a**2)*z**5 + (-0.8*a**4 - 0.64*a**3)*z**4 + 0.16*a**4*z**3, z, domain='RR[a]')

In [27]:

    

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


    

    




{s3: 0.0, s4: 0.0, r4: 0.455360670818637*a**4 - 2.32869913240506*a**3 + 5.38239533770142*a**2 - 5.12669902254084*a + 1.75143708540282, s0: 25.447156161159*a**4 - 156.844826966403*a**3 + 336.13593650405*a**2 - 306.22482671798*a + 101.635942011644, r3: -4.0*a**3 + 8.52841654896102*a**2 - 8.38442099709322*a + 2.92018873285642, s2: 0.0, s1: -17.3171761421304*a**4 + 124.324906890288*a**3 - 287.356056389879*a**2 + 273.704906641866*a - 93.5059619926153, r1: -4.0*a + 1.42140275816017, r2: 6.0*a**2 - 5.68561103264068*a + 2.0961052492733}

In [24]:

    

print len(dioph)
print dioph


    

    




11
[1.00000000000000, 1.0*r1 - 2.22140275816017, -2.22140275816017*r1 + 1.0*r2 + 0.22140275816017, 1.22140275816017*r1 - 2.22140275816017*r2 + 1.0*r3 + 0.8, 1.22140275816017*r2 - 2.22140275816017*r3 + 1.0*r4 + 0.0214027581601699*s0 - 0.16, 1.22140275816017*r3 - 2.22140275816017*r4 + 1.0*r5 + 0.022877793471864*s0 + 0.0214027581601699*s1, 1.22140275816017*r4 - 2.22140275816017*r5 + 0.022877793471864*s1 + 0.0214027581601699*s2, 1.22140275816017*r5 + 0.022877793471864*s2 + 0.0214027581601699*s3, 0.022877793471864*s3 + 0.0214027581601699*s4, 0.022877793471864*s4 + 0.0214027581601699*s5, 0.022877793471864*s5]

In [28]:

    

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


    

In [29]:

    

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


    

In [30]:

    

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


    

    




function [out1, out2, out3, out4, out5, out6, out7, out8, out9, out10, out11] = RST_example_aliasing(a)
  %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 = -4.0*a + 1.42140275816017;
  out2 = 6.0*a.^2 - 5.68561103264068*a + 2.0961052492733;
  out3 = -4.0*a.^3 + 8.52841654896102*a.^2 - 8.38442099709322*a + 2.92018873285642;
  out4 = 0.455360670818637*a.^4 - 2.32869913240506*a.^3 + 5.38239533770142*a.^2 - 5.12669902254084*a + 1.75143708540282;
  out5 = 25.447156161159*a.^4 - 156.844826966403*a.^3 + 336.13593650405*a.^2 - 306.22482671798*a + 101.635942011644;
  out6 = -17.3171761421304*a.^4 + 124.324906890288*a.^3 - 287.356056389879*a.^2 + 273.704906641866*a - 93.5059619926153;
  out7 = 0.0;
  out8 = 0.0;
  out9 = 0.0;
  out10 = 8.1299800190286;
  out11 = a;

end

In [20]:

    

sy.expand((z-a)**3)


    

    
Out[20]:





-a**3 + 3*a**2*z - 3*a*z**2 + z**3

In [ ]: