algebraic expressions in python

Question

I'm working with some relatively simple algebraic expressions in python and wondered if there's a better way of plotting algebraic expressions than what I'm currently doing:

I have the equation of a notch filter (below image - left side from a book; the right side graph is the fig generated by my code), and so far my code works though is crude.

Is there a better way of plotting the magnitude? As shown the constants are; R = 50.0, C = 470e-12, L = 54e-6, and the desired frequency range is from 0 to 2MHz.

import matplotlib.pyplot as plt
import numpy as np
import math

R = 50.0
C = 470e-12
L = 54e-6

FREQ = []
DATA = []
for i in range(1, 200):
    f = i*10000.0
    w = 2*np.pi*f

    Ztop = w*L - 1.0/(w*C)
    Zbot = math.sqrt( (math.pow(R,2) + math.pow((w*L) -(1.0/(w*C)),2)) )
    Zout = abs(Ztop / Zbot)
    FREQ.append( f/1e6 )
    DATA.append( Zout )

plt.figure(1)
plt.plot(FREQ,DATA, '-k')
plt.xlabel('Frequency (MHz)')
plt.ylabel('Mag.')
plt.grid()
plt.show()

enter image description here

score 4 · Accepted Answer · answered Apr 28 '13 at 00:48

The main improvement I'd suggest is using Numpy's array broadcasting facility. Basically this will perform a mathematical operation on a whole array at once. Other than that, there's nothing really wrong with your code, though you could calculate the entire magnitude in one expression. That's your call, depending on whether you think it makes the code clearer.

# creates 201 values evenly spaced from 0 to 2e6, inclusive
f = np.linspace(0, 2e6, 201)
w = 2 * np.pi * f
# now calculate the magnitude for all 201 values at once
magnitude = (w * L - 1 / (w * C)) / np.sqrt(R**2 + (w * L - 1 / (w * C))**2)
# now you can plot magnitude vs. frequency
plt.plot(f / 1e6, magnitude, '-k')

Thanks David, yes 'array broadcasting' is far better - I didn't know python/numpy could do this. — Harry Lime, Apr 28 '13 at 00:52

score 3 · Answer 2 · answered Apr 28 '13 at 02:27

You could also use sympy to plot symbolic expressions:

from __future__ import division
from sympy import pi, pprint, sqrt, symbols # $ pip install sympy
from sympy.plotting import plot

f = symbols("f")
L, C, R = symbols("L, C, R", real=True)
w = 2*pi*f
W = w * L - 1 / (w * C)
mag = abs(W / sqrt(R**2 + W**2))
pprint(mag)
plot(mag.subs({R: 50.0, C: 470e-12, L: 54e-6}), (f, 1, 2e6))

Output

│                    1           │
│       2⋅π⋅L⋅f - ───────        │
│                 2⋅π⋅C⋅f        │
│────────────────────────────────│
│     ___________________________│
│    ╱                         2 │
│   ╱   2   ⎛             1   ⎞  │
│  ╱   R  + ⎜2⋅π⋅L⋅f - ───────⎟  │
│╲╱         ⎝          2⋅π⋅C⋅f⎠  │

magnitude

algebraic expressions in python

2 Answers2

Output

Linked