*----------------------------------------------------------------------- * harmony.f - key and chord music analysis * by Andy Allinger, 2011-2015 * andy_a@users.sourceforge.net * * This program is free software, released to the public domain. * It may be used by any person for any purpose. *----------------------------------------------------------------------- *----------------------------------------------------------------------- * segment - split a note list into minimal partitions *----------------------------------------------------------------------- * * I/O LIST *_Name______________Type________In/Out____Description___________________ * note_time[n] REAL In times of note on/off events * (assumed already sorted) * pitch[n] INTEGER In MIDI note # * offon[n] INTEGER In 0=OFF, 1=ON * n INTEGER In number of note events * pall INTEGER Out # of minimal partitions * partit[n] REAL Out times of partitions * seg[0:11,0:n] INTEGER Out notes each pitch class, each segment * cumseg[0:11,0:n] INTEGER Out cumulative seg[] * fault INTEGER Out error code, 0=success * *----------------------------------------------------------------------- SUBROUTINE SEGMENT (note_time, pitch, offon, n, & pall, partit, seg, cumseg, fault) IMPLICIT NONE INTEGER pitch, offon, n, pall, seg, cumseg, fault REAL note_time, partit DIMENSION note_time(n), pitch(n), offon(n), partit(n), & seg(0:11,n), cumseg(0:11,n) * local variables INTEGER & i, j, k, ! counters & npc, ! pitch % 12 & state(0:11) ! pitch vector REAL & pip, ! very brief duration & t, tprev ! event times PARAMETER (pip = 0.035) ! negligible time *----------------------------------------------------------------------- * initialize, validate *----------------------------------------------------------------------- fault = 0 ! clear error codes IF (n < 2) THEN fault = 1 ! at least one NoteOff and NoteOn RETURN END IF *----------------------------------------------------------------------- * Find minimal segments -initialize states, times, and partition times * seg[i] the pitch vector for the segment beginning at partit[i] * cumseg[i] the accumulation of seg up to BUT NOT INCLUDING seg[i] * partit[1] time of the first note on * partit[pall] time of the last note off * In the case of non-overlapping monophony, n = pall *----------------------------------------------------------------------- DO j = 0, 11 state(j) = 0 END DO *----------------------------------------------------------------------- * consider each NoteOff or NoteOn as a possible partition point *----------------------------------------------------------------------- tprev = note_time(1) ! indifferent to beginning silence t = tprev ! avoid compiler warning i = 0 DO j = 1, n npc = MOD(pitch(j), 12) t = note_time(j) IF (t - tprev > pip) THEN ! different time? i = i + 1 DO k = 0, 11 seg(k, i) = state(k) END DO partit(i) = tprev tprev = t END IF IF (offon(j) .EQ. 1) THEN state(npc) = state(npc) + 1 ! increment pitch vector ELSE IF (offon(j) .EQ. 0) THEN state(npc) = state(npc) - 1 ! decrement pitch vector ELSE fault = 2 ! bad input RETURN END IF END DO * last note off (bogus silent segment) i = i + 1 pall = i partit(i) = t DO j = 0, 11 seg(j,i) = state(j) IF (state(j) .NE. 0) PRINT *, 'segment: note sounding at end!!' END DO *----------------------------------------------------------------------- * Accumulate the pitch vectors. *----------------------------------------------------------------------- DO j = 0, 11 cumseg(j,1) = 0 END DO DO i = 2, pall ! over all minimal segments DO j = 0, 11 cumseg(j,i) = cumseg(j,i-1) + seg(j,i-1) END DO END DO RETURN END ! of SEGMENT *XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX * keyal - find key changes using the chi-squared test and cross-entropy *XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX * * TRANSFER LIST *__Name________________Type_____I/O___Description_______________________ * partit[n] REAL In times of minimal partitions * cumseg[0:11,n] INTEGER In cumulative pitch vector * n INTEGER In length of partit and cumseg * key_t[n] REAL Out times of key changes * sharps[n] INTEGER Out # sharps in key signature * majmin[n] INTEGER Out major or minor key * k INTEGER Out # key changes detected * plist[n] INTEGER Work 0=Not a partition, -1=is a partition * cumw[n] INTEGER Work sum of cumseg each partition * uprofile[0:11,0:1] REAL In custom pitch profile * usecfg LOGICAL In user has set a custom profile * fault INTEGER Out 0=no errors * *XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX SUBROUTINE KEYAL (partit, cumseg, n, key_t, sharps, majmin, k, & plist, cumw, uprofile, usecfg, fault) IMPLICIT NONE INTEGER cumseg, n, sharps, majmin, k, plist, cumw, fault REAL partit, key_t, uprofile LOGICAL usecfg DIMENSION partit(n), cumseg(0:11,n), key_t(n), sharps(n), & majmin(n), plist(n), cumw(n), uprofile(0:11,0:1) * local variables INTEGER & ba(0:11,2), ! before-after contingency table & i0, i, j, ! counters & iearly(32), ilate(32), ! indices of partition points & L, m, ! count NPC's (Neutral Pitch Classes), modes & Ldir, ! -1 descending tree, +1 ascending tree & root, scale, ! detected root and scale & wmin, ! limit for performing chi-square test & wt(2) ! cumw before and after proposed split REAL & default, ! small probability & e, ! expected contingency cell value & rmachine, ! external function & obs(0:11), ! real-valued pitch vector & profile(0:11,0:11,0:1), ! key-finding template & s, sbest, ! entropy & tprev, & wtot, ! total of contingency table entries & xcrit, ! cutoff value for chi-square test & x2, x2big ! chi-square statistic LOGICAL first SAVE first, profile PARAMETER (xcrit = 33.0, ! p > .9995 with 11 degrees of freedom & wmin = 60, ! 5 objects x 12 cells & default = 0.000 001) ! 1 millionth DATA first / .TRUE. / *XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX * The Kostka-Payne key profiles * * Refer to: * David Temperley * Music and Probability, ch. 6: "A Polyphonic Key-Finding Model" * MIT Press, 2007 * * profile[npc, root, mode] *XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX DATA profile( 0,0,0) / 0.748 / DATA profile( 1,0,0) / 0.060 / DATA profile( 2,0,0) / 0.488 / DATA profile( 3,0,0) / 0.082 / DATA profile( 4,0,0) / 0.670 / DATA profile( 5,0,0) / 0.460 / DATA profile( 6,0,0) / 0.096 / DATA profile( 7,0,0) / 0.715 / DATA profile( 8,0,0) / 0.104 / DATA profile( 9,0,0) / 0.366 / DATA profile(10,0,0) / 0.057 / DATA profile(11,0,0) / 0.400 / DATA profile( 0,0,1) / 0.712 / DATA profile( 1,0,1) / 0.084 / DATA profile( 2,0,1) / 0.474 / DATA profile( 3,0,1) / 0.618 / DATA profile( 4,0,1) / 0.049 / DATA profile( 5,0,1) / 0.460 / DATA profile( 6,0,1) / 0.105 / DATA profile( 7,0,1) / 0.747 / DATA profile( 8,0,1) / 0.404 / DATA profile( 9,0,1) / 0.067 / DATA profile(10,0,1) / 0.133 / DATA profile(11,0,1) / 0.330 / *XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX * begin *XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX IF (n < 2) THEN fault = -1 RETURN ELSE fault = 0 END IF * prepare key profiles IF (first) THEN first = .FALSE. IF (usecfg) THEN ! copy user's data DO m = 0, 1 DO i = 0, 11 IF (uprofile(i,m) > 0.) THEN profile(i,0,m) = uprofile(i,m) ELSE profile(i,0,m) = default END IF END DO END DO END IF * log probability DO m = 0, 1 DO i = 0, 11 profile(i,0,m) = LOG(profile(i,0,m)) END DO END DO * transpose profiles DO m = 0, 1 ! major, minor DO j = 1, 11 ! roots 1...11 DO i = 0, 11 ! each pitch class profile(i,j,m) = profile(MOD(i+11,12), j-1, m) END DO END DO END DO END IF ! first ? *XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX * initial values *XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX DO i = 1, n L = 0 DO j = 0, 11 L = L + cumseg(j,i) END DO cumw(i) = L plist(i) = 0 END DO *XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX * recursively partition *XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX L = 1 ! at trunk of tree Ldir = +1 iearly(L) = 1 ilate(L) = n plist(1) = -1 plist(n) = -1 DO WHILE (L > 0) ! until back down from tree IF (Ldir < 0) THEN ! already considered this area ? IF (iearly(L) .EQ. iearly(L+1)) THEN ! came from the early side? iearly(L+1) = ilate(L+1) ! do the late sub-region ilate(L+1) = ilate(L) L = L + 1 ! ascend Ldir = +1 ELSE ! came from late sub-region L = L - 1 ! continue slide END IF ELSE ! find optimal split x2big = 0. i0 = 1 ! avoid compiler warning DO 40 i = iearly(L)+1, ilate(L)-1 wt(1) = cumw(i) - cumw(iearly(L)) wt(2) = cumw(ilate(L)) - cumw(i) wtot = FLOAT(wt(1) + wt(2)) IF (wt(1) < wmin .OR. wt(2) < wmin) GO TO 40 ! can't apply test DO m = 1, 2 ! clear contingency table DO j = 0, 11 ba(j,m) = 0 END DO END DO DO j = 0, 11 ! count ba(j,1) = cumseg(j,i) - cumseg(j,iearly(L)) END DO DO j = 0, 11 ba(j,2) = cumseg(j,ilate(L)) - cumseg(j,i) END DO * chi-square x2 = 0. DO m = 1, 2 DO j = 0, 11 e = FLOAT(wt(m) * (ba(j,1) + ba(j,2))) / wtot ! expected count x2 = x2 + (ba(j,m) - e) **2 / e END DO END DO IF (x2 > x2big) THEN x2big = x2 ! maximum chi-square i0 = i ! best split point END IF 40 END DO ! next i IF (x2big > xcrit) THEN ! split? plist(i0) = -1 ! true iearly(L+1) = iearly(L) ! consider subdividing this part also ilate(L+1) = i0 L = L + 1 ! ascend ELSE ! leave intact Ldir = -1 L = L - 1 ! descend END IF END IF ! ascend or descend ? END DO ! while loop *----------------------------------------------------------------------- * Simple key-finding based on the Krumhansl-Schmuckler algorithm. * Cross-entropy is used instead of the correlation coefficient. *----------------------------------------------------------------------- i0 = 1 k = 0 tprev = 0. DO i = 2, n IF (plist(i) .EQ. 0) GO TO 60 ! not a partition * get an observation vector DO j = 0, 11 obs(j) = FLOAT(cumseg(j,i) - cumseg(j,i0)) END DO *----------------------------------------------------------------------- * measure cross-entropy S = -SUM (No[i] @ ln(Pm[i]) * where No observed frequency * Pm model probability *----------------------------------------------------------------------- sbest = rmachine(3) ! very large number root = 0 ! avoid uninitialized warnings scale = 0 DO m = 0, 1 ! major, minor DO L = 0, 11 ! each root s = 0. DO j = 0, 11 ! each pitch class s = s - obs(j) * profile(j,L,m) END DO IF (s < sbest) THEN sbest = s root = L scale = m END IF END DO END DO * convert root to number of sharps in key signature j = MOD(7 * root + 5, 12) - 5 IF (k > 0) THEN ! check for same key IF (j .EQ. sharps(k) .AND. scale .EQ. majmin(k)) GO TO 60 END IF k = k + 1 key_t(k) = tprev sharps(k) = j majmin(k) = scale i0 = i tprev = partit(i) 60 END DO ! next i RETURN END ! of KEYAL *----------------------------------------------------------------------- * chordal - detect chord changes * * based on the algorithm published in: * Bryan Pardo & William P. Birmingham, * "The Chordal Analysis of Tonal Music" * University of Michigan technical report, 28 March 2001 * * and: * Bryan Pardo & William P. Birmingham * "Algorithms for Chordal Analysis" * Computer Music Journal, v. 26 n. 2 pp. 27-49, Summer 2002 * *----------------------------------------------------------------------- * * I/O LIST *_Name______________Type________In/Out____Description___________________ * partit[n] Real In times of minimal partitions * seg[0:11,n] Integer In pitch vector each partition * cumseg[0:11,n] Integer In cumulative seg[] * n Integer In length of partit, seg, and cumseg * key_time[k] Real In times of key change events * (assumed already sorted) * sharps[k] Integer In # sharps in key signature {-7...7} * (matches MIDI representation) * majmin[k] Integer In 0=major key, 1=minor key * k Integer In number of key change events * job Character In {L, T, E} * L => label chords * T => label and train * E => label and evaluate * c Integer Out # of chord segments detected * ch_t[n] Real Out times of chord changes [seconds] * q[n] Integer Out "qualities" of chords * (as in the templates, see below) * r[n] Integer Out roots of chords {0...11} * as pitch % 12 * p Real Out log probability of chord sequence * ccarr[num_q,0:11,0:11,0:1] * Real In log probability of chord given chord * in Roman numeral notation. The * table is given twice, once * for major keys, and again for * minor keys. num_q is the number * of chord qualities. A separate * routine is needed to define ccarr[]. * prox[0:6] Real Out circle-of-fifth proximity of roots * qual[num_q, 0:1] Real Out count of each quality, each mode * root[0:11, 0:1] Real Out count of each root, each mode * fault Integer Out error code, 0=no errors * *----------------------------------------------------------------------- SUBROUTINE CHORDAL (partit, seg, cumseg, n, key_time, sharps, & majmin, k, job, c, ch_t, q, r, p, ccarr, prox, qual, root, fault) IMPLICIT NONE INTEGER num_q PARAMETER (num_q = 15) ! how many chord qualities INTEGER seg, cumseg, n, sharps, majmin, k, c, q, r, fault REAL partit, key_time, ch_t, p, ccarr, prox, qual, root CHARACTER job DIMENSION partit(n), seg(0:11,n), cumseg(0:11,n), & key_time(k), sharps(k), majmin(k), ch_t(n), q(n), r(n), & ccarr(num_q, 0:11, 0:11, 0:1), & prox(0:6), qual(num_q, 0:1), root(0:11, 0:1) * local variables INTEGER & C5(0:23), ! half-steps to circle of fifths & hsteps(0:23), ! shortcut for % 12 & h, i, j, ! counters & prev, pprev, ! index of prev seg, past prev & qdel, qkeep, qold, qnew, ! qualities & rdel, rkeep, rold, rnew, ! roots of chords & romold, romnew, ! Roman Numeral version of roots & rpp, ! an earlier root value & scale, ! 0=major or 1=minor & skip, & stot(-7:7,0:1), ! convert sharps to tonic -12...-1 & template(0:11, num_q), ! chord templates & tm12, ! tonic - 12 & weight(0:11) ! pitch vector to analyze REAL & sdel, skeep, sprev ! possible segment scores LOGICAL & ALMEQ, ! an almost-equal function & shifted ! whether templates have been transposed SAVE hsteps, C5, shifted, stot, template *----------------------------------------------------------------------- * chord templates: bits 0...11 of an integer represent pitches as * MIDI note number % 12 * If a note is present in the chord, the corresponding bit is set * in the template. Template data is included for a root note of C * and transposed to the other keys. For example, C Major has the notes * C, E, G and the template {0, 4, 7}, which as an integer is: * 2^0 + 2^4 + 2^7 = 145 * * The order of the array is significant, lower numbered qualities * have precedence in case of ties during template-matching. * * It is possible to redefine which templates are included. There is no * authority as to what is or is not a real chord. This list should be * adequate for most purposes. *----------------------------------------------------------------------- * integer name template symbol DATA template(0, 1) / 145/ ! major {0,4,7} C DATA template(0, 2) / 137/ ! minor {0,3,7} Cm DATA template(0, 3) /1169/ ! seventh {0,4,7,10} C7 DATA template(0, 4) /2193/ ! major seventh {0,4,7,11} Cmaj7 DATA template(0, 5) /1161/ ! minor seventh {0,3,7,10} Cm7 DATA template(0, 6) / 657/ ! sixth {0,4,7,9} C6 DATA template(0, 7) / 649/ ! minor sixth {0,3,7,9} Cm6 DATA template(0, 8) / 73/ ! diminished {0,3,6} Cdim DATA template(0, 9) / 273/ ! augmented {0,4,8} Caug DATA template(0,10) / 585/ ! diminished seventh {0,3,6,9} Cdim7 DATA template(0,11) /1105/ ! seventh flat fifth {0,4,6,10} C7-5 DATA template(0,12) /1185/ ! seventh suspended fourth {0,5,7,10} C7sus4 DATA template(0,13) /1173/ ! ninth {0,2,4,7,10} C9 DATA template(0,14) /1205/ ! eleventh {0,2,4,5,7,10} C11 DATA template(0,15) /1685/ ! thirteenth {0,2,4,7,9,10} C13 DATA shifted /.FALSE./ DATA C5 / 0,5,2,3,4,1,6,1,4,3,2,5,0,5,2,3,4,1,6,1,4,3,2,5 / DATA hsteps /0,1,2,3,4,5,6,7,8,9,10,11,0,1,2,3,4,5,6,7,8,9,10,11/ DATA stot / -1,-6,-11,-4,-9,-2,-7,-12,-5,-10,-3,-8,-1,-6,-11, ! major & -4,-9,-2,-7,-12,-5,-10,-3,-8,-1,-6,-11,-4,-9,-2 / ! minor *----------------------------------------------------------------------- * initialize, validate *----------------------------------------------------------------------- fault = 0 ! clear error codes * transpose chord templates into all keys IF (.NOT. shifted) THEN DO i = 1, num_q DO j = 1, 11 template(j,i) = ISHFTC(template(j-1,i), 1, 12) END DO END DO shifted = .TRUE. END IF * check for valid job IF (job .NE. 'L' .AND. job .NE. 'E' .AND. job .NE. 'T') THEN fault = 3 RETURN END IF *----------------------------------------------------------------------- * HarmAn search. Consider each partition point in order, as * a possible point for a chord change, based on whether keeping the * partition or deleting it yields a better match to a chord template. *----------------------------------------------------------------------- h = 1 tm12 = stot(sharps(h), majmin(h)) prev = 0 DO j = 0, 11 weight(j) = seg(j,1) END DO CALL score_seg (weight, template, tm12, prev, 0, & q(1), r(1), sprev) *----------------------------------------------------------------------- * main loop *----------------------------------------------------------------------- pprev = 0 prev = 1 rpp = 0 DO i = 2, n-1 IF (h < k) THEN ! change key IF (partit(i) .GE. key_time(h+1)) THEN h = h + 1 tm12 = stot(sharps(h), majmin(h)) END IF END IF *----------------------------------------------------------------------- * score segment i *----------------------------------------------------------------------- DO j = 0, 11 weight(j) = seg(j,i) END DO CALL score_seg (weight, template, tm12, prev, r(prev), & qkeep, rkeep, skeep) *----------------------------------------------------------------------- * score segment formed by deleting partition i *----------------------------------------------------------------------- DO j = 0, 11 weight(j) = cumseg(j,i+1) - cumseg(j,prev) END DO IF (pprev > 0) rpp = r(pprev) CALL score_seg (weight, template, tm12, pprev, rpp, & qdel, rdel, sdel) *----------------------------------------------------------------------- * keep or delete partition point ? *----------------------------------------------------------------------- IF (skeep + sprev > sdel) THEN ! keep partition q(i) = qkeep r(i) = rkeep sprev = skeep pprev = prev prev = i ELSE ! delete partition q(i) = -1 ! delete labeling r(i) = -1 ! " " q(prev) = qdel r(prev) = rdel sprev = sdel END IF END DO *----------------------------------------------------------------------- * remove values of q[] and r[] at deleted partition points *----------------------------------------------------------------------- c = 1 qold = q(1) rold = r(1) DO 30 i = 2, n-1 IF (q(i) < 0) GO TO 30 ! not a partition point qnew = q(i) rnew = r(i) IF ( qold .EQ. qnew .AND. rold .EQ. rnew ) GO TO 30 ! skip duplicates c = c + 1 ! copy data to lower position ch_t(c) = partit(i) q(c) = qnew r(c) = rnew qold = qnew rold = rnew 30 END DO ! next i IF (job .EQ. 'L') RETURN IF (job .EQ. 'T') GO TO 70 ! count chord change probabilities *----------------------------------------------------------------------- * If JOB is set to evaluate, find probability of the sequence of chords, * figured as a Markov chain of events E[1]...E[n] : * P(E) = P(E[n]) * PROD ( P(E[i-1] | E[i]) ) 2...i...n *----------------------------------------------------------------------- p = 0. h = 1 rold = r(1) scale = majmin(1) tm12 = stot(sharps(1), scale) romold = hsteps(r(1) - tm12) DO i = 2, c IF (h < k) THEN ! change key IF (ch_t(i) .GE. key_time(h+1)) THEN h = h + 1 scale = majmin(h) tm12 = stot(sharps(h), scale) END IF END IF qnew = q(i) rnew = r(i) romnew = hsteps(rnew - tm12) * Markov chain of log-probability p = p + ccarr(qnew, romnew, romold, scale) rold = rnew romold = romnew END DO RETURN ! Done evaluating probability *----------------------------------------------------------------------- 70 CONTINUE ! count chord proximity, quality and root *----------------------------------------------------------------------- h = 1 rold = r(1) scale = majmin(1) tm12 = stot(sharps(1), scale) romold = hsteps(r(1) - tm12) skip = 0 DO i = 2, c IF (h < k) THEN ! change key IF (ch_t(i) .GE. key_time(h+1)) THEN h = h + 1 scale = majmin(h) tm12 = stot(sharps(h),scale) IF (ALMEQ(ch_t(i), key_time(h))) THEN skip = 1 ! skip one transition on a neat key change ELSE skip = 2 ! skip two transitions on a mid-chord key change END IF END IF END IF qnew = q(i) rnew = r(i) romnew = hsteps(rnew - tm12) IF (skip .LE. 0) THEN j = C5(romnew - romold + 12) ! proximity prox(j) = prox(j) + 1. qual(qnew,scale) = qual(qnew,scale) + 1. ! quality root(romnew,scale) = root(romnew,scale) + 1. ! root END IF rold = rnew romold = romnew skip = skip - 1 END DO RETURN END ! of CHORDAL *----------------------------------------------------------------------- * score_seg - score a segment * *__Name__________________Type______In/Out___Description_________________ * weight[0:11] Integer In observed pitch vector * template[0:11,num_q] Integer In chord template data * tm12 Integer In tonic key as NPC, minus 12 (!) * prev Integer In index of previous segment * rprev Integer In root chosen at previous segment * qbest Integer Out quality this segment * rbest Integer Out root this segment * sbest Real Out score this segment labeling * bigger numbers are better *----------------------------------------------------------------------- SUBROUTINE score_seg (weight, template, tm12, prev, rprev, & qbest, rbest, sbest) IMPLICIT NONE INTEGER num_q PARAMETER (num_q = 15) INTEGER weight, template, tm12, prev, rprev, qbest, rbest REAL sbest DIMENSION weight(0:11), template(0:11,num_q) * local variables INTEGER C5(0:23), ! half steps to circle of fifths dist. & ctPP, ctPA, ctAP, ctAA, ! contingency table entries & h, i, j, ! count qualities, roots, pitches & o, n, ! temporary integers & wsum ! total of weight[] REAL rmachine, score SAVE C5 DATA C5 / 0,5,2,3,4,1,6,1,4,3,2,5,0,5,2,3,4,1,6,1,4,3,2,5 / *----------------------------------------------------------------------- wsum = 0 DO j = 0, 11 wsum = wsum + weight(j) END DO IF (wsum .EQ. 0) THEN qbest = 1 rbest = 0 sbest = 0. RETURN END IF sbest = -rmachine(3) ! real-valued scores DO h = 1, num_q ! for each quality DO i = 0, 11 ! for each root score = 0. ctPP = 0 ctPA = 0 ctAP = 0 ctAA = 0 DO j = 0, 11 ! for each pitch class IF ( BTEST(template(i,h), j) ) THEN IF (weight(j) > 0) THEN ctPP = ctPP + weight(j) ! in template ELSE !! ctPA = ctPA + wsum ! don't know why, but it helps... ctPA = ctPA + 1 END IF ELSE ! not in template IF (weight(j) > 0) THEN ctAP = ctAP + weight(j) ELSE !! ctAA = ctAA + wsum ! ...to add wsum rather than 1 ! ctAA = ctAA + 1 END IF END IF END DO *----------------------------------------------------------------------- * contingency table formula * Don't know why it's the fourth root of the denominator, but it helps *----------------------------------------------------------------------- score = FLOAT(ctPP * ctAA - ctPA * ctAP) * SQRT(FLOAT(wsum)) & / SQRT(SQRT(FLOAT( & (ctPP + ctPA) * (ctPA + ctAA) * (ctAA + ctAP) * (ctAP + ctPP) ))) IF (score - sbest) 50, 30, 40 30 CONTINUE ! tiebreakers * prefer more root pitch IF (weight(i) > weight(rbest)) THEN qbest = h rbest = i sbest = score GO TO 50 ELSE IF (weight(i) < weight(rbest)) THEN GO TO 50 END IF * prefer clearer harmonies IF (h < qbest) THEN qbest = h rbest = i sbest = score GO TO 50 ELSE IF (h > qbest) THEN GO TO 50 END IF * nearest root on the circle of fifths IF (prev > 0) THEN o = C5(rbest - rprev + 12) ! old n = C5(i - rprev + 12) ! new IF (n < o) THEN qbest = h rbest = i sbest = score GO TO 50 ELSE IF (n > o) THEN GO TO 50 END IF END IF * prefer roots near to the tonic on the circle of fifths o = C5(rbest - tm12) ! old n = C5(i - tm12) ! new IF (n < o) THEN qbest = h rbest = i sbest = score GO TO 50 ELSE IF (n > o) THEN GO TO 50 END IF * resolve aug, dim7, 7b5 based on prior chord * This rule is intended to make Bdim7 follow C IF ( (h .EQ. 9) .OR. (h .EQ. 10) .OR. (h .EQ. 11) & .AND. (prev > 0) ) THEN o = MOD(rprev - rbest + 12, 12) ! half steps DESCENDED n = MOD(rprev - i + 12, 12) IF (n < o) THEN qbest = h rbest = i sbest = score GO TO 50 ELSE IF (n > o) THEN GO TO 50 END IF END IF ! dim7 resolution ! PRINT *, 'chordal: warning - unbreakable tie: ', ! & h, i, ' vs ', qbest, rbest GO TO 50 40 CONTINUE ! new preferred chord qbest = h rbest = i sbest = score 50 END DO ! next i END DO ! next h RETURN END ! of score_seg ************************************************************************ * finish making the user custom chord profile ************************************************************************ SUBROUTINE MKPQR (prox, qual, root) IMPLICIT NONE INTEGER num_q PARAMETER (num_q = 15) REAL prox, qual, root DIMENSION prox(0:6), qual(num_q,0:1), root(0:11,0:1) INTEGER i, L REAL seventh, twelfth, tot LOGICAL safdiv ! safe-to-divide function SAVE seventh, twelfth DATA seventh, twelfth / 0.14285714, 0.083333333 / * proximity tot = 0. DO i = 0, 6 tot = tot + prox(i) END DO DO i = 0, 6 IF (SAFDIV(prox(i), tot)) THEN prox(i) = prox(i) / tot ELSE prox(i) = seventh END IF END DO * quality DO L = 0, 1 tot = 0. DO i = 1, num_q tot = tot + qual(i,L) END DO DO i = 1, num_q IF (SAFDIV(qual(i,L), tot)) THEN qual(i,L) = qual(i,L) / tot ELSE qual(i,L) = twelfth END IF END DO END DO * root DO L = 0, 1 tot = 0. DO i = 0, 11 tot = tot + root(i,L) END DO DO i = 0, 11 IF (SAFDIV(root(i,L), tot)) THEN root(i,L) = root(i,L) / tot ELSE root(i,L) = twelfth END IF END DO END DO RETURN END ! of MKPQR *----------------------------------------------------------------------- * Make a chord change probability matrix * The chord data here has been collected from 10 tunes, 5 in major keys * and 5 in minor keys, selected at whim from: * Songs of the Dragon * WNGGA, 1993. * * TRANSFER LIST *_Name______________Type_____I/O__Description______________________ * uprox[0:6] Real In probability circle-of-fifths root proximity * uqual[num_q, 0:1] Real In probability each chord quality, each mode * uroot[0:11, 0:1] Real In probability of Roman numeral root, each mode * usecfg Logical In user supplies values for prox, qual, root * * ccarr[num_q, 0:11, 0:11, 0:1] * Real Out log probability each successive quality, * each succesive (Roman numeral) root; * given each preceding root, * for major and minor modes *----------------------------------------------------------------------- SUBROUTINE MKCCARR (uprox, uqual, uroot, usecfg, ccarr) IMPLICIT NONE INTEGER num_q PARAMETER (num_q = 15) REAL ccarr, uqual, uprox, uroot LOGICAL usecfg DIMENSION ccarr(num_q, 0:11, 0:11, 0:1), & uprox(0:6), uqual(num_q, 0:1), uroot(0:11, 0:1) * local variables INTEGER i, j, k, m, ! counters & C5(0:23), ! circle of fifths distance & dif REAL & default, ! arbitrary small probability & prox(0:6), ! prob. of root proximity & qual(num_q, 0:1), ! prob. of each quality, each mode & root(0:11, 0:1) ! prob. of each root, each mode DOUBLE PRECISION tot SAVE C5, prox, qual, root DATA C5 / 0,5,2,3,4,1,6,1,4,3,2,5,0,5,2,3,4,1,6,1,4,3,2,5 / * Casual inspection showed no difference between major and minor modes, * so the data was pooled. A Cauchy distribution with average 1.161 * and dispersion .471 extrapolated the rare/unobserved values. DATA prox / .102, .644, .173, .044, .019, .011, .007 / * The quality probabilities also take into account chords noted in * the Kostka-Payne data set, as reported by Pardo & Birmingham. * This array must be arranged to match the chord templates! DATA qual / .600, .056, .222, .003, .003, .003, .035, .019, ! Major & .003, .041, .003, .003, .003, .003, .003, & .186, .526, .134, .003, .003, .003, .068, .019, & .003, .041, .003, .003, .003, .003, .003 / * probability of each root, nominal value for unobserved roots DATA root / 0.387, 0.001, 0.114, 0.001, 0.001, 0.114, ! major & 0.001, 0.334, 0.001, 0.030, 0.015, 0.001, & 0.441, 0.001, 0.015, 0.052, 0.001, 0.150, ! minor & 0.001, 0.269, 0.007, 0.001, 0.060, 0.001 / DATA default / 0.000 001 / ! one in a million *----------------------------------------------------------------------- * multiply *----------------------------------------------------------------------- IF (usecfg) THEN ! use defaults DO i = 0, 6 prox(i) = uprox(i) END DO DO j = 0, 1 DO i = 1, num_q qual(i,j) = uqual(i,j) END DO END DO DO j = 0, 1 DO i = 0, 11 root(i,j) = uroot(i,j) END DO END DO END IF ! default chord data? DO m = 0, 1 ! major and minor modes DO k = 0, 11 ! each prior root DO j = 0, 11 ! each present root dif = C5(ABS(j - k)) DO i = 1, num_q ! each present quality ccarr(i,j,k,m) = prox(dif) * qual(i,m) * root(j,m) END DO END DO END DO END DO *----------------------------------------------------------------------- * normalize *----------------------------------------------------------------------- DO m = 0, 1 ! each mode DO k = 0, 11 ! each prior root * add tot = 0. ! sum over priors DO j = 0, 11 ! each present root DO i = 1, num_q ! each present quality tot = tot + ccarr(i,j,k,m) END DO END DO IF (tot .LE. 0.) tot = default * divide and take logarithm DO j = 0, 11 DO i = 1, num_q ccarr(i,j,k,m) = SNGL(DBLE(ccarr(i,j,k,m)) / tot) IF (ccarr(i,j,k,m) > default) THEN ccarr(i,j,k,m) = LOG(ccarr(i,j,k,m)) ELSE ccarr(i,j,k,m) = LOG(default) END IF END DO END DO END DO END DO RETURN END ! of mkccarr *+++++++++++++++++++++++ End of file harmony.f +++++++++++++++++++++++++