stereoscope.py

#!/usr/bin/env python

# a script to analyse a set of files to illustrate how the panning tends to vary by frequency, via a simple FFT analysis.

# REQUIREMENTS:
#  - Python (was tested on 2.7)
#  - Python modules:
#     - numpy
#     - scikits.audiolab
#     - matplotlib

# written by Dan Stowell 2014. public domain.
# modified by Daniel Jones with fixed FFT and wave for audio reading
# comments: http://mcld.co.uk/blog/blog.php?417


import sys
import wave, struct

import numpy as np
import os, glob

import matplotlib.pyplot as plt
import matplotlib.cm as cm

from scikits.audiolab import Sndfile
from scikits.audiolab import Format

USE_WAVE = True
USE_STEREO_FFT = True

#########################################
# config

# this is the pattern to search for files. 
globpattern = sys.argv[1] if len(sys.argv) > 1 else 'bette.1024.mono.wav'

framelen = 1024
fs = 44100
stereobins = 41  # odd number, we want to have a central bin
panlimit = 10.0  # +- extreme of the values we expect generally

# derived config
window = np.hamming(framelen)
window = np.vstack((window, window)).T
panmul = ((stereobins-1)/2) / panlimit
bintofreq = fs/framelen

#########################################
# functions

def analyse_file(audiopath):
	print "opening file %s" % audiopath
	if USE_WAVE:
		fd = wave.open(file(audiopath))
	else:
		sf = Sndfile(audiopath, "r")
		if (sf.channels != 2):
			raise ValueError("non-stereo file: %i channels." % sf.channels)
		if sf.samplerate != fs:
			raise ValueError("wanted sample rate %g - got %g." % (fs, sf.samplerate))

	histo = np.zeros((framelen/2, stereobins))	

	while(True):
		try:
			if USE_WAVE:
				frames = fd.readframes(framelen)
				if len(frames) < framelen * 4:
					break
				# unpack int16 to float32
				chunk = np.array([ struct.unpack("<hh", frames[n*4:(n+1)*4]) for n in range(1024) ]) / 32768.0
			else:
				# appears to give false results when tested with mono files...
				chunk = sf.read_frames(framelen, dtype=np.float32)
				if len(chunk) != framelen:
					print("Not read sufficient samples (%i) - returning." % (len(chunk)))
					break

			# uncomment to test with mono input by assigning right channel to left
			# chunk[:, 0] = chunk[:, 1]

			if USE_STEREO_FFT:
				# split our time series into two spectra and analyse separately
				specL = np.fft.fft(window[:, 0] * chunk[:, 0])
				powspecL = np.abs(specL[:framelen/2]) ** 2
				powspecL = np.maximum(1e-3, powspecL)
				specR = np.fft.fft(window[:, 1] * chunk[:, 1])
				powspecR = np.abs(specR[:framelen/2]) ** 2
				powspecR = np.maximum(1e-3, powspecR)
				panpos = np.clip(np.log(powspecR / powspecL), -panlimit, panlimit)
			else:
				# incorrect code for stereo FFT
				framespectrum = np.fft.fft(window * chunk)
				powspec = np.maximum(1e-3, abs(framespectrum[:framelen/2]) ** 2)
				panpos = np.clip(np.log(powspec[:, 0] / powspec[:, 1]), -panlimit, panlimit)

			histo[range(len(panpos)), map(int, (panpos + panlimit) * panmul)] += 1
		except RuntimeError:
			break

	if not USE_WAVE:
		sf.close()

	return histo

def plot_histo(histo, plotpath):
	histo = np.log(histo + 1) # log-scale the results
	plt.imshow(histo, origin='lower', interpolation='nearest', aspect='auto', cmap=cm.hot)
	plt.xlabel("Pan pos", fontsize='small')
	plt.ylabel("Frequency (kHz)", fontsize='small')
	plt.xticks([(stereobins-1)/2], ['0'], fontsize='small')
	yticks = np.arange(0, histo.shape[0], histo.shape[0]/10)
	plt.yticks(yticks, np.round(yticks * bintofreq * 0.001), fontsize='small')
	plt.title("Histogram of pan positions per frequency", fontsize='small')
	plt.savefig(plotpath)
	plt.clf()

#########################################
# Let's go
if __name__ == '__main__':
	bighisto = None
	print("Searching for files matching '%s'" % globpattern)
	for matched in glob.iglob(globpattern):
		try:
			histo = analyse_file(matched)
			print("Analysed: %s" % matched)
		except ValueError, e:
			print("Skipping a file due to ValueError (presumably wasn't stereo): %s, %s" % (matched, e))
			continue
		if bighisto == None:
			bighisto = histo
		else:
			bighisto += histo
	plot_histo(histo, "stereoscope.png")