From e64488e0b91fbb5459b6a08b8fd210406141d214 Mon Sep 17 00:00:00 2001
From: Bryan Steinbach <bsteinbach@berkeley.edu>
Date: Thu, 23 Sep 2010 18:04:10 -0700
Subject: [PATCH] week3 problem 2 imported

---
 week3/problem2/explore.py | 108 ++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 108 insertions(+)
 create mode 100755 week3/problem2/explore.py

diff --git a/week3/problem2/explore.py b/week3/problem2/explore.py
new file mode 100755
index 0000000..a8c9d4b
--- /dev/null
+++ b/week3/problem2/explore.py
@@ -0,0 +1,108 @@
+#!/usr/bin/env python
+
+import sys
+import aifc
+
+import pylab as pl
+import numpy as np
+from scipy import signal
+from numpy import pi
+
+
+fn = sys.argv[1]
+
+sf = aifc.open(fn)
+
+nchannels = sf.getnchannels()
+sampwidth = sf.getsampwidth()
+framerate = sf.getframerate()
+nframes = sf.getnframes()
+
+if nchannels != 2:
+	print "Expected stereo sound file"
+	exit(1)
+if sampwidth != 2:
+	print "Expected 16b samples"
+	exit(1)
+
+# I played with the sign and endianness until the histogram wasn't uniformly distributed
+data = np.fromstring(sf.readframes(nframes),dtype=np.int16).byteswap()/2.0**15
+
+left = data[::nchannels]
+right = data[1::nchannels]
+
+if left.size != right.size:
+	print "left/right channels don't have the same number of samples, is this really stereo?"
+	exit(1)
+
+if left.size != nframes:
+	print "Number of samples in left channel doesn't match frame count"
+	exit(1)
+
+#pl.hist(left,bins=30)
+#pl.show()
+
+#pl.plot(left)
+#pl.plot(right)
+#pl.show()
+
+#pl.psd(left,Fs=framerate,NFFT=nframes/16)
+#pl.show()
+
+# Glancing at the first file, the 3 strongest tones are 261.5Hz, 391Hz, and 523.0Hz, or root fifth octave
+# The tuning appears to be pretty good.
+
+# The hints suggest fast fourier transforming, but given that we know the tuning is good, and that modern computers are fast
+# We can afford to directly set up homodyne AM receivers on the western frequencies over a few octaves
+# The shortest notes in the score are 16th notes, and the quarter note is 152/minute
+# That gives a shortest note time of 0.0987 seconds
+# So we should set the bandwidth of the receivers at no less than 2*1/0.0987s, or in demodulated half bandwidth, 10Hz.
+
+root_note = 55.0		# A pretty low A, frequency in Hz
+notes_per_octave = 12
+noctave = 6				# 55,110,220,440,880,1660
+
+if_bw = 20.0		# Hz
+decimation = 1
+nyquist = framerate/(2.0*decimation)
+nfilter = 8192
+if_filter = signal.firwin(nfilter,if_bw/nyquist)
+
+t = np.arange(nframes,dtype=np.float64)/framerate	# Time in seconds
+
+# Try C first
+def extract_note(sound,note):
+	f = root_note*2**(note/12.0)
+	print f
+	k = 2*pi*(root_note*2**(note/12.0))
+	c = np.cos(k*t)
+	s = np.sin(k*t)
+	I = signal.decimate(c*sound,decimation)
+	I = signal.fftconvolve(I,if_filter)
+
+#	pl.psd(I,NFFT=nframes/16,Fs=framerate)
+#	pl.show()
+#	exit()
+
+	Q = signal.decimate(s*sound,decimation)
+	Q = signal.fftconvolve(Q,if_filter)
+
+	pl.plot(I)
+	pl.plot(Q)
+	pl.show()
+	exit()
+
+	P = I*I + Q*Q
+
+	P = signal.fftconvolve(P,if_filter)
+
+	return P
+
+#C = extract_note(left,12*2+3)
+A = extract_note(left,12*2+3+9)
+#C2 = extract_note(left,12*3+3)
+
+#pl.plot(C)
+pl.plot(t,A)
+#pl.plot(C2)
+pl.show()
-- 
2.11.4.GIT