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.
# fixed thanks to Daniel Jones - see blog and comments: http://mcld.co.uk/blog/blog.php?417


import numpy as np
import os, sys, glob

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

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

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

# this is the pattern to search for files. 
globpattern = sys.argv[1] if len(sys.argv) > 1 else '/home/dan/Music/Rumour_Cubes/The_Narrow_State/*.wav'

framelen = 2048
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)
panmul = ((stereobins-1)/2) / panlimit
bintofreq = fs/framelen

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

def analyse_file(audiopath):

	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:
			chunk = sf.read_frames(framelen, dtype=np.float32)
			if len(chunk) != framelen:
				print("Not read sufficient samples (%i) - returning." % (len(chunk)))
				break


			################################################################################
			# NOTE: This is weird - why do I need to flip the data around to get it correct?
			# Is this a bug in audiolab? I'm using scikits.audiolab-0.11.0 on ubuntu.
			chunk = chunk.T.reshape((framelen, 2))
			################################################################################

			# uncomment this to force the data into pure-mono, for testing
			#chunk[:,1] = chunk[:,0]

			framespectrum = np.array([np.fft.fft(window * chunk[:,channel].flatten()) for channel in [1,0]])

			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    # use pow rather than just 1? undecided.
		except RuntimeError:
			break
	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:
			print("Skipping a file due to ValueError (presumably wasn't stereo): %s" % matched)
			continue
		if bighisto == None:
			bighisto = histo
		else:
			bighisto += histo
	plot_histo(histo, "stereoscope.png")