pythonsyntaxsympy

Syntax improvement with working SymPy statement (novice level question)


My question may be about how to avoid putting an array into another unneeded array in SymPy. There may be still more to the question I'm not aware of, though. But within my limitations, that is the question I have at hand.

To make this question explicit and concrete, see the following...

I want to compute the magnitude of a complex-valued expression. In effect:

enter image description here

where H here is a complex expression and the asterisk means conjugation.

The actual computation I want to perform is:

enter image description here

I have a function in SymPy that computes the Butterworth terms. Here's some examples of what it produces (so that it is very clear what I'm starting with):

Butterworth(2)
[s**2 + sqrt(2)*s + 1]
Butterworth(3)
[s + 1, s**2 + s + 1]
Butterworth(4)
[s**2 + s*(sqrt(2)*sqrt(2 - sqrt(2)) + sqrt(2)*sqrt(sqrt(2) + 2) + 2*sqrt(sqrt(2) + 2))/4 + 1,
 s**2 - s*(-sqrt(2)*sqrt(sqrt(2) + 2) - 2*sqrt(2 - sqrt(2)) + sqrt(2)*sqrt(2 - sqrt(2)))/4 + 1]
Butterworth(6)
[s**2 + s*(sqrt(2) + sqrt(6))/2 + 1,
 s**2 + sqrt(2)*s + 1,
 s**2 - s*(-sqrt(6) + sqrt(2))/2 + 1]

The function creates an array of 1st and 2nd order factors. I combine these using prod() and expand() in order to compose the fuller expression. The complex variable s is then replaced with I*omega, where omega is declared to be a real, positive variable.

What I use right now is the following, in SymPy:

-20*ln(sqrt(expand(prod([[i,conjugate(i)] for i in [expand(prod(Butterworth(4))).subs(s,I*omega)]][0])).subs(omega,1.8)),10).n()
-20.4610324752877

And that is the correct answer for a low-pass 4th order Butterworth filter at a normalized frequency of 1.8. But to get there I had to use what I could figure out through some guesswork.

Let me break down my thinking in the above statement:

# expand out the product of the array of terms and replace complex `s` with `I*omega`
expand(prod(Butterworth(4))).subs(s,I*omega)

# put that into an array of just 1 element so I can apply a "for" statement
[expand(prod(Butterworth(4))).subs(s,I*omega)]

# construct an array of 2 elements, with the 2nd one being the conjugate
[[i,conjugate(i)] for i in [expand(prod(Butterworth(4))).subs(s,I*omega)]

# the unfortunate side effect of the above method is that it creates an array
# within an array so the following extracts out the array I want.
  [[i,conjugate(i)] for i in [expand(prod(Butterworth(4))).subs(s,I*omega)][0]
# ^^^^^^^^^^^^^^^^^                                                        ^^^
#      can I do                                      that would allow me to avoid
# something above here                               having to add this [0] here?

# ... the above may be where I want some help

# at this point I can just apply the prod(), expand(), sqrt(), etc:
-20*ln(sqrt(expand(prod([[i,conjugate(i)] for i in [expand(prod(Butterworth(4))).subs(s,I*omega)]][0])).subs(omega,1.8)),10).n()
-20.4610324752877

This question is about syntax. I'm in a learning-mode and trying to expand my knowledge of SymPy syntax and I would like to know if there is a simpler syntax for that single-line expression I've already worked out.


Solution

  • This is how I would approach it:

    from sympy import *
    init_printing()
    var("s")
    omega = symbols("omega", real=True, positive=True)
    
    # Butterworth(4)
    b4 = [
        s**2 + s*(sqrt(2)*sqrt(2 - sqrt(2)) + sqrt(2)*sqrt(sqrt(2) + 2) + 2*sqrt(sqrt(2) + 2))/4 + 1,
        s**2 - s*(-sqrt(2)*sqrt(sqrt(2) + 2) - 2*sqrt(2 - sqrt(2)) + sqrt(2)*sqrt(2 - sqrt(2)))/4 + 1
    ]
    
    H = prod(b4)
    -20 * log(sqrt(H * conjugate(H)).subs(s, I*omega).subs(omega, 1.8), 10).n()
    # -20.4610324752877
    

    EDIT to satisfy comment:

    By removing the outer list, you create an iterator. Then, you can call next to retrieve the first and only element.

    -20 * log(sqrt(prod(next([i,conjugate(i)] for i in [expand(prod(b4)).subs(s,I*1.8)]))), 10).n()