How to make a .wav file with Python, revisited

I've noticed that my two-year old post about making a .wav file with Python is being read a lot so I forced myself to modify the code a little. Now it should be nicer to read and modify further according to your needs/wishes. I've also incorporated the comments of helpful readers Mike Axiak and Fingon. Thanks guys ;-)
The link to the source of this idea, code by Andrea Valle, does not work as I'm typing this so I can only thank him by saying his name. My contribution to all this is minor.

import numpy as N
import wave

def get_signal_data(frequency=440, duration=1, volume=32768, samplerate=44100):
    """Outputs a numpy array of intensities"""
    samples = duration * samplerate
    period = samplerate / float(frequency)
    omega = N.pi * 2 / period
    t = N.arange(samples, dtype=N.float)
    y = volume * N.sin(t * omega)
    return y

def numpy2string(y):
    """Expects a numpy vector of numbers, outputs a string"""
    signal = "".join((wave.struct.pack('h', item) for item in y))
    # this formats data for wave library, 'h' means data are formatted
    # as short ints
    return signal

class SoundFile:
    def  __init__(self, signal, filename, duration=1, samplerate=44100):
        self.file = wave.open(filename, 'wb')
        self.signal = signal
        self.sr = samplerate
        self.duration = duration
  
    def write(self):
        self.file.setparams((1, 2, self.sr, self.sr*self.duration, 'NONE', 'noncompressed'))
        # setparams takes a tuple of:
        # nchannels, sampwidth, framerate, nframes, comptype, compname
        self.file.writeframes(self.signal)
        self.file.close()

if __name__ == '__main__':
    duration = 2
    myfilename = 'test.wav'
    data = get_signal_data(440, duration)
    signal = numpy2string(data)
    f = SoundFile(signal, myfilename, duration)
    f.write()
    print 'file written'

Creating NumPy arrays in Pyrex

In a previous post I have described basic usage of the Pyrex language which may be used to interface C code into Python or to speed up your Python code by adding C-like declarations. At the time I did not knew how to use numpy C-API to create new numpy arrays. I have mentioned it in the post and luckily for me a very helpful comment by diffusing thoughts showed me how to do this. It also got me started so that I could decipher how to use other C-API calls.
I currently do not have any production code in which the array creation is a bottleneck. Nevertheless I was curious if these C-API calls are actually faster than using the Python calls to NumPy. Here I present a howto describing what I learned together with some benchmarking results. I'm certain that I will forget how to use these C-API calls and hopefully this post will one day save me some time.

Prerequisites

• numpy and pyrex have to be installed

• c_numpy.pxd file from the numpy installation must be accessible (I just put it into my current working directory)
  

I use the following header of the pyrex (*.pyx) file:

import numpy
cimport c_numpy
cdef extern from "C:\Python25\Lib\site-packages\numpy\core\include\numpy\arrayobject.h":
    cdef object PyArray_SimpleNewFromData(int nd,
                                          c_numpy.npy_intp *dims,
                                          int typenum,
                                          void *data)
    cdef object PyArray_ZEROS(int nd,
                              c_numpy.npy_intp *dims,
                              int typenum,
                              int fortran)
    cdef object PyArray_SimpleNew(int nd,
                                  c_numpy.npy_intp *dims,
                                  int typenum)
    cdef object PyArray_Arange(double start,
                               double stop,
                               double step,
                               int typenum)

c_numpy.import_array()

The documentation for C-API of numpy is available for download on the numpy homepage in the form of "Guide to NumPy" pdf file (numpybook).

numpy.zeros vs. PyArray_ZEROS

example:

cdef int length
cdef c_numpy.ndarray newarr
length = 10
newarr = PyArray_ZEROS(1, &length, c_numpy.NPY_DOUBLE, 0)

notes:

• for multidimensional arrays the first two variables within PyArray_ZEROS have to be changed accordingly, see numpybook (I have not tested this)
  

• the type may also be changed if desirable (I however only need doubles)

equivalent numpy code:

newarr = numpy.zeros(length)

Benchmarking these two ways shows that they have the same speed for creating arrays larger than ~ 100 000 values. The C-API is faster on arrays smaller than ~ 50 000 values and about 50% faster on arrays of length 1 000.

numpy.arange vs. PyArray_Arange

example:

cdef double start, stop, step
cdef c_numpy.ndarray newarr
start = 0
stop = 10
step = 1
newarr = PyArray_Arange(start, stop, step, c_numpy.NPY_DOUBLE)

equivalent numpy code:

newarr = numpy.arange(start, stop, step)

Here the C-API is only faster on small arrays (length less than 1 000).

numpy.empty vs. PyArray_SimpleNew

example:

cdef int length
cdef c_numpy.ndarray newarr
length = 10
newarr = PyArray_SimpleNew(1, &length, c_numpy.NPY_DOUBLE)

equivalent numpy code:

newarr = numpy.empty(length)

This is the only case where using C-API is always faster than the numpy way. PyArray_SimpleNew is about 65% faster on arrays of length less than 50 000. It is ~20% faster on arrays of length 500 000. It is still somewhat faster in creating arrays of length 50 000 000.

PyArray_SimpleNewFromData
This call creates a new numpy array from malloc-ed C array.

example:

cdef extern from "stdlib.h":
    ctypedef int size_t
    void *malloc(size_t)

cdef c_numpy.npy_intp size
cdef c_numpy.ndarray newarr
cdef double *arrsource

size = 10
arrsource = <double *>malloc(sizeof(double) * size)
newarr = PyArray_SimpleNewFromData(1, &size, c_numpy.NPY_DOUBLE, <void *>arrsource)
newarr.flags = newarr.flags|(c_numpy.NPY_OWNDATA) # sets the ownership bit

Notes

• I have seen some discussion posts which discourage from doing this. See e.g. here.

• It may also be helpful to read this post by Travis Oliphant

Technical notes
The benchmarking was done by repeated allocation of a lot of arrays from within Pyrex code. The length of the arrays was also varied to assess whether there is some dependency on the amount of memory being allocated. The largest arrays tested were the largest which didn't caused MemoryError in Python. The number of allocations was selected to assure that the benchmark run at least several seconds (generally thousands of calls or more). I have also tested that it is safe to read and write values into the created arrays.
The benchmarking was done on a PC with Windows XP, Python 2.5.2, Pyrex 0.9.8.2, MinGW 5.1.4 and NumPy 1.2.0.

Pyrex for dummies

Pyrex is a very powerful tool if you need to speed up some Python code or connect it with some C code. My main motivation to use it was the ability to distribute only the compiled files and not also a C compiler etc.

Unfortunately, the learning curve is a bit steep if you want to do anything non-trivial (if you are not a skilled programmer, that is). I find it very difficult to read and understand just about anything about Python C API and related topics. So any attempt to learn new tricks means spending hours with Google and trying to find an example of usage which is simple enough for me to understand.

What follows is a few steps which I generally follow when rewriting Python into Pyrex code.

• Profile your code. Start with rewriting the slowest parts of your code (unless you want to know how fast you can make everything).


• Try to compile unchanged Python code with Pyrex. This way you will start with a working code and if something goes wrong you probably broke it with later changes.


• Declare loop control variables as C types:

  cdef int i

  and change Python for loops into Pyrex for loops

  for i from 0 <= i < 10:



• Numpy array data should be accessed directly. It's fairly easy to do now, that I finally know how to do this (see further).
  


• Additional speedup may be achieved by eliminating or minimising Python function calls and/or replacing them with C function calls.
  


• To help understanding your functions, add at least one line to the docstring such as:

  """myfunc( N, array x) -> array y"""

  Pyrex (or C) function parameters are not accessible as in Python functions by using help(myfunc) therefore you must explicitly write it to the docstring.
  

These steps have so far done the trick for my purposes. I essentially only need fast maths code so I have no idea about other areas. What I may have to learn later is some string stuff but so far I had neither guts nor reason to try it. This means I'm using pure Python strings in my Pyrex modules.

Now some selected details as promised:
To use some standard C function you have to declare it before use. So, if I e.g. need asin() from the math.h library I put this at the beginning of the Pyrex module:

cdef extern from "math.h":
   double asin(double)

Using Numpy C API to access array data directly was tough to learn, this is my current way:

• put to the beginning of the script:

  cimport c_numpy

  for this to work you have to copy 'c_numpy.pxd' file from 'Python\Lib\site-packages\numpy\doc\pyrex' into the directory with your script (there is a warning about future removal, I hope the same-named file in '\doc\cython' will work as well).


• initialize numpy by :

  c_numpy.import_array()



• declare new numpy array:

  cdef c_numpy.ndarray x



• create numpy array:

  x = numpy.zeros(10)

  There is another (possibly faster) way of creating new arrays but this is what i use now (I also do not know the other way, should have made a note...).


• declare pointer to store the adress of the numpy array data:

  cdef double *xdata



• copy data adress to your pointer:

  xdata = <double *>x.data

  The x.data is a char poiter to the first number in your array. I have no idea what this means (why char?).


• you may now index xdata to get desired element value:

  xdata[6] = 12.54

  or

  tempvar = xdata[1]



• you may declare the numpy array in the same way during function declaration:

  def myfunc(double step, c_numpy.ndarray x):
  ...

Using the Numpy C API is more cumbersome than just indexing numpy arrays but the code speedup is often significant.
Looking through my Pyrex modules this should in essence be all that is needed to get started with Pyrex.

How to make a simple .wav file with Python

EDIT: I've edited this code a little and incorporated reader's comments, here's the new version

I've recently had an idea to test some vibration-resonance problems by playing a set of sounds and detecting the vibrations in question. The glitch is, how to create a sound of desired frequency? I believe there is a lot of ways to solve this but having only a hammer (read Python) as my tool it looked like a nail ;-)
The idea is then to create a simple sine wave and play it back with a computer. My first try was with pygame module but the sound quality was poor and unsatisfactory.
If not playing sound directly one can write the wave to the file and then play the file with whatever means are suitable. This turned out to be the better way after all.
A search of the internet turned up this discussion of creating .wav files with Python. I've adopted the "minimal example" from there (written by Andrea Valle) to be usable with modern modules (it used obsolete Numeric and RandomArray modules). The code is not neat but it works well.

import numpy as N
import wave

class SoundFile:
   def  __init__(self, signal):
       self.file = wave.open('test.wav', 'wb')
       self.signal = signal
       self.sr = 44100

   def write(self):
       self.file.setparams((1, 2, self.sr, 44100*4, 'NONE', 'noncompressed'))
       self.file.writeframes(self.signal)
       self.file.close()

# let's prepare signal
duration = 4 # seconds
samplerate = 44100 # Hz
samples = duration*samplerate
frequency = 440 # Hz
period = samplerate / float(frequency) # in sample points
omega = N.pi * 2 / period

xaxis = N.arange(int(period),dtype = N.float) * omega
ydata = 16384 * N.sin(xaxis)

signal = N.resize(ydata, (samples,))

ssignal = ''
for i in range(len(signal)):
   ssignal += wave.struct.pack('h',signal[i]) # transform to binary

f = SoundFile(ssignal)
f.write()
print 'file written'

This code should work as intended and produce .wav files of given length and frequency. A few comments:
a) if you need different sampling frequencies, just replace all 44100 numbers with desired sampling frequency;
b) the setparams command is relatively well described in wave doc in Python help, its fourth parameter (nframes) may probably be whatever number you want, the procedure is designed to write all data you sent to it;
c) the 16384 number seems to be a volume setting for .wav file, this value is near the maximum so only lowering the volume is possible from here;
d) wave.struct.pack('h', #) - this was tough to find but without this command the resulting file will not be correct (at least on Windows machines), I do not yet understand what it does...