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)
wn, zeta, a,h,r1,r2,r3, s0,s1, s2,s3 = sy.symbols("wn, zeta,a,h,r1,r2,r3, s0,s1,s2,s3", real=True)


    

In [5]:

    

pc1 = -wn*zeta + I*wn*sy.sqrt(1-zeta**2)
pc2 = -wn*zeta - I*wn*sy.sqrt(1-zeta**2)
hn = 0.14
rd = exp(-wn*zeta*h)
argd = wn*sy.sqrt(1-zeta**2)*h
reald = rd*sy.cos(argd)
imd = rd*sy.sin(argd)
Ac = sy.poly( ((z**2 - 2*reald*z + rd**2)*z**2).subs(h,hn).subs(wn, 2*np.sqrt(2)).subs(zeta, 0.707), z ) 
ad = sy.exp(-a*h).subs(h,hn)
Ao = sy.poly( ((z-ad)**3).subs(h, hn), z )
Acl = Ac*Ao

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

print Ao
print Ac


    

    




Poly(z**3 - 3*exp(-0.14*a)*z**2 + 3*exp(-0.28*a)*z - exp(-0.42*a), z, domain='ZZ(exp(0.14*a))')
Poly(1.0*z**4 - 1.45274391357953*z**3 + 0.571257370977114*z**2, z, domain='RR')

In [4]:

    

print B
print A


    

    




Poly(-0.104216258544275*z + 0.138406534740473, z, domain='RR')
Poly(1.0*z**4 - 1.75578374145573*z**3 + 0.755783741455725*z**2, z, domain='RR')

In [6]:

    

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 - 1.75578374145573)*z**6 + (-1.75578374145573*r1 + 1.0*r2 + 0.755783741455725)*z**5 + (0.755783741455725*r1 - 1.75578374145573*r2 + 1.0*r3 - 0.104216258544275*s0)*z**4 + (0.755783741455725*r2 - 1.75578374145573*r3 + 0.138406534740473*s0 - 0.104216258544275*s1)*z**3 + (0.755783741455725*r3 + 0.138406534740473*s1 - 0.104216258544275*s2)*z**2 + (0.138406534740473*s2 - 0.104216258544275*s3)*z + 0.138406534740473*s3, z, domain='RR[r1,r2,r3,s0,s1,s2,s3]')
Poly(1.0*z**7 + 1.0*(-1.45274391357953*exp(0.14*a) - 3.0)*exp(-0.14*a)*z**6 + 1.0*(4.35823174073858*exp(0.14*a) + 0.571257370977114*exp(0.28*a) + 3.0)*exp(-0.28*a)*z**5 + 1.0*(-4.35823174073858*exp(0.14*a) - 1.71377211293134*exp(0.28*a) - 1.0)*exp(-0.42*a)*z**4 + 1.0*(1.71377211293134*exp(0.14*a) + 1.45274391357953)*exp(-0.42*a)*z**3 - 0.571257370977114*exp(-0.42*a)*z**2, z, domain='RR(exp(0.14*a))')

In [7]:

    

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


    

    




{s3: 0.0, s0: -11.1379297663082*exp(-0.42*a) + 35.6524726341564*exp(-0.28*a) - 37.3444331011921*exp(-0.14*a) + 12.8743964377876, r3: -2.16075336817355*exp(-0.42*a) + 4.62468678941254*exp(-0.28*a) - 4.93453519222313*exp(-0.14*a) + 1.72290452579443, s2: 0.0, s1: 7.67163845244122*exp(-0.42*a) - 25.253598692555*exp(-0.28*a) + 26.9455591595906*exp(-0.14*a) - 9.40810512392042, r1: -3.0*exp(-0.14*a) + 0.303039827876199, r2: 3.0*exp(-0.28*a) - 0.90911948362861*exp(-0.14*a) + 0.347546032319962}

In [8]:

    

# Reorganize solution expression for matlab code generation
sol_expr = ('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), ad])


    

In [9]:

    

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


    

In [10]:

    

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


    

    




function [out1, out2, out3, out4, out5, out6, out7, out8, out9] = 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 = -3.0*exp(-0.14*a) + 0.303039827876199;
  out2 = 3.0*exp(-0.28*a) - 0.90911948362861*exp(-0.14*a) + 0.347546032319962;
  out3 = -2.16075336817355*exp(-0.42*a) + 4.62468678941254*exp(-0.28*a) - 4.93453519222313*exp(-0.14*a) + 1.72290452579443;
  out4 = -11.1379297663082*exp(-0.42*a) + 35.6524726341564*exp(-0.28*a) - 37.3444331011921*exp(-0.14*a) + 12.8743964377876;
  out5 = 7.67163845244122*exp(-0.42*a) - 25.253598692555*exp(-0.28*a) + 26.9455591595906*exp(-0.14*a) - 9.40810512392042;
  out6 = 0.0;
  out7 = 0.0;
  out8 = 3.46629131386677;
  out9 = exp(-0.14*a);

end