Bodeplots in Python
===================


DIY Python
----------

Consider the (angular) frequency reponse function of a low-pass filter:

.. math::
   H(\omega) = \frac{1}{1+j \omega / \omega_c}

where $\omega_c$ is the cut-off frequency. In a Bode magnitude plot we
plot the magnitude (in decibels) of the transfer function (frequency
response), i.e.

.. math::
   20 \log | H(\omega) | = 20 \log \frac{1}{|1+j \omega / \omega_c|}\\
   = - 20 \log |1+j \omega / \omega_c|\\
   = - 10 \log\left( 1 + \frac{\omega^2}{\omega_c^2} \right)

Be sure you can do these steps yourself, especcially the last step is not trivial!

Most often in plots we plot 'real' frequencies and not angular
frequencies. Remember that $\omega = 2\pi f$. The cut-off frequency is
taken at 2KHz leading to $\omega_c = 4000\pi$. We may write a simple
python function to represent the transfer function:

.. ipython:: python

   from pylab import *

   def H(w):
       wc = 4000*pi
       return 1.0 / (1.0 + 1j * w / wc)

and then plot it:

.. ipython:: python

   f = logspace(1,5) # frequencies from 10**1 to 10**5
   @savefig bodeplot_lowpass.png height=8cm align=center
   plot(f, 20*log10(abs(H(2*pi*f)))); xscale('log')

Observe that the corner (or cut-off) frequency is at around 2000 kHz
and that the decay above that frequency is 20 dB per octave (as it
should be).


Scipy Signal Processing Package
-------------------------------

Scipy also contains functions to represent continuous time linear
systems. An LTI system is specified in the $s$-domain. For the
low-pass filter we have used in the previous section the transfer
function is:

.. math::
   H(s) = \frac{1}{\frac{1}{\omega_c}s + 1}

.. ipython:: python

   from pylab import *
   from scipy import signal

   system = signal.lti([1], [1/(4000*pi),1])
   f = logspace(1, 5)
   w = 2 * pi * f
   w, mag, phase = signal.bode(system,w)
   @savefig bode_scipy.png height=8cm align=center
   semilogx( f, mag);