#!/usr/pkg/bin/python
"""determine primary frequency of wave file"""
import sys, os, wave, pwd
sys.path.append(os.sep.join([pwd.getpwuid(os.geteuid())[5], 'lib', 'python']))
from com.jcomeau.midi import midi
from power import *
from struct import *
self = sys.argv.pop(0).split(os.sep)[-1].split('.')[0]
verbose, debugging = True, False
if len(sys.argv) < 1:
    sys.stderr.write("Usage: %s INFILE[...]\n" % self)
    sys.exit(0)

def verboseprint(*args):
 "helpful printout"
 global verbose
 if verbose:
  sys.stderr.write("%s\n" % repr(args))

def dbgprint(*args):
 "debugging printout"
 global debugging
 if debugging:
  sys.stderr.write("%s\n" % repr(args))

def soundchunk(wavefile, startposition = None, fraction = 1.0):
 "get a fractional-second chunk from sum of channels and return as list"
 channels = wavefile.getnchannels()
 samplewidth = wavefile.getsampwidth()
 framerate = wavefile.getframerate()
 dbgprint("getting sound chunk")
 format = {1: 'b', 2: 'h', 4: 'l'}
 if fraction > 1: fraction = 1
 chunksize = int(fraction * framerate)
 if startposition is not None:
  skip = int(startposition * framerate)
  wavefile.setpos(skip)
 chunk = wavefile.readframes(chunksize)
 chunk = unpack(format[samplewidth] * (len(chunk) / samplewidth), chunk)
 #dbgprint(chunk)
 total = array([0] * (len(chunk) / channels))
 dbgprint('total', total, 'channels', channels)
 for channel in range(0, channels):
  total = total + array(list(chunk)[channel::channels])
 dbgprint('total', total)
 return list(total)

def dcfilter(array):
 "eliminate DC component from signal"
 if len(array) == 0: return array
 dc_component = average(array)
 dbgprint("DC component is %d" % dc_component)
 dbgprint('before: ', array[0:10])
 for index in range(0, len(array)):
  array[index] = array[index] - dc_component
 dbgprint('after: ', array[0:10])
 return array
 
def average(array):
 "compute average of samples in array"
 return float(sum(array)) / len(array)

def averagepower(wavefile):
 length, total = 0, 0
 framerate = wavefile.getframerate()
 while True:
  try:
   chunk = soundchunk(wavefile)
   if len(chunk) == 0: raise Exception
  except:
   break
  chunkfft = powerspectrum(dcfilter(chunk), framerate).tolist()
  total = sum(chunkfft)
  length = length + len(chunkfft)
 wavefile.setpos(0)
 averagepower = float(total) / length
 verboseprint('average power: %.2f' % averagepower)
 return averagepower
 
def nearestnotes():
 nearest = dict()
 global notes
 twelfth = 1.0 / 12
 for note in notes.keys():
  # calculate half a semitone above and below
  minimum = notes[note] * pow(2, -twelfth / 2)
  maximum = notes[note] * pow(2, twelfth / 2)
  for frequency in range(int(minimum), int(maximum)):
   nearest[frequency] = note
 return nearest
 
notes = midi.NoteMapping('all', 'frequency')
midinotes = nearestnotes()
while len(sys.argv) > 0:
 filename = sys.argv.pop(0)
 wavefile = wave.open(filename, "rb")
 averagepower = averagepower(wavefile)
 channels = wavefile.getnchannels()
 samplewidth = wavefile.getsampwidth()
 framerate = wavefile.getframerate()
 nyquist = framerate / 2
 comptype = wavefile.getcomptype()
 length = .04 # seconds
 dbgprint("stats for %s (%s): channels=%d, samplewidth=%d, framerate=%d" % \
  (filename, comptype, channels, samplewidth, framerate))
 while True:
  try:
   chunk = soundchunk(wavefile, None, length)
   if len(chunk) == 0: raise Exception
  except:
   if debugging: raise
   else: break
  # weed DC out of signal
  chunkfft = powerspectrum(dcfilter(chunk), framerate).tolist()
  notelist = [] # make a list of strongest frequencies in each chunk
  peakfrequency = 0 # strength of strongest signal
  while True:
   # make sure to include nyquist frequency in check for max
   # ignore anything less than 30 Hz
   maximum = max(chunkfft[30:nyquist + 1])
   if maximum < 20 * averagepower: break
   #if maximum < 20 * peakfrequency: break
   if maximum > peakfrequency: peakfrequency = maximum
   frequency = chunkfft.index(maximum)
   try: note = frequency
   #try: note = midinotes[frequency]
   except: break
   if averagepower > 0: power = maximum / averagepower
   chunkfft[frequency] = 0
   #notesonly = map(lambda note: note[0], notelist)
   notesonly = notelist
   if note in notesonly:
    # add this frequency's strength to the primary frequency for that note
    notelist[notesonly.index(note)][1] += maximum
    if notelist[notesonly.index(note)][1] > peakfrequency:
     peakfrequency = notelist[notesonly.index(note)][1]
   else:
    try: notelist.append([note, power])
    except: raise
  print repr(notelist)
 wavefile.close()