8 from scipy
import signal
16 nchannels
= sf
.getnchannels()
17 sampwidth
= sf
.getsampwidth()
18 framerate
= sf
.getframerate()
19 nframes
= sf
.getnframes()
22 print "Expected stereo sound file"
25 print "Expected 16b samples"
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
31 left
= data
[::nchannels
]
32 right
= data
[1::nchannels
]
34 if left
.size
!= right
.size
:
35 print "left/right channels don't have the same number of samples, is this really stereo?"
38 if left
.size
!= nframes
:
39 print "Number of samples in left channel doesn't match frame count"
42 #pl.hist(left,bins=30)
49 #pl.psd(left,Fs=framerate,NFFT=nframes/16)
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.
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.
61 root_note
= 55.0 # A pretty low A, frequency in Hz
63 noctave
= 6 # 55,110,220,440,880,1660
67 nyquist
= framerate
/(2.0*decimation
)
69 if_filter
= signal
.firwin(nfilter
,if_bw
/nyquist
)
71 t
= np
.arange(nframes
,dtype
=np
.float64
)/framerate
# Time in seconds
74 def extract_note(sound
,note
):
75 f
= root_note
*2**(note
/12.0)
77 k
= 2*pi
*(root_note
*2**(note
/12.0))
80 I
= signal
.decimate(c
*sound
,decimation
)
81 I
= signal
.fftconvolve(I
,if_filter
)
83 # pl.psd(I,NFFT=nframes/16,Fs=framerate)
87 Q
= signal
.decimate(s
*sound
,decimation
)
88 Q
= signal
.fftconvolve(Q
,if_filter
)
97 P
= signal
.fftconvolve(P
,if_filter
)
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)