week3/problem2/explore.py

   1 #!/usr/bin/env python
   2
   3 import sys
   4 import aifc
   5
   6 import pylab as pl
   7 import numpy as np
   8 from scipy import signal
   9 from numpy import pi
  10
  11
  12 fn = sys.argv[1]
  13
  14 sf = aifc.open(fn)
  15
  16 nchannels = sf.getnchannels()
  17 sampwidth = sf.getsampwidth()
  18 framerate = sf.getframerate()
  19 nframes = sf.getnframes()
  20
  21 if nchannels != 2:
  22         print "Expected stereo sound file"
  23         exit(1)
  24 if sampwidth != 2:
  25         print "Expected 16b samples"
  26         exit(1)
  27
  28 # I played with the sign and endianness until the histogram wasn't uniformly distributed
  29 data = np.fromstring(sf.readframes(nframes),dtype=np.int16).byteswap()/2.0**15
  30
  31 left = data[::nchannels]
  32 right = data[1::nchannels]
  33
  34 if left.size != right.size:
  35         print "left/right channels don't have the same number of samples, is this really stereo?"
  36         exit(1)
  37
  38 if left.size != nframes:
  39         print "Number of samples in left channel doesn't match frame count"
  40         exit(1)
  41
  42 #pl.hist(left,bins=30)
  43 #pl.show()
  44
  45 #pl.plot(left)
  46 #pl.plot(right)
  47 #pl.show()
  48
  49 #pl.psd(left,Fs=framerate,NFFT=nframes/16)
  50 #pl.show()
  51
  52 # Glancing at the first file, the 3 strongest tones are 261.5Hz, 391Hz, and 523.0Hz, or root fifth octave
  53 # The tuning appears to be pretty good.
  54
  55 # The hints suggest fast fourier transforming, but given that we know the tuning is good, and that modern computers are fast
  56 # We can afford to directly set up homodyne AM receivers on the western frequencies over a few octaves
  57 # The shortest notes in the score are 16th notes, and the quarter note is 152/minute
  58 # That gives a shortest note time of 0.0987 seconds
  59 # So we should set the bandwidth of the receivers at no less than 2*1/0.0987s, or in demodulated half bandwidth, 10Hz.
  60
  61 root_note = 55.0                # A pretty low A, frequency in Hz
  62 notes_per_octave = 12
  63 noctave = 6                             # 55,110,220,440,880,1660
  64
  65 if_bw = 20.0            # Hz
  66 decimation = 1
  67 nyquist = framerate/(2.0*decimation)
  68 nfilter = 8192
  69 if_filter = signal.firwin(nfilter,if_bw/nyquist)
  70
  71 t = np.arange(nframes,dtype=np.float64)/framerate       # Time in seconds
  72
  73 # Try C first
  74 def extract_note(sound,note):
  75         f = root_note*2**(note/12.0)
  76         print f
  77         k = 2*pi*(root_note*2**(note/12.0))
  78         c = np.cos(k*t)
  79         s = np.sin(k*t)
  80         I = signal.decimate(c*sound,decimation)
  81         I = signal.fftconvolve(I,if_filter)
  82
  83 #       pl.psd(I,NFFT=nframes/16,Fs=framerate)
  84 #       pl.show()
  85 #       exit()
  86
  87         Q = signal.decimate(s*sound,decimation)
  88         Q = signal.fftconvolve(Q,if_filter)
  89
  90         pl.plot(I)
  91         pl.plot(Q)
  92         pl.show()
  93         exit()
  94
  95         P = I*I + Q*Q
  96
  97         P = signal.fftconvolve(P,if_filter)
  98
  99         return P
 100
 101 #C = extract_note(left,12*2+3)
 102 A = extract_note(left,12*2+3+9)
 103 #C2 = extract_note(left,12*3+3)
 104
 105 #pl.plot(C)
 106 pl.plot(t,A)
 107 #pl.plot(C2)
 108 pl.show()