Repository: NickleDave/songdkl (language: python)

Issue:
make sure all functions have clear docstrings



Discussion / hints:

Fix patch:
diff --git a/.gitignore b/.gitignore
index 2142d19..3afd33b 100644
--- a/.gitignore
+++ b/.gitignore
@@ -7,7 +7,10 @@
 # jupyter
 .ipynb_checkpoints/
 
-# python
+# Mac
+.DS_Store
+
+#python
 *__pycache__*
 *egg-info/
 dist/
diff --git a/pyproject.toml b/pyproject.toml
index d9735b2..8ce8355 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -36,6 +36,7 @@ dev = [
     "twine",
 ]
 doc = [
+    "jupyterlab >=3.0.3",
     "Sphinx >= 3.4.1",
 ]
 
diff --git a/src/pcb-scripts/Song_D_KL_calc.py b/src/pcb-scripts/Song_D_KL_calc.py
index ba7f018..eae2dd8 100644
--- a/src/pcb-scripts/Song_D_KL_calc.py
+++ b/src/pcb-scripts/Song_D_KL_calc.py
@@ -33,202 +33,200 @@
 
 IMPORTANT! Throughout the paper we calculated PSDs for the raw wave forms of syllables.  This is the default setting for this script. If your song is contaminated with low frequency noise this noise may be incorporated into the model for song potentially causing over estimates of song D_KL if there is low frequency noise in the tutor song but not the tutee song.  This could occur if the birds were recorded under different conditions or in different sound recording boxes.  If you uncomment line 99 below, the script will calculate the song D_KL using filtered syllable data. We only advise this if you have low frequency noise that is differential between the tutor and tutee song and you don't want that noise incorporated into the song D_KL calculations."""
 
-#functions
+
+# functions
 
 def impwav(a):
- """Imports a wave file as an array where a[1] 
- is the sampling frequency and a[0] is the data"""
- out=[]
- wav = sc.io.wavfile.read(a)
- out=[wav[1],wav[0]]
- return out
+    """Imports a wave file as an array where a[1]
+    is the sampling frequency and a[0] is the data"""
+    out = []
+    wav = sc.io.wavfile.read(a)
+    out = [wav[1], wav[0]]
+    return out
+
 
 def norm(a):
- """normalizes a string by it's average and sd"""
- a=(np.array(a)-np.average(a))/np.std(a)
- return a
+    """normalizes a string by it's average and sd"""
+    a = (np.array(a) - np.average(a)) / np.std(a)
+    return a
+
 
 def filtersong(a):
- """highpass iir filter for song."""
- out=[]
- b=sc.signal.iirdesign(wp=0.04, ws=0.02, gpass=1, gstop=60,ftype='ellip')
- out.append(sc.signal.filtfilt(b[0],b[1],a[0]))
- out.append(a[1])
- return (out)
-
-def threshold(a,thresh=None):
- """Returns a thresholded array of the same length as input
- with everything below a specific threshold set to 0.
- By default threshold is sigma."""
- if thresh==None: thresh = sc.std(a)
- out=np.where(abs(a)>thresh,a,np.zeros(a.shape))
- return out
+    """highpass iir filter for song."""
+    out = []
+    b = sc.signal.iirdesign(wp=0.04, ws=0.02, gpass=1, gstop=60, ftype='ellip')
+    out.append(sc.signal.filtfilt(b[0], b[1], a[0]))
+    out.append(a[1])
+    return (out)
 
-def findobject(file):
- """finds objects.  Expects a smoothed rectified amplitude envelope"""
- value=(otsu(np.array(file,dtype=np.uint32)))/2 #calculate a threshold
- #value=(np.average(file))/2 #heuristically, this also usually works  for establishing threshold
- thresh=threshold(file,value) #threshold the envelope data
- thresh=threshold(sc.ndimage.convolve(thresh,np.ones(512)),0.5) #pad the threshold
- label=(sc.ndimage.label(thresh)[0]) #label objects in the threshold
- objs=sc.ndimage.find_objects(label) #recover object positions
- return(objs)
-
-def smoothrect(a,window=None,freq=None):
- """smooths and rectifies a song.  Expects (data,samprate)"""
- if freq== None: freq=32000 #baseline values if none are provided
- if window == None: window=2 #baseline if none are provided
- le=int(round(freq*window/1000)) #calculate boxcar kernel length
- h=np.ones(le)/le #make boxcar 
- smooth= np.convolve(h,abs(a)) # convovlve boxcar with signal
- offset = int(round((len(smooth)-len(a))/2)) #calculate offset imposed by convolution
- smooth=smooth[(1+offset):(len(a)+offset)] #correct for offset
- return smooth
 
-def getsyls(a):
- """takes a file read in with impwav and returns a list of sylables"""
- fa=filtersong(a) #filter song input
- frq=a[1] #get sampling frequency from data (a)
- a=a[0] #get data from data (a)
- frqs=frq/1000 # calcualte length of a ms in samples
- objs=findobject(smoothrect(fa[0],10,frq)) #get syllable positions
- sylables=[x for x in [a[y] for y in objs] if int(len(x))>(10*frqs)] #get syllable data if of sufficient duration
- '''uncomment the next line to recover syllables that have been high pass filtered as opposed to raw data.
- Using data filtered prior to PSD calculation helps if you data are contaminated
- with low frequency noise''' 
- #sylables=[x for x in [fa[0][y] for y in objs] if int(len(x))>(10*frqs)] #get syllable data if of sufficient duration.
- objs=[y for y in objs if int(len(a[y]))>10*frqs] #get objects of sufficient duraiton
- return sylables,objs,frq
+def threshold(a, thresh=None):
+    """Returns a thresholded array of the same length as input
+    with everything below a specific threshold set to 0.
+    By default threshold is sigma."""
+    if thresh == None: thresh = sc.std(a)
+    out = np.where(abs(a) > thresh, a, np.zeros(a.shape))
+    return out
 
 
-#main program
+def findobject(file):
+    """finds objects.  Expects a smoothed rectified amplitude envelope"""
+    value = (otsu(np.array(file, dtype=np.uint32))) / 2  # calculate a threshold
+    # value=(np.average(file))/2 #heuristically, this also usually works  for establishing threshold
+    thresh = threshold(file, value)  # threshold the envelope data
+    thresh = threshold(sc.ndimage.convolve(thresh, np.ones(512)), 0.5)  # pad the threshold
+    label = (sc.ndimage.label(thresh)[0])  # label objects in the threshold
+    objs = sc.ndimage.find_objects(label)  # recover object positions
+    return (objs)
+
+
+def smoothrect(a, window=None, freq=None):
+    """smooths and rectifies a song.  Expects (data,samprate)"""
+    if freq == None: freq = 32000  # baseline values if none are provided
+    if window == None: window = 2  # baseline if none are provided
+    le = int(round(freq * window / 1000))  # calculate boxcar kernel length
+    h = np.ones(le) / le  # make boxcar
+    smooth = np.convolve(h, abs(a))  # convovlve boxcar with signal
+    offset = int(round((len(smooth) - len(a)) / 2))  # calculate offset imposed by convolution
+    smooth = smooth[(1 + offset):(len(a) + offset)]  # correct for offset
+    return smooth
 
-path1=sys.argv[1]
-path2=sys.argv[2]
 
-fils1=[x for x in os.listdir(path1) if x[-4:]== '.wav']
-fils2=[x for x in os.listdir(path2) if x[-4:]== '.wav']
+def getsyls(a):
+    """takes a file read in with impwav and returns a list of sylables"""
+    fa = filtersong(a)  # filter song input
+    frq = a[1]  # get sampling frequency from data (a)
+    a = a[0]  # get data from data (a)
+    frqs = frq / 1000  # calcualte length of a ms in samples
+    objs = findobject(smoothrect(fa[0], 10, frq))  # get syllable positions
+    sylables = [x for x in [a[y] for y in objs] if
+                int(len(x)) > (10 * frqs)]  # get syllable data if of sufficient duration
+    '''uncomment the next line to recover syllables that have been high pass filtered as opposed to raw data.
+    Using data filtered prior to PSD calculation helps if you data are contaminated
+    with low frequency noise'''
+    # sylables=[x for x in [fa[0][y] for y in objs] if int(len(x))>(10*frqs)] #get syllable data if of sufficient duration.
+    objs = [y for y in objs if int(len(a[y])) > 10 * frqs]  # get objects of sufficient duraiton
+    return sylables, objs, frq
 
-fils1=fils1[:120]
-fils2=fils2[:120]
 
-filename1=path1.split('/')[-2]
-filename2=path2.split('/')[-2]
+# main program
 
-datano=10000 #variable to constrain the amount of data used
+path1 = sys.argv[1]
+path2 = sys.argv[2]
 
-k=int(sys.argv[3]) #number of syllable classes in song 1
-k2=int(sys.argv[4]) #number of syllable classes in song 2
+fils1 = [x for x in os.listdir(path1) if x[-4:] == '.wav']
+fils2 = [x for x in os.listdir(path2) if x[-4:] == '.wav']
 
+fils1 = fils1[:120]
+fils2 = fils2[:120]
 
+filename1 = path1.split('/')[-2]
+filename2 = path2.split('/')[-2]
 
+datano = 10000  # variable to constrain the amount of data used
 
+k = int(sys.argv[3])  # number of syllable classes in song 1
+k2 = int(sys.argv[4])  # number of syllable classes in song 2
 
-#get syllables from all of the songs in folder 1
-syls1=[]
-objss1=[]
+# get syllables from all of the songs in folder 1
+syls1 = []
+objss1 = []
 for file in fils1:
- song=impwav(path1+file)
- if len(song[0])<1: break
- syls,objs,frq=(getsyls(song))
- objss1.append(objs)
- syls1.append([frq]+syls)
-
-
- 
-#get syllables from all of the songs in folder 2
-syls2=[]
-objss2=[]
+    song = impwav(path1 + file)
+    if len(song[0]) < 1: break
+    syls, objs, frq = (getsyls(song))
+    objss1.append(objs)
+    syls1.append([frq] + syls)
+
+# get syllables from all of the songs in folder 2
+syls2 = []
+objss2 = []
 for file in fils2:
- song=impwav(path2+file)
- if len(song[0])<1: break
- syls,objs,frq=(getsyls(song))
- objss2.append(objs)
- syls2.append([frq]+syls)
-
-
-#convert syllables into PSDs
-segedpsds1=[]
+    song = impwav(path2 + file)
+    if len(song[0]) < 1: break
+    syls, objs, frq = (getsyls(song))
+    objss2.append(objs)
+    syls2.append([frq] + syls)
+
+# convert syllables into PSDs
+segedpsds1 = []
 for x in syls1[:datano]:
- fs=x[0]
- nfft=int(round(2**14/32000.0*fs))
- segstart=int(round(600/(fs/float(nfft))))
- segend=int(round(16000/(fs/float(nfft))))
- psds=[psd(norm(y),NFFT=nfft,Fs=fs) for y in x[1:]]
- spsds=[norm(n[0][segstart:segend]) for n in psds]
- for n in spsds: segedpsds1.append(n)
- if len(segedpsds1)>datano: break
-
-#convert syllables to PSDs
-segedpsds2=[]
+    fs = x[0]
+    nfft = int(round(2 ** 14 / 32000.0 * fs))
+    segstart = int(round(600 / (fs / float(nfft))))
+    segend = int(round(16000 / (fs / float(nfft))))
+    psds = [psd(norm(y), NFFT=nfft, Fs=fs) for y in x[1:]]
+    spsds = [norm(n[0][segstart:segend]) for n in psds]
+    for n in spsds: segedpsds1.append(n)
+    if len(segedpsds1) > datano: break
+
+# convert syllables to PSDs
+segedpsds2 = []
 for x in syls2[:datano]:
- fs=x[0]
- nfft=int(round(2**14/32000.0*fs))
- segstart=int(round(600/(fs/float(nfft))))
- segend=int(round(16000/(fs/float(nfft))))
- psds=[psd(norm(y),NFFT=nfft,Fs=fs) for y in x[1:]]
- spsds=[norm(n[0][segstart:segend]) for n in psds]
- for n in spsds: segedpsds2.append(n)
- if len(segedpsds2)>datano: break
+    fs = x[0]
+    nfft = int(round(2 ** 14 / 32000.0 * fs))
+    segstart = int(round(600 / (fs / float(nfft))))
+    segend = int(round(16000 / (fs / float(nfft))))
+    psds = [psd(norm(y), NFFT=nfft, Fs=fs) for y in x[1:]]
+    spsds = [norm(n[0][segstart:segend]) for n in psds]
+    for n in spsds: segedpsds2.append(n)
+    if len(segedpsds2) > datano: break
 
-#establish basis set:
-basis_set=segedpsds1[:50] #select the first 50 syllables of the reference song as the basis set
-#basis_set=[segedpsds1[rnd.randint(0,len(segedpsds1))] for x in range(50)] #select a random set of 50 syllablse as the basis set
+# establish basis set:
+basis_set = segedpsds1[:50]  # select the first 50 syllables of the reference song as the basis set
+# basis_set=[segedpsds1[rnd.randint(0,len(segedpsds1))] for x in range(50)] #select a random set of 50 syllablse as the basis set
 
 
-lone=len(segedpsds1)
-ltwo=len(segedpsds2)
-lone_half=int(lone/2)
-ltwo_half=int(ltwo/2)
+lone = len(segedpsds1)
+ltwo = len(segedpsds2)
+lone_half = int(lone / 2)
+ltwo_half = int(ltwo / 2)
 
-#calculate distance matrices:
-d1=sc.spatial.distance.cdist(segedpsds1[:lone_half],basis_set,'sqeuclidean')
+# calculate distance matrices:
+d1 = sc.spatial.distance.cdist(segedpsds1[:lone_half], basis_set, 'sqeuclidean')
 
-d1_2=sc.spatial.distance.cdist(segedpsds1[lone_half:lone],basis_set,'sqeuclidean')
+d1_2 = sc.spatial.distance.cdist(segedpsds1[lone_half:lone], basis_set, 'sqeuclidean')
 
-d2=sc.spatial.distance.cdist(segedpsds2[:ltwo_half],basis_set,'sqeuclidean')
+d2 = sc.spatial.distance.cdist(segedpsds2[:ltwo_half], basis_set, 'sqeuclidean')
 
-d2_2=sc.spatial.distance.cdist(segedpsds2[ltwo_half:ltwo],basis_set,'sqeuclidean')
+d2_2 = sc.spatial.distance.cdist(segedpsds2[ltwo_half:ltwo], basis_set, 'sqeuclidean')
 
-mx=np.max([np.max(d1),np.max(d2),np.max(d1_2),np.max(d2_2)])
+mx = np.max([np.max(d1), np.max(d2), np.max(d1_2), np.max(d2_2)])
 
-#convert to similarity matrices:
-s1=1-(d1/mx)
-s1_2=1-(d1_2/mx)
-s2=1-(d2/mx)
-s2_2=1-(d2_2/mx)
+# convert to similarity matrices:
+s1 = 1 - (d1 / mx)
+s1_2 = 1 - (d1_2 / mx)
+s2 = 1 - (d2 / mx)
+s2_2 = 1 - (d2_2 / mx)
 
-#estimate GMMs:
-mod1=mixture.GMM(n_components=k,n_iter=100000,n_init=5,covariance_type='full')
+# estimate GMMs:
+mod1 = mixture.GMM(n_components=k, n_iter=100000, n_init=5, covariance_type='full')
 mod1.fit(s1)
 
-mod2=mixture.GMM(n_components=k2,n_iter=100000,n_init=5,covariance_type='full')
+mod2 = mixture.GMM(n_components=k2, n_iter=100000, n_init=5, covariance_type='full')
 mod2.fit(s2)
 
-len2=len(s2)
-len1=len(d1)
-
-#calculate likelihoods for held out data:
-score1_1=mod1.score(s1_2)
-score2_1=mod2.score(s1_2)
-
-score1_2=mod1.score(s2_2)
-score2_2=mod2.score(s2_2)
-
-len2=float(len(basis_set))
-len1=float(len(basis_set))
-
-#calculate song divergence (DKL estimate):
-score1= np.log2(np.e)*((np.mean(score1_1))-(np.mean(score2_1)))
-score2= np.log2(np.e)*((np.mean(score2_2))-(np.mean(score1_2)))
-
-score1=score1/len1
-score2=score2/len2
+len2 = len(s2)
+len1 = len(d1)
 
+# calculate likelihoods for held out data:
+score1_1 = mod1.score(s1_2)
+score2_1 = mod2.score(s1_2)
 
-#output
-print filename1+'\t'+filename2+'\t'+str(k)+'\t'+str(k2)+'\t'+str(len2)+'\t'+str(score1)+'\t'+str(score2)+'\t'+str(len(segedpsds1))+'\t'+str(len(segedpsds2))
+score1_2 = mod1.score(s2_2)
+score2_2 = mod2.score(s2_2)
 
+len2 = float(len(basis_set))
+len1 = float(len(basis_set))
 
+# calculate song divergence (DKL estimate):
+score1 = np.log2(np.e) * ((np.mean(score1_1)) - (np.mean(score2_1)))
+score2 = np.log2(np.e) * ((np.mean(score2_2)) - (np.mean(score1_2)))
 
+score1 = score1 / len1
+score2 = score2 / len2
 
+# output
+print
+filename1 + '\t' + filename2 + '\t' + str(k) + '\t' + str(k2) + '\t' + str(len2) + '\t' + str(score1) + '\t' + str(
+    score2) + '\t' + str(len(segedpsds1)) + '\t' + str(len(segedpsds2))
diff --git a/src/pcb-scripts/calc_sylno.py b/src/pcb-scripts/calc_sylno.py
index d60859f..cf7fac0 100644
--- a/src/pcb-scripts/calc_sylno.py
+++ b/src/pcb-scripts/calc_sylno.py
@@ -30,156 +30,168 @@
 foldername_bird	number_of_syllables
 e.g. y34br6 9"""
 
-#functions
+
+# functions
 def impwav(a):
- """Imports a wave file as an array where a[1] 
- is the sampling frequency and a[0] is the data"""
- out=[]
- wav = sc.io.wavfile.read(a)
- out=[wav[1],wav[0]]
- return out
+    """Imports a wave file as an array where a[1]
+    is the sampling frequency and a[0] is the data"""
+    out = []
+    wav = sc.io.wavfile.read(a)
+    out = [wav[1], wav[0]]
+    return out
+
 
 def norm(a):
- """normalizes a string by it's average and sd"""
- a=(np.array(a)-np.average(a))/np.std(a)
- return a
+    """normalizes a string by it's average and sd"""
+    a = (np.array(a) - np.average(a)) / np.std(a)
+    return a
+
 
 def filtersong(a):
- """highpass iir filter for song."""
- out=[]
- b=sc.signal.iirdesign(wp=0.04, ws=0.02, gpass=1, gstop=60,ftype='ellip')
- out.append(sc.signal.filtfilt(b[0],b[1],a[0]))
- out.append(a[1])
- return (out)
-
-def threshold(a,thresh=None):
- """Returns a thresholded array of the same length as input
- with everything below a specific threshold set to 0.
- By default threshold is sigma."""
- if thresh==None: thresh = sc.std(a)
- out=np.where(abs(a)>thresh,a,np.zeros(a.shape))
- return out
+    """highpass iir filter for song."""
+    out = []
+    b = sc.signal.iirdesign(wp=0.04, ws=0.02, gpass=1, gstop=60, ftype='ellip')
+    out.append(sc.signal.filtfilt(b[0], b[1], a[0]))
+    out.append(a[1])
+    return (out)
+
+
+def threshold(a, thresh=None):
+    """Returns a thresholded array of the same length as input
+    with everything below a specific threshold set to 0.
+    By default threshold is sigma."""
+    if thresh == None: thresh = sc.std(a)
+    out = np.where(abs(a) > thresh, a, np.zeros(a.shape))
+    return out
+
 
 def findobject(file):
- """finds objects.  Expects a smoothed rectified amplitude envelope"""
- value=(otsu(np.array(file,dtype=np.uint32)))/2 #calculate a threshold
- #value=(np.average(file))/2 #heuristically, this also usually works  for establishing a threshold
- thresh=threshold(file,value) #threshold the envelope data
- thresh=threshold(sc.ndimage.convolve(thresh,np.ones(512)),0.5) #pad the threshold
- label=(sc.ndimage.label(thresh)[0]) #label objects in the threshold
- objs=sc.ndimage.find_objects(label) #recover object positions
- return(objs)
-
-def smoothrect(a,window=None,freq=None):
- """smooths and rectifies a song.  Expects (data,samprate)"""
- if freq== None: freq=32000 #baseline values if none are provided
- if window == None: window=2 #baseline if none are provided
- le=int(round(freq*window/1000)) #calculate boxcar kernel length
- h=np.ones(le)/le #make boxcar 
- smooth= np.convolve(h,abs(a)) # convovlve boxcar with signal
- offset = int(round((len(smooth)-len(a))/2)) #calculate offset imposed by convolution
- smooth=smooth[(1+offset):(len(a)+offset)] #correct for offset
- return smooth
+    """finds objects.  Expects a smoothed rectified amplitude envelope"""
+    value = (otsu(np.array(file, dtype=np.uint32))) / 2  # calculate a threshold
+    # value=(np.average(file))/2 #heuristically, this also usually works  for establishing a threshold
+    thresh = threshold(file, value)  # threshold the envelope data
+    thresh = threshold(sc.ndimage.convolve(thresh, np.ones(512)), 0.5)  # pad the threshold
+    label = (sc.ndimage.label(thresh)[0])  # label objects in the threshold
+    objs = sc.ndimage.find_objects(label)  # recover object positions
+    return (objs)
+
+
+def smoothrect(a, window=None, freq=None):
+    """smooths and rectifies a song.  Expects (data,samprate)"""
+    if freq == None: freq = 32000  # baseline values if none are provided
+    if window == None: window = 2  # baseline if none are provided
+    le = int(round(freq * window / 1000))  # calculate boxcar kernel length
+    h = np.ones(le) / le  # make boxcar
+    smooth = np.convolve(h, abs(a))  # convovlve boxcar with signal
+    offset = int(round((len(smooth) - len(a)) / 2))  # calculate offset imposed by convolution
+    smooth = smooth[(1 + offset):(len(a) + offset)]  # correct for offset
+    return smooth
+
 
 def getsyls(a):
- """takes a file read in with impwav and returns a list of sylables"""
- fa=filtersong(a) #filter song input
- frq=a[1] #get sampling frequency from data (a)
- a=a[0] #get data from data (a)
- frqs=frq/1000 # calcualte length of a ms in samples
- objs=findobject(smoothrect(fa[0],10,frq)) #get syllable positions
- sylables=[x for x in [a[y] for y in objs] if int(len(x))>(10*frqs)] #get syllable dat$
- objs=[y for y in objs if int(len(a[y]))>10*frqs] #get objects of sufficient duraiton
- return sylables,objs,frq
+    """takes a file read in with impwav and returns a list of sylables"""
+    fa = filtersong(a)  # filter song input
+    frq = a[1]  # get sampling frequency from data (a)
+    a = a[0]  # get data from data (a)
+    frqs = frq / 1000  # calcualte length of a ms in samples
+    objs = findobject(smoothrect(fa[0], 10, frq))  # get syllable positions
+    sylables = [x for x in [a[y] for y in objs] if int(len(x)) > (10 * frqs)]  # get syllable dat$
+    objs = [y for y in objs if int(len(a[y])) > 10 * frqs]  # get objects of sufficient duraiton
+    return sylables, objs, frq
+
 
 def chunks(l, n):
- """ Yield successive n-sized chunks from l.
- """
- out=[]
- lens=(len(l)/n)
- for i in range(0, lens):
-  out.append(l[i*n:i*n+n])
- out.append(l[lens*n:])
- return out[:-1]
+    """ Yield successive n-sized chunks from l.
+    """
+    out = []
+    lens = (len(l) / n)
+    for i in range(0, lens):
+        out.append(l[i * n:i * n + n])
+    out.append(l[lens * n:])
+    return out[:-1]
+
 
 def EMofgmmcluster(array_of_syls):
- """takes an array of segmented sylables, fits a series of
- gaussian mixture models with an increasing number of mixtures, 
- and identifies the best number of mixtures to describe the data
- by BIC. If data is limiting, we don't recommend using cross validation 
- for BIC.  If you have sufficient data 3 fold cross validated BIC values
- compare well with un-cross validated values and human assesment of syllable
- number.  If data are limiting, cross validated values underestimate
- the number of syllables relative to human assessment or non-cross
- validated BIC values.
- """
- nfft=2**14
- fs=32000
- segstart=int(600/(fs/nfft))
- segend=int(16000/(fs/nfft))
-
- welchpsds=[psd(x,NFFT=nfft,Fs=fs) for x in array_of_syls]
- segedpsds1=[norm(x[0][segstart:segend]) for x in welchpsds]
-
- models=range(2,22)
-
- basis_set=segedpsds1[:50] #select the first 50 syllables of the reference song as the basis set
- #basis_set=[segedpsds1[rnd.randint(0,len(segedpsds1))] for x in range(50)] #select a random set of 50 syllablse as the basis set
- d=sc.spatial.distance.cdist(segedpsds1,basis_set,'sqeuclidean')
- s=1-d/np.max(d)*1000
- bic=[]
- #xv=3
- for x in models:
-  bics=[]
-  #this commented section implements cross validation for the BIC values
-  '''ss=chunks(s,len(s)/xv)
-  for y in range(len(ss)):
-   testset=ss[y]
-   trainset=[]
-   [trainset.extend(ss[n]) for n in [z for z in range(len(ss)) if z!=y]]
-   gmm = mixture.GMM(n_components=x, n_iter=100000,n_init=5, covariance_type='full')
-   gmm.fit(np.array(s))
-   bics.append(gmm.bic(np.array(s)))
-  bic.append(np.mean(bics))'''
-  gmm = mixture.GMM(n_components=x, n_iter=100000,n_init=5, covariance_type='full')
-  gmm.fit(np.array(s))
-  bic.append(gmm.bic(np.array(s)))
- return segedpsds1,bic.index(min(bic))+2
-  
-#main program
-path1=sys.argv[1]
-
-fils1=[x for x in os.listdir(path1) if x[-4:]== '.wav']
-
-fils1=fils1[:120]
-
-filename1=path1.split('/')[-2]
-
-datano=15000 #maximum number of syllables to be analyzed
-
-#get syllables from all of the songs in folder 
-syls1=[]
-objss1=[]
+    """takes an array of segmented sylables, fits a series of
+    gaussian mixture models with an increasing number of mixtures,
+    and identifies the best number of mixtures to describe the data
+    by BIC. If data is limiting, we don't recommend using cross validation
+    for BIC.  If you have sufficient data 3 fold cross validated BIC values
+    compare well with un-cross validated values and human assesment of syllable
+    number.  If data are limiting, cross validated values underestimate
+    the number of syllables relative to human assessment or non-cross
+    validated BIC values.
+    """
+    nfft = 2 ** 14
+    fs = 32000
+    segstart = int(600 / (fs / nfft))
+    segend = int(16000 / (fs / nfft))
+
+    welchpsds = [psd(x, NFFT=nfft, Fs=fs) for x in array_of_syls]
+    segedpsds1 = [norm(x[0][segstart:segend]) for x in welchpsds]
+
+    models = range(2, 22)
+
+    basis_set = segedpsds1[:50]  # select the first 50 syllables of the reference song as the basis set
+    # basis_set=[segedpsds1[rnd.randint(0,len(segedpsds1))] for x in range(50)] #select a random set of 50 syllablse as the basis set
+    d = sc.spatial.distance.cdist(segedpsds1, basis_set, 'sqeuclidean')
+    s = 1 - d / np.max(d) * 1000
+    bic = []
+    # xv=3
+    for x in models:
+        bics = []
+        # this commented section implements cross validation for the BIC values
+        '''ss=chunks(s,len(s)/xv)
+        for y in range(len(ss)):
+         testset=ss[y]
+         trainset=[]
+         [trainset.extend(ss[n]) for n in [z for z in range(len(ss)) if z!=y]]
+         gmm = mixture.GMM(n_components=x, n_iter=100000,n_init=5, covariance_type='full')
+         gmm.fit(np.array(s))
+         bics.append(gmm.bic(np.array(s)))
+        bic.append(np.mean(bics))'''
+        gmm = mixture.GMM(n_components=x, n_iter=100000, n_init=5, covariance_type='full')
+        gmm.fit(np.array(s))
+        bic.append(gmm.bic(np.array(s)))
+    return segedpsds1, bic.index(min(bic)) + 2
+
+
+# main program
+path1 = sys.argv[1]
+
+fils1 = [x for x in os.listdir(path1) if x[-4:] == '.wav']
+
+fils1 = fils1[:120]
+
+filename1 = path1.split('/')[-2]
+
+datano = 15000  # maximum number of syllables to be analyzed
+
+# get syllables from all of the songs in folder
+syls1 = []
+objss1 = []
 for file in fils1:
- song=impwav(path1+file)
- if len(song[0])<1: break
- syls,objs,frq=(getsyls(song))
- objss1.append(objs)
- syls1.append([frq]+syls)
-
-#convert syllables so PSDs
-segedpsds1=[]
+    song = impwav(path1 + file)
+    if len(song[0]) < 1: break
+    syls, objs, frq = (getsyls(song))
+    objss1.append(objs)
+    syls1.append([frq] + syls)
+
+# convert syllables so PSDs
+segedpsds1 = []
 for x in syls1[:datano]:
- fs=x[0]
- nfft=int(round(2**14/32000.0*fs))
- segstart=int(round(600/(fs/float(nfft))))
- segend=int(round(16000/(fs/float(nfft))))
- psds=[psd(norm(y),NFFT=nfft,Fs=fs) for y in x[1:]]
- spsds=[norm(n[0][segstart:segend]) for n in psds]
- for n in spsds: segedpsds1.append(n)
- if len(segedpsds1)>datano: break
-
-
-segedpsds,sylno_bic=EMofgmmcluster(segedpsds1)
-
-print filename1+'\t'+str(sylno_bic)
+    fs = x[0]
+    nfft = int(round(2 ** 14 / 32000.0 * fs))
+    segstart = int(round(600 / (fs / float(nfft))))
+    segend = int(round(16000 / (fs / float(nfft))))
+    # not clear to me why norm is applied twice: once here
+    psds = [psd(norm(y), NFFT=nfft, Fs=fs) for y in x[1:]]
+    # and then applied here again. Might have no effect at this point?
+    spsds = [norm(n[0][segstart:segend]) for n in psds]
+    for n in spsds: segedpsds1.append(n)
+    if len(segedpsds1) > datano: break
+
+# we pass in segedpsds1 but then these get run through the code again that makes PSDs
+segedpsds, sylno_bic = EMofgmmcluster(segedpsds1)
+
+print filename1 + '\t' + str(sylno_bic)
diff --git a/src/songdkl/__main__.py b/src/songdkl/__main__.py
index f21f22f..8d46802 100644
--- a/src/songdkl/__main__.py
+++ b/src/songdkl/__main__.py
@@ -9,13 +9,26 @@ def main():
     args = parser.parse_args()
 
     if args.command == 'calculate':
-        calculate(path1=args.path1,
-                  path2=args.path2,
-                  k=args.n_syl1,
-                  k2=args.n_syl2)
+        score1, score2, n_psds_ref, n_psds_compare = calculate(ref_dir_path=args.ref_dir_path,
+                                                               compare_dir_path=args.compare_dir_path,
+                                                               k_ref=args.k_ref,
+                                                               k_compare=args.k_compare,
+                                                               max_wavs=args.max_wavs,
+                                                               max_num_psds=args.max_num_psds,
+                                                               n_basis=args.n_basis,
+                                                               basis=args.basis)
+        print(
+            f'{args.ref_dir_path}\t{args.compare_dir_path}\t'
+            f'{args.k_ref}\t{args.k_compare}\t'
+            f'{args.n_basis}\t{score1}\t{score2}\t'
+            f'{n_psds_ref}\t{n_psds_compare}'
+        )
 
     elif args.command == 'numsyls':
-        numsyls(path1=args.path)
+        n_syls = numsyls(dir_path=args.dir_path)
+        print(
+            f'{args.dir_path}\t{n_syls}'
+        )
 
     elif args.command is None:
         parser.print_help()
diff --git a/src/songdkl/argparser/argparser.py b/src/songdkl/argparser/argparser.py
index c8e345a..462773c 100644
--- a/src/songdkl/argparser/argparser.py
+++ b/src/songdkl/argparser/argparser.py
@@ -5,7 +5,8 @@
 
 def get():
     """creates argparser, used by __main__.main function"""
-    parser = argparse.ArgumentParser(description='a similarity metric for birdsong',
+    parser = argparse.ArgumentParser(prog='songdkl',
+                                     description='a similarity metric for birdsong',
                                      formatter_class=argparse.RawTextHelpFormatter,
                                      epilog=PARSER_EPILOG)
 
@@ -17,18 +18,32 @@ def get():
     calculate_subparser = subparser.add_parser('calculate',
                                                help='calculate the song divergence between two birds',
                                                epilog=CALCULATE_EPILOG)
-    calculate_subparser.add_argument('path1', type=str,
-                                     help='path to directory containing songs from bird 1')
-    calculate_subparser.add_argument('path2', type=str,
-                                     help='path to directory containing songs from bird 1')
-    calculate_subparser.add_argument('n_syl1', type=int,
-                                     help='number of syllable_classes for bird 1')
-    calculate_subparser.add_argument('n_syl2', type=int,
-                                     help='number of syllable_classes for bird 2')
+    calculate_subparser.add_argument('ref_dir_path', type=str,
+                                     help=('Path to directory with .wav files of songs from bird '
+                                           'that should be used as reference.'))
+    calculate_subparser.add_argument('compare_dir_path', type=str,
+                                     help=('Path to directory with .wav files of songs from bird '
+                                           'that should be compared with reference.'))
+    calculate_subparser.add_argument('k_ref', type=int,
+                                     help=('Number of syllable classes in song of bird used as reference.'
+                                           'Also the number of components $k_{ref}$ used for Gaussian Mixture '
+                                           'Model fit to the reference distances.'))
+    calculate_subparser.add_argument('k_compare', type=int,
+                                     help=('Number of syllable classes in song of bird compared with reference.'
+                                           'Also the number of components $k_{compare}$ used for Gaussian Mixture '
+                                           'Model fit to the comparison distances.'))
+    calculate_subparser.add_argument('--max_wavs', type=int, default=120,
+                                     help='Maximum number of .wav files to use. Default  is 120.')
+    calculate_subparser.add_argument('--max_num_psds', type=str, default=10000,
+                                     help='Maximum number of power spectral densities (PSDs) to use. Default is 10000.')
+    calculate_subparser.add_argument('--n_basis', type=int, default=50,
+                                     help='Number of PSDs to use for the basis set. Default is 50.')
+    calculate_subparser.add_argument('--basis', type=str, default='first', choices={'first', 'random'},
+                                     help="How to select PSDs for basis set. Either 'first' (default) or 'random'")
 
     numsyls_subparser = subparser.add_parser('numsyls',
                                              help="estimate the number of syllable types in a song.",
                                              epilog=NUMSYLS_EPILOG)
-    numsyls_subparser.add_argument('path', type=str, help='path to directory containing songs from bird')
+    numsyls_subparser.add_argument('dir_path', type=str, help='path to directory containing songs from bird')
 
     return parser
diff --git a/src/songdkl/audio.py b/src/songdkl/audio.py
index dcf77f4..7cdfa6c 100644
--- a/src/songdkl/audio.py
+++ b/src/songdkl/audio.py
@@ -1,72 +1,144 @@
+from __future__ import annotations
+import pathlib
+
 import numpy as np
-import scipy.signal as signal
+import scipy
+import scipy.signal
 from scipy.io import wavfile
 from scipy import ndimage
 
 
-def impwav(a):
-    """Imports a wave file as an array where a[1] 
-    is the sampling frequency and a[0] is the data"""
-    wav = wavfile.read(a)
-    out = [wav[1], wav[0]]
-    return out
+def load_wav(wav_path: str | pathlib.Path) -> (int, np.array):
+    """Load .wav file from path.
 
+    Parameters
+    ----------
+    wav_path : str, pathlib.Path
+        Path to a .wav file.
 
-def filtersong(a):
-    """highpass iir filter for song."""
-    out = []
-    b = signal.iirdesign(wp=0.04, ws=0.02, gpass=1, gstop=60, ftype='ellip')
-    out.append(signal.filtfilt(b[0], b[1], a[0]))
-    out.append(a[1])
-    return out
+    Returns
+    -------
+    rate : int
+        Sampling rate in Hz.
+    data : np.ndarray
+        Data from .wav file.
+    """
+    return wavfile.read(wav_path)
 
 
-def smoothrect(a, window=None, freq=None):
-    """smooths and rectifies a song.  Expects (data,samprate)"""
-    if freq == None: freq = 32000  # baseline values if none are provided
-    if window == None: window = 2  # baseline if none are provided
-    le = int(round(freq * window / 1000))  # calculate boxcar kernel length
-    h = np.ones(le) / le  # make boxcar
-    smooth = np.convolve(h, abs(a))  # convovlve boxcar with signal
-    offset = int(round((len(smooth) - len(a)) / 2))  # calculate offset imposed by convolution
-    smooth = smooth[(1 + offset):(len(a) + offset)]  # correct for offset
-    return smooth
+def filtersong(data: np.ndarray) -> np.ndarray:
+    """Apply highpass iir filter to ``data`` to remove low-frequency noise.
+    """
+    b, a = scipy.signal.iirdesign(wp=0.04, ws=0.02, gpass=1, gstop=60, ftype='ellip')
+    return scipy.signal.filtfilt(b, a, data)
 
 
-def getsyls(a):
-    """takes a file read in with impwav and returns a list of sylables"""
-    fa = filtersong(a)  # filter song input
-    frq = a[1]  # get sampling frequency from data (a)
-    a = a[0]  # get data from data (a)
-    frqs = frq / 1000  # calcualte length of a ms in samples
-    objs = findobject(smoothrect(fa[0], 10, frq))  # get syllable positions
-    sylables = [x for x in [a[y] for y in objs] if
-                int(len(x)) > (10 * frqs)]  # get syllable data if of sufficient duration
-    '''uncomment the next line to recover syllables that have been high pass filtered as opposed to raw data.
-    Using data filtered prior to PSD calculation helps if you data are contaminated
-    with low frequency noise'''
-    # sylables=[x for x in [fa[0][y] for y in objs] if int(len(x))>(10*frqs)] #get syllable data if of sufficient duration.
-    objs = [y for y in objs if int(len(a[y])) > 10 * frqs]  # get objects of sufficient duraiton
-    return sylables, objs, frq
-
-
-def threshold(a, thresh=None):
-    """Returns a thresholded array of the same length as input
-    with everything below a specific threshold set to 0.
-    By default threshold is sigma."""
-    if thresh == None: thresh = sc.std(a)
-    out = np.where(abs(a) > thresh, a, np.zeros(a.shape))
-    return out
+def smoothrect(data: np.ndarray,
+               window: int = 2,
+               rate: int = 32000) -> np.ndarray:
+    """Smooth and rectify audio.
+
+    Parameters
+    ----------
+    data : np.ndarray
+        Audio data.
+    window : int
+        Default is 2.
+    rate : int
+        Sampling rate.
+        Default is 32000.
+
+    Returns
+    -------
+    smooth : np.ndarray
+        Smoothed rectified audio.
+    """
+    le = int(round(rate * window / 1000))  # calculate boxcar kernel length
+    h = np.ones(le) / le  # make boxcar
+    smooth = np.convolve(h, abs(data))  # convolve boxcar with signal
+    offset = int(round((len(smooth) - len(data)) / 2))  # calculate offset imposed by convolution
+    smooth = smooth[(1 + offset):(len(data) + offset)]  # correct for offset
+    return smooth
 
 
 # TODO: add option to use scikit-image implementation of Otsu's method
 # https://github.com/NickleDave/songdkl/issues/37
-def findobject(file):
-    """finds objects.  Expects a smoothed rectified amplitude envelope"""
+def findobject(arr: np.ndarray) -> list[tuple[slice]]:
+    """Segment audio into syllables
+    using ``scipy.ndimage.find_objects``.
+    Expects a smoothed rectified amplitude envelope,
+    e.g. as returned by ``songdkl.audio.smoothrect``.
+
+    Parameters
+    ----------
+    arr : numpy.ndarray
+        Containing smoothed rectified amplitude envelope
+        from a .wav file.
+
+    Returns
+    -------
+    objs : list
+        Of tuples of slices;
+        the segments identified by
+        ``scipy.ndimage.find_objects``.
+    """
     # heuristic way of establishing threshold
-    value = (np.average(file)) / 2
-    thresh = threshold(file, value)  # threshold the envelope data
+    value = (np.average(arr)) / 2
+    thresh = threshold(arr, value)  # threshold the envelope data
     thresh = threshold(ndimage.convolve(thresh, np.ones(512)), 0.5)  # pad the threshold
     label = (ndimage.label(thresh)[0])  # label objects in the threshold
     objs = ndimage.find_objects(label)  # recover object positions
     return objs
+
+
+def getsyls(data: np.ndarray,
+            rate: int,
+            min_syl_dur=10,
+            syls_filtered=False) -> tuple[list[np.array], list[tuple[slice]]]:
+    """Return a ``list`` of syllables segmented out of an array of audio.
+
+    Parameters
+    ----------
+    data : np.ndarray
+        Audio data.
+    rate : int
+        Sampling rate, in Hz.
+    min_syl_dur : int
+        Minimum syllable duration, in milliseconds.
+    syls_filtered : bool
+        If True, use audio data to which high-pass filter
+        has been applied. Default is False.
+        Using data filtered prior to PSD calculation helps
+        if data are contaminated with low frequency noise.
+
+    Returns
+    -------
+    syllables : list
+        of ``numpy.ndarray``.
+    slices : list
+        of tuples of slices.
+    """
+    data_filtered = filtersong(data)
+    slices = findobject(smoothrect(data_filtered, 10, rate))
+
+    # get objects of sufficient duration
+    frqs = rate / 1000  # calculate length of a ms in samples
+    # use name ``slice_`` to not clobber ``slice`` function
+    slices = [slice_ for slice_ in slices if data[slice_].shape[0] > min_syl_dur * frqs]
+
+    if syls_filtered:
+        syllables = [x for x in [data_filtered[slice_] for slice_ in slices]]
+    else:
+        syllables = [x for x in [data[slice_] for slice_ in slices]]
+
+    return syllables, slices
+
+
+def threshold(a: np.ndarray, thresh: int | float | None = None) -> np.ndarray:
+    """Returns a thresholded array of the same length as input
+    with everything below a specific threshold set to 0.
+
+    By default threshold is sigma."""
+    if thresh is None:
+        thresh = scipy.std(a)
+    return np.where(abs(a) > thresh, a, np.zeros(a.shape))
diff --git a/src/songdkl/numsyls.py b/src/songdkl/numsyls.py
index cf28919..4906a5e 100644
--- a/src/songdkl/numsyls.py
+++ b/src/songdkl/numsyls.py
@@ -1,111 +1,130 @@
 """functions to estimate the number of syllables in a bird's song"""
-import sys as sys
-import os as os
-import scipy as sc
+import pathlib
 
-from pylab import psd
 import numpy as np
-from sklearn.mixture import GMM
-
-from .audio import impwav, getsyls
-from .songdkl import norm
-
-
-def chunks(l, n):
-    """ Yield successive n-sized chunks from l.
-    """
-    out = []
-    lens = (len(l) / n)
-    for i in range(0, lens):
-        out.append(l[i * n:i * n + n])
-    out.append(l[lens * n:])
-    return out[:-1]
-
-
-def EMofgmmcluster(array_of_syls):
-    """takes an array of segmented sylables, fits a series of
-    gaussian mixture models with an increasing number of mixtures, 
+from sklearn.mixture import GaussianMixture
+import scipy.spatial
+
+from .syllables import (
+    convert_syl_to_psd,
+    get_all_syls,
+)
+
+
+def em_of_gmm_cluster(segedpsds: list,
+                      n_basis: int = 50,
+                      basis: str = 'first',
+                      min_components: int = 2,
+                      max_components: int = 22,
+                      n_splits: int = 1,
+                      ) -> int:
+    """Identifies the best number of mixtures to describe the data.
+
+    Takes an array of segmented syllables, fits a series of
+    gaussian mixture models with an increasing number of mixtures,
     and identifies the best number of mixtures to describe the data
-    by BIC. If data is limiting, we don't recommend using cross validation 
-    for BIC.  If you have sufficient data 3 fold cross validated BIC values
-    compare well with un-cross validated values and human assesment of syllable
-    number.  If data are limiting, cross validated values underestimate
-    the number of syllables relative to human assessment or non-cross
-    validated BIC values.
+    by Bayesian Information Criterion (BIC).
+
+    Parameters
+    ----------
+    segedpsds : list
+        Of ``numpy.ndarray``, PSDs of segmented syllables,
+        returned by ``songdkl.songdkl.convert_syls_to_psd``.
+    n_basis : int
+        Number of syllables to use as basis set. Default is 50.
+    basis : str
+        One of {'first', 'random'}.
+        Controls which syllables are used as the basis set.
+        If 'first', use the first `n_basis` syllables.
+        If `random`, grab a random set of size `n_basis`.
+        Default is 'first'.
+    min_components : int
+        Minimum number of components to consider
+        when fitting Gaussian Mixture Models. Default is 2.
+    max_components : int
+        Maximum number of components to consider
+        when fitting Gaussian Mixture Models. Default is 22.
+    n_splits : int
+        Number of splits to use when estimating components.
+        Default is 1, in which case ``segedpsds`` is not split.
+        If sufficient data are available,
+        a good rule of thumb is to use 3 splits.
+        See Notes below for more detail.
+
+    Returns
+    -------
+    n_syls : int
+        The number of components that gave the minimum
+        Bayesian Information Criterion, plus two.
+
+    Notes
+    -----
+    If data is limiting, we don't recommend using
+    splits to estimate BIC.
+    When the data are limiting,
+    values from splits underestimate
+    the number of syllables relative to
+    human assessment or BIC values
+    computed without splits.
     """
-    nfft = 2 ** 14
-    fs = 32000
-    segstart = int(600 / (fs / nfft))
-    segend = int(16000 / (fs / nfft))
-
-    welchpsds = [psd(x, NFFT=nfft, Fs=fs) for x in array_of_syls]
-    segedpsds1 = [norm(x[0][segstart:segend]) for x in welchpsds]
-
-    models = range(2, 22)
-
-    # select the first 50 syllables of the reference song as the basis set
-    basis_set = segedpsds1[:50]
-    # select a random set of 50 syllablse as the basis set
-    # basis_set=[segedpsds1[rnd.randint(0,len(segedpsds1))] for x in range(50)]
-    d = sc.spatial.distance.cdist(segedpsds1, basis_set, 'sqeuclidean')
-    s = 1 - d / np.max(d) * 1000
-    bic = []
-    # xv=3
-    for x in models:
-        bics = []
-        # this commented section implements cross validation for the BIC values
-        '''ss=chunks(s,len(s)/xv)
-        for y in range(len(ss)):
-        testset=ss[y]
-        trainset=[]
-        [trainset.extend(ss[n]) for n in [z for z in range(len(ss)) if z!=y]]
-        gmm = mixture.GMM(n_components=x, n_iter=100000,n_init=5, covariance_type='full')
-        gmm.fit(np.array(s))
-        bics.append(gmm.bic(np.array(s)))
-        bic.append(np.mean(bics))'''
-        gmm = GMM(n_components=x, n_iter=100000, n_init=5, covariance_type='full')
-        gmm.fit(np.array(s))
-        bic.append(gmm.bic(np.array(s)))
-    return segedpsds1, bic.index(min(bic)) + 2
-
-
-def numsyls(path1):
-    fils1 = [x for x in os.listdir(path1) if x[-4:] == '.wav']
-
-    fils1 = fils1[:120]
-
-    filename1 = path1.split('/')[-2]
-
-    datano = 15000  # maximum number of syllables to be analyzed
-
-    # get syllables from all of the songs in folder
-    syls1 = []
-    objss1 = []
-    for file in fils1:
-        song = impwav(path1 + file)
-        if len(song[0]) < 1: break
-        syls, objs, frq = (getsyls(song))
-        objss1.append(objs)
-        syls1.append([frq] + syls)
-
-    # convert syllables so PSDs
-    segedpsds1 = []
-    for x in syls1[:datano]:
-        fs = x[0]
-        nfft = int(round(2 ** 14 / 32000.0 * fs))
-        segstart = int(round(600 / (fs / float(nfft))))
-        segend = int(round(16000 / (fs / float(nfft))))
-        psds = [psd(norm(y), NFFT=nfft, Fs=fs) for y in x[1:]]
-        spsds = [norm(n[0][segstart:segend]) for n in psds]
-        for n in spsds: segedpsds1.append(n)
-        if len(segedpsds1) > datano: break
-
-    segedpsds, sylno_bic = EMofgmmcluster(segedpsds1)
-
-    print(filename1 + '\t' + str(sylno_bic))
-
+    if basis == 'first':
+        # select the first `n_basis` syllables of the reference song as the basis set
+        basis_set = segedpsds[:n_basis]
+    elif basis == 'random':
+        # select a random set of `n_basis` syllables as the basis set
+        basis_set = [segedpsds[ind]
+                     for ind in np.random.randint(0, len(segedpsds), size=n_basis)]
+
+    D = scipy.spatial.distance.cdist(segedpsds, basis_set, 'sqeuclidean')
+    s = 1 - D / np.max(D) * 1000
+    bics = []
+    n_components_list = list(range(min_components, max_components))
+    for n_components in n_components_list:
+        if n_splits > 1:
+            split_bics = []
+            splits = np.array_split(s, n_splits)
+            for split in splits:
+                gmm = GaussianMixture(n_components=n_components, max_iter=100000, n_init=5, covariance_type='full')
+                gmm.fit(split)
+                split_bics.append(gmm.bic(split))
+            bics.append(np.mean(split_bics))
+        else:
+            gmm = GaussianMixture(n_components=n_components, max_iter=100000, n_init=5, covariance_type='full')
+            gmm.fit(np.array(s))
+            bics.append(gmm.bic(np.array(s)))
+    lowest_bic_ind = np.argmin(bics)
+    n_syls = n_components_list[lowest_bic_ind]
+    return n_syls
+
+
+def numsyls(dir_path,
+            max_wavs: int = 120,
+            max_num_psds: int = 15000,
+            ) -> int:
+    """Determine number of syllable classes in a bird's song,
+    by fitting Gaussian Mixture Models to PSDs of segmented
+    syllable renditions, and selecting the number of components
+    that maximizes the Bayesian Information Criterion.
+
+    Parameters
+    ----------
+    max_wavs : int
+        Maximum number of wav files to use. Default is 120.
+    max_num_psds : int
+        Maximum number of power spectral densities (PSDs) to calculate.
+
+    Returns
+    -------
+    sylno_bic : int
+        The estimated number of syllables, equal to
+        the number of components
+        that produced the fit Gaussian Mixture Model
+        with the lowest Bayesian Information Criterion.
+    """
+    wav_paths = sorted(pathlib.Path(dir_path).glob('*.wav'))
+    wav_paths = wav_paths[:max_wavs]
 
-if __name__ == '__main__':
-    # main program
-    path1 = sys.argv[1]
-    numsyls(path1)
+    syls_from_wavs = get_all_syls(wav_paths)
+    segedpsds = convert_syl_to_psd(syls_from_wavs, max_num_psds)
+    sylno_bic = em_of_gmm_cluster(segedpsds)
+    return sylno_bic
diff --git a/src/songdkl/songdkl.py b/src/songdkl/songdkl.py
index 44cd5bf..2ef95ca 100644
--- a/src/songdkl/songdkl.py
+++ b/src/songdkl/songdkl.py
@@ -1,171 +1,133 @@
 """functions to compute song divergence"""
-import sys
-import os
-from glob import glob
+from __future__ import annotations
+import pathlib
+from typing import Any, Tuple, Union
 
 import scipy.spatial as spatial
-from matplotlib.pylab import psd
 import numpy as np
-from sklearn.mixture import GMM
+from sklearn.mixture import GaussianMixture
 
-from . import audio
+from .syllables import (
+    convert_syl_to_psd,
+    get_all_syls,
+)
 
 
-def norm(a):
-    """normalizes a string by its average and syd"""
-    a=(np.array(a) - np.average(a)) / np.std(a)
-    return a
-
-
-def _get_all_syls(files):
-    """get all syllables from all .wav files in a folder
-
-    Parameters
-    ----------
-    files : list
-        of str, absolute paths to .wav files
-        (output of using glob to find all *.wav files in path)
-
-    Returns
-    -------
-    syls : 
-    objss :
-    """
-    syls = []
-    objss = []
-    for file in files:
-        song = audio.impwav(file)
-        if len(song[0]) < 1:
-            break
-        syls_this_song, objs, frq = audio.getsyls(song)
-        objss.append(objs)
-        syls.append([frq] + syls_this_song)
-    return syls, objss
-
-
-def convert_syl_to_psd(syls, max_num_psds):
-    """convert syllable segments to power spectral density
+def calculate(ref_dir_path: str,
+              compare_dir_path: str,
+              k_ref: int,
+              k_compare: int,
+              max_wavs: int = 120,
+              max_num_psds: int = 10000,
+              n_basis: int = 50,
+              basis:str = 'first') -> Tuple[Union[float, Any], Union[float, Any], int, int]:
+    """Calculate :math:`\text{Song }D_{KL}` metric.
 
     Parameters
     ----------
-    syls : list
-        of ndarray, syllable segments extracted from .wav files
-    max_num_psds : int
-        maximum number of psds to calculate
-
-    Returns
-    -------
-    segedpsds
-    """
-    segedpsds = []
-    for x in syls[:max_num_psds]:
-        fs = x[0]
-        nfft = int(round(2**14/32000.0*fs))
-        segstart = int(round(600/(fs/float(nfft))))
-        segend = int(round(16000/(fs/float(nfft))))
-        psds = [psd(norm(y), NFFT=nfft, Fs=fs) for y in x[1:]]
-        spsds = [norm(n[0][segstart:segend]) for n in psds]
-        for n in spsds:
-            segedpsds.append(n)
-    return segedpsds
-
-
-def calculate(path1, path2, k, k2, max_wavs=120, max_num_psds=10000):
-    """calculate songdkl metric
-
-    Parameters
-    ----------
-    path1 : str
-        path to folder with .wav files of songs from bird 1
-    path2 : str
-        path to folder with .wav files of songs from bird 2
-    k : int
-        number of syllable classes in song 1
-    k2 : int
-        number of syllable classes in song 2
+    ref_dir_path : str
+        Path to directory with .wav files of songs from bird
+        that should be used as reference.
+    compare_dir_path : str
+        Path to directory with .wav files of songs from bird
+        that should be compared with reference.
+    k_ref : int
+        Number of syllable classes in song of bird used as reference.
+    k_compare : int
+        Number of syllable classes in song of bird compared with reference.
     max_wavs : int
-        maximum number of wav files to use. Default is 120.
+        Maximum number of wav files to use. Default is 120.
     max_num_psds : int
-        maximum number of psds to calculate. Default is 10000
+        Maximum number of power spectral densities (PSDs) to calculate.
+        Default is 10000.
+    n_basis : int
+        Number of syllables to use as basis set. Default is 50.
+    basis : str
+        One of {'first', 'random'}.
+        Controls which syllables are used as the basis set.
+        If 'first', use the first `n_basis` syllables.
+        If `random`, grab a random set of size `n_basis`.
+        Default is 'first'.
 
     Returns
     -------
-    None
-
-    prints estimate of song D(KL) to stdout
+    DKL_PQ : float
+        Calculated value for :math:`D_{KL}(\hat{P}||\hat{Q}`,
+        the information lost when encoding P with Q,
+        as computed according to Equation 6 in Mets Brainard 2018.
+        The main value of interest calculated by this function.
+    DKL_QP : float
+        Same computation as ``DKL_PQ``,
+        but in the opposite direction:
+        Q with respect to P.
+    n_psds_ref : int
+        Number of PSDs used from reference data set.
+    n_psd_compare : int
+        Number of PDSs used from comparison data set.
     """
-    fils1 = glob(os.path.join(path1, '*.wav'))
-    fils2 = glob(os.path.join(path2, '*.wav'))
+    wav_paths_ref = sorted(pathlib.Path(ref_dir_path).glob('*.wav'))
+    wav_paths_compare = sorted(pathlib.Path(compare_dir_path).glob('*.wav'))
 
-    fils1 = fils1[:max_wavs]
-    fils2 = fils2[:max_wavs]
+    wav_paths_ref = wav_paths_ref[:max_wavs]
+    wav_paths_compare = wav_paths_compare[:max_wavs]
 
-    syls1, objss1 = _get_all_syls(fils1) 
-    syls2, objss2 = _get_all_syls(fils2)
+    syls_from_wavs_ref = get_all_syls(wav_paths_ref)
+    syls_from_wavs_compare = get_all_syls(wav_paths_compare)
 
-    segedpsds1 = convert_syl_to_psd(syls1, max_num_psds)
-    segedpsds2 = convert_syl_to_psd(syls2, max_num_psds)
+    segedpsds_ref = convert_syl_to_psd(syls_from_wavs_ref, max_num_psds)
+    segedpsds_compare = convert_syl_to_psd(syls_from_wavs_compare, max_num_psds)
 
-    # select the first 50 syllables of the reference song as the basis set
-    basis_set = segedpsds1[:50] 
-    # select a random set of 50 syllablse as the basis set
-    # basis_set=[segedpsds1[rnd.randint(0,len(segedpsds1))] for x in range(50)]
+    if basis == 'first':
+        # select the first `n_basis` syllables of the reference song as the basis set
+        basis_set = segedpsds_ref[:n_basis]
+    elif basis == 'random':
+        # select a random set of `n_basis` syllables as the basis set
+        basis_set = [segedpsds_ref[ind]
+                     for ind in np.random.randint(0, len(segedpsds_ref), size=n_basis)]
+    else:
+        raise ValueError(
+            f"Invalid value for basis: {basis}. Must be one of {{'first', 'random'}}"
+        )
 
-    lone = len(segedpsds1)
-    ltwo = len(segedpsds2)
-    lone_half = int(lone / 2)
-    ltwo_half = int(ltwo / 2)
+    len_ref_half = int(len(segedpsds_ref) / 2)
+    len_compare_half = int(len(segedpsds_compare) / 2)
 
     # calculate distance matrices
-    d1 = spatial.distance.cdist(segedpsds1[:lone_half], basis_set, 'sqeuclidean')
-    d1_2 = spatial.distance.cdist(segedpsds1[lone_half:lone], basis_set, 'sqeuclidean')
-    d2 = spatial.distance.cdist(segedpsds2[:ltwo_half], basis_set, 'sqeuclidean')
-    d2_2 = spatial.distance.cdist(segedpsds2[ltwo_half:ltwo], basis_set, 'sqeuclidean')
+    D_ref = spatial.distance.cdist(segedpsds_ref[:len_ref_half], basis_set, 'sqeuclidean')
+    D_ref_2 = spatial.distance.cdist(segedpsds_ref[len_ref_half:], basis_set, 'sqeuclidean')
+    D_compare = spatial.distance.cdist(segedpsds_compare[:len_compare_half], basis_set, 'sqeuclidean')
+    D_compare_2 = spatial.distance.cdist(segedpsds_compare[len_compare_half:], basis_set, 'sqeuclidean')
 
-    mx = np.max([np.max(d1), np.max(d2), np.max(d1_2), np.max(d2_2)])
+    mx = np.max([np.max(D_ref), np.max(D_compare), np.max(D_ref_2), np.max(D_compare_2)])
 
     # convert to similarity matrices
-    s1 = 1 - (d1 / mx)
-    s1_2 = 1 - (d1_2 / mx)
-    s2 = 1 - (d2 / mx)
-    s2_2 = 1 - (d2_2 / mx)
+    s_ref = 1 - (D_ref / mx)
+    s_ref_2 = 1 - (D_ref_2 / mx)
+    s_compare = 1 - (D_compare / mx)
+    s_compare_2 = 1 - (D_compare_2 / mx)
 
     # estimate GMMs
-    mod1 = GMM(n_components=k, n_iter=100000, n_init=5, covariance_type='full')
-    mod1.fit(s1)
-
-    mod2 = GMM(n_components=k2, n_iter=100000, n_init=5, covariance_type='full')
-    mod2.fit(s2)
+    P = GaussianMixture(n_components=k_ref, max_iter=100000, n_init=5, covariance_type='full')
+    P.fit(s_ref)
 
-    len2 = len(s2)
-    len1 = len(d1)
+    Q = GaussianMixture(n_components=k_compare, max_iter=100000, n_init=5, covariance_type='full')
+    Q.fit(s_compare)
 
     # calculate likelihoods for held out data
-    score1_1 = mod1.score(s1_2)
-    score2_1 = mod2.score(s1_2)
+    p_hat_p = P.score(s_ref_2)
+    q_hat_p = Q.score(s_ref_2)
 
-    score1_2 = mod1.score(s2_2)
-    score2_2 = mod2.score(s2_2)
-
-    len2 = float(len(basis_set))
-    len1 = float(len(basis_set))
+    p_hat_q = P.score(s_compare_2)
+    q_hat_q = Q.score(s_compare_2)
 
     # calculate song divergence (DKL estimate)
-    score1 = np.log2(np.e) * ((np.mean(score1_1)) - (np.mean(score2_1)))
-    score2 = np.log2(np.e) * ((np.mean(score2_2)) - (np.mean(score1_2)))
-
-    score1 = score1 / len1
-    score2 = score2 / len2
+    DKL_PQ = np.log2(np.e) * ((np.mean(p_hat_p)) - (np.mean(q_hat_p)))
+    DKL_QP = np.log2(np.e) * ((np.mean(q_hat_q)) - (np.mean(p_hat_q)))
 
-    # output
-    print(path1 + '\t' + path2 + '\t' + str(k) + '\t' + str(k2) 
-          + '\t' + str(len2) + '\t' + str(score1) + '\t' + str(score2)
-          + '\t' + str(len(segedpsds1)) + '\t' + str(len(segedpsds2)))
+    DKL_PQ = DKL_PQ / len(basis_set)
+    DKL_QP = DKL_QP / len(basis_set)
 
+    n_psds_ref = len(segedpsds_ref)
+    n_psds_compare = len(segedpsds_compare)
 
-if __name__ == '__main__':
-    path1 = sys.argv[1]
-    path2 = sys.argv[2]
-    k = int(sys.argv[3])
-    k2 = int(sys.argv[4])
-    calculate(path1, path2, k, k2)
+    return DKL_PQ, DKL_QP, n_psds_ref, n_psds_compare
diff --git a/src/songdkl/syllables.py b/src/songdkl/syllables.py
new file mode 100644
index 0000000..9f750b2
--- /dev/null
+++ b/src/songdkl/syllables.py
@@ -0,0 +1,128 @@
+"""Functions to segment syllables or other vocalizations
+out of audio, and to convert those segments to
+power spectral densities (PSDs).
+Used by both ``songdkl`` and ``numsyls`` modules.
+"""
+from __future__ import annotations
+import dataclasses
+import pathlib
+from typing import Optional
+
+import matplotlib.mlab
+import numpy as np
+
+from . import audio
+
+
+def norm(arr: np.ndarray) -> np.ndarray:
+    """Normalize an array by
+    subtracting off the mean
+    and dividing by the standard deviation.
+
+    Parameters
+    ----------
+    arr : numpy.ndarray
+        A numpy array
+
+    Returns
+    -------
+    normed : numpy.ndarray
+        ``arr`` with
+    """
+    # replace use of ``np.average`` with ``np.mean``
+    # since ``np.average`` calls ``np.mean`` anyways,
+    # and we are not using weighting that ``np.average`` provides
+    return (np.array(arr) - arr.mean()) / np.std(arr)
+
+
+@dataclasses.dataclass
+class SyllablesFromWav:
+    """
+    Dataclass representing
+    all the syllables
+    found by segmenting a
+    .wav file.
+
+    Attributes
+    ----------
+    syls : list
+        Of ``numpy.ndarray``, syllables
+        returned by ``songdkl.audio.getsyls``
+    rate : int
+        Sampling rate of .wav file.
+    """
+    syls: list
+    rate: int
+
+
+def get_all_syls(wav_paths: list[str] | list[pathlib.Path]) -> list[SyllablesFromWav]:
+    """Get all syllables from a list of .wav files.
+
+    Parameters
+    ----------
+    wav_paths : list
+        Of str, absolute paths to .wav files
+
+    Returns
+    -------
+    syls_from_wavs : list
+        of ``SyllablesFromWav`` instances,
+        same length as ``wav_paths``.
+        Each ``SyllablesFromWav`` has an
+        attribute ``syls``, a list of
+        ``numpy.ndarray``s, and an attribute
+        ``rate``, the sampling rate
+        for the .wav file the syllables
+        are taken from.
+    """
+    syls_from_wavs = []
+    for wav_path in wav_paths:
+        rate, data = audio.load_wav(wav_path)
+        syls_this_wav, slices_this_wav = audio.getsyls(data, rate)
+        syls_from_wavs.append(
+            SyllablesFromWav(syls=syls_this_wav, rate=rate)
+        )
+    return syls_from_wavs
+
+
+def convert_syl_to_psd(syls_from_wavs: list[SyllablesFromWav],
+                       max_num_psds: Optional[int] = None) -> list[np.ndarray]:
+    """Convert syllable segments to power spectral densities (PSDs).
+
+    Parameters
+    ----------
+    syls_from_wavs : list
+        Of ``SyllablesFromWav`` instances,
+        syllable segments extracted from .wav files.
+    max_num_psds : int
+        Maximum number of PSDs to calculate.
+        Default is None, in which case
+        PSDs will be computed for all syllables.
+
+    Returns
+    -------
+    segedpsds : list
+        Of ``numpy.ndarray``,
+        PSDs from segmented syllables.
+    """
+    segedpsds = []
+    for syls_from_wav in syls_from_wavs:
+        fs = syls_from_wav.rate
+        nfft = int(round(2 ** 14 / 32000.0 * fs))
+        segstart = int(round(600 / (fs / float(nfft))))
+        segend = int(round(16000 / (fs / float(nfft))))
+        psds = []
+        for syl in syls_from_wav.syls:
+            Pxx, _ = matplotlib.mlab.psd(norm(syl), NFFT=nfft, Fs=fs)
+            psds.append(Pxx)
+        spsds = [norm(psd[segstart:segend]) for psd in psds]
+        segedpsds.extend(spsds)
+        if max_num_psds:
+            if len(segedpsds) > max_num_psds:
+                break
+
+    if max_num_psds:
+        # since we are likely over `max_num_psds` even after `break` above
+        segedpsds = segedpsds[:max_num_psds]
+
+    return segedpsds


Test patch:
diff --git a/tests/data-for-tests/.gitkeep b/tests/data-for-tests/.gitkeep
new file mode 100644
index 0000000..e69de29
diff --git a/tests/data-for-tests/generated/.gitkeep b/tests/data-for-tests/generated/.gitkeep
new file mode 100644
index 0000000..e69de29
diff --git a/tests/scripts/calc-sylno-pcb.sh b/tests/scripts/calc-sylno-pcb.sh
new file mode 100644
index 0000000..ea43b3d
--- /dev/null
+++ b/tests/scripts/calc-sylno-pcb.sh
@@ -0,0 +1,12 @@
+FIRST_N=$1
+
+PCB_SCRIPT_SYL_COUNTS_FILE=tests/data-for-tests/generated/pcb-script-syl-counts.txt
+
+rm -f $PCB_SCRIPT_SYL_COUNTS_FILE
+touch $PCB_SCRIPT_SYL_COUNTS_FILE
+
+BIRD_FOLDERS=(./tests/data-for-tests/song_data/*/)
+for BIRD_FOLDER in ${BIRD_FOLDERS[@]:1:$FIRST_N};do
+  echo running script for $BIRD_FOLDER
+  conda run -n songdkl-pcb python ./src/pcb-scripts/calc_sylno.py $BIRD_FOLDER >> $PCB_SCRIPT_SYL_COUNTS_FILE
+done
diff --git a/tests/scripts/calc-sylno-py-pkg.sh b/tests/scripts/calc-sylno-py-pkg.sh
new file mode 100644
index 0000000..e9af45a
--- /dev/null
+++ b/tests/scripts/calc-sylno-py-pkg.sh
@@ -0,0 +1,13 @@
+# this script runs the `songdkl` package command `numsyls` on the Bengalese finch song dataset from the PLOS Comp Bio paper
+FIRST_N=$1
+
+PY_PKG_SYL_COUNTS_FILE=tests/data-for-tests/generated/py-pkg-syl-counts.txt
+
+rm -f $PY_PKG_SYL_COUNTS_FILE
+touch $PY_PKG_SYL_COUNTS_FILE
+
+BIRD_FOLDERS=(./tests/data-for-tests/song_data/*/)
+for BIRD_FOLDER in ${BIRD_FOLDERS[@]:1:$FIRST_N};do
+  echo running songdkl for $BIRD_FOLDER
+  songdkl numsyls $BIRD_FOLDER >> $PY_PKG_SYL_COUNTS_FILE
+done