This page describes how to perform some basic sound processing functions in python. We'll be using the pylab interface, which gives access to numpy and matplotlib, both these packages need to be installed. We'll also use the Audiolab package (formerly Pyaudiolab) to import wav files. Download and installation instructions are given on the author's page.

We'll begin by importing the necessary packages, assuming they've been already installed correctly 

  >>> from pylab import*
  >>> import pyaudiolab
please notice that the packages that are actually imported by the from pylab import* statement depend on the contents of your matplotlibrc file, which also determines whether interactive plotting will be available or not. For this tutorial I have set the following values on my matplotlibrc file (you can download it here on an Ubuntu install it goes into the '~/.matplotlib' directory) 
backend      : TkAgg
numerix      :  numpy #, Numeric or numarray
interactive  : True      # see http://matplotlib.sourceforge.net/interactive.html
Next we read in a wav file. You can download it here 440_sine.wav it contains a complex tone with a 440 Hz fundamental frequency (F0) plus noise. 
 
  >>> (snd, sampFreq, nBits) = pyaudiolab.wavread('440_sine.wav')
the wav file has two channels and 5060 sample points

  >>> snd.shape
  (5060, 2)
considering the sample rate (sampFreq = 44110) this corresponds to about 114 ms duration 
  >>> 5060.0 / sampFreq
  0.11473922902494331
we'll select and work only with one of the channels from now onwards

>> s1 = snd[:,0]

Unfortunately there is not an easy way of listening to the sound directly from python. If you're interested in having sound playback from python however check out the pygame library and have also a look at the psychopy package that uses the pygame library for sound playback in psychology experiments.

Plotting the Tone

A time representation of the sound can be obtained by plotting the pressure values against the time axis. However we need to create an array containing the time points first: 

  >>> timeArray = arange(0, 5060.0, 1)
  >>> timeArray = timeArray / sampFreq
  >>> timeArray = timeArray * 1000  #scale to milliseconds
now we can plot the tone 
  >>> plot(timeArray, s1, color='k')
  >>> ylabel('Amplitude')
  >>> xlabel('Time (ms)')


440 Hz Tone, Time Plot

Plotting the Frequency Content

Another useful graphical representation is that of the frequency content of the tone. We can obtain the frequency content of the sound using the fft function, that implements a Fast Fourier Transform algorithm. We'll follow closely the following technical document http://www.mathworks.com/support/tech-notes/1700/1702.html to obtain the power spectrum of our sound. 

 n = len(s1) 
 p = fft(s1) # take the fourier transform
notice that compared to the technical document, we didn't specify the number of points on which to take the fft, by default then the fft is computed on the number of points of the signal (n). Since we're not using a power of two the computation will be a bit slower, but for signals of this duration this is negligible. 
  nUniquePts = ceil((n+1)/2.0)
  p = p[0:nUniquePts]
  p = abs(p)
the fourier transform of the tone returned by the fft function contains both magnitude and phase information and is given in a complex representation (i.e. returns complex numbers). By taking the absolute value of the fourier transform we get the information about the magnitude of the frequency components. 
 
  p = p / float(n) # scale by the number of points so that
                   # the magnitude does not depend on the length 
                   # of the signal or on its sampling frequency  
  p = p**2  # square it to get the power 

 # multiply by two (see technical document for details)
 # odd nfft excludes Nyquist point
  if n % 2 > 0: # we've got odd number of points fft
      p[1:len(p)] = p[1:len(p)] * 2
  else:
      p[1:len(p) -1] = p[1:len(p) - 1] * 2 # we've got even number of points fft

  freqArray = arange(0, nUniquePts, 1.0) * (sampFreq / n);
  plot(freqArray/1000, 10*log10(p), color='k')
  xlabel('Frequency (kHz)')
  ylabel('Power (dB)')
The resulting plot can bee seen below, notice that we're plotting the power in decibels by taking 10*log10(p), we're also scaling the frequency array to kilohertz by dividing it by 1000

440 Hz Tone, Spectrum Plot

To confirm that the value we have computed is indeed the power of the signal, we'll also compute the root mean square (rms) of the signal. Loosely speaking the rms can be seen as a measure of the amplitude of a waveform. If you just took the average amplitude of a sinusoidal signal oscillating around zero, it would be zero since the negative parts would cancel out the positive parts. To get around this problem you can square the amplitude values before averaging, and then take the square root (notice that squaring also gives more weight to the extreme amplitude values): 
  >>> rms_val = sqrt(mean(s1**2))
  >>> rms_val
  0.0615000626299
since the rms is equal to the square root of the overall power of the signal, summing the power values calculated previously with the fft over all frequencies and taking the square root of this sum should give a very similar value 
  >>> sqrt(sum(p))
  0.0615000626299

References