Python Programming for Biology: Bioinformatics and Beyond

Figure 19.1.  Time series for various signal situations. A pure, single frequency signal

is shown in various combinations with an exponential decay function, random (normally

distributed) noise and a second pure frequency.

These considerations lead to the following Python code to create a simulated signal:

from numpy import arange, zeros, exp, pi, sqrt

from numpy.random import standard_normal as normal

I = 1j # Square root of -1; complex(0, 1)

def createSignal(parameters, tStep, nPoints, noise):

sig = zeros(nPoints, dtype=complex)

t = tStep * arange(nPoints, dtype=float)

for amplitude, frequency, decay in parameters:

sig += amplitude * exp(2*pi*I*frequency*t) * exp(-decay*t)

noise *= sqrt(0.5)

sig += noise*(normal(nPoints) + I*normal(nPoints))

return sig

The  function  starts  by  initialising  the  signal  array  sig  to  be  filled  with  zeros  for  the

correct  number  of  time  points.  The  NumPy  data  type  is  explicitly  set  to  be  complex

because  otherwise  it  would  be  regular  floating  point  numbers  and  the  sums  would  not

work  (the  imaginary  parts  would  be  dropped).  The  function  then  loops  over  each

component of the signal, each of which is specified by its amplitude, frequency and decay

from the parameters list passed to the function. Then a noise contribution is added, which

is  simply  sampling  a  normal  distribution  (via  normal())  for  both  real  and  imaginary

---

components  of  the  complex  number.  The  factor  of  sqrt(0.5)  in  the  noise  calculation  is

because  the  real  and  imaginary  components  both  contribute  to  the  standard  deviation.

Finally, the function passes back the NumPy array containing the signal.