*======================================================================= * affine.f - search for self-affinity in a musical sequence * by Andy Allinger, 2012, released to the public domain * This program may be used by any person for any purpose. * * refer to: * Dimitrios Rafailidis and Yannis Manolopoulos * "The Power of Music: Seaching for Power-Laws in Symbolic Musical Data" * draft report, 2006. * *======================================================================= * I/O LIST * Name Type In/Out Description * note_time[m] REAL in times of note onsets * pitch[m] INTEGER in MIDI note numbers * m INTEGER in how many notes in the sequence * key_time[k] REAL in times of key change events * sharps[k] INTEGER in MIDI key signature * majmin[k] INTEGER in major or minor keys * k INTEGER in # of key changes * beat[m] REAL in phase of each note event * n INTEGER in number of beats * aff REAL out estimated self-affinity * fault INTEGER out error code 0 = no errors * *======================================================================= SUBROUTINE affine(note_time, pitch, m, key_time, sharps, majmin, & k, beat, n, aff, fault) IMPLICIT NONE INTEGER pitch, m, sharps, majmin, k, n, fault REAL note_time, key_time, beat, aff DIMENSION note_time(m), pitch(m), key_time(k), sharps(k), & majmin(k), beat(m) * constants INTEGER L PARAMETER (L = 5) ! length of a segment to compare * local variables INTEGER & D(m), ! pitch, in scale degrees & E(n), ! pitch, at each beat & G(0:11, 0:1), ! translate note # to Do Re Mi & h, i, j, o, ! counters & ind(n-L+1), ! permutation vector for IORDER & JSW(2*(n-L+1)), ! workspace for radix sort & p, ! ocurrences of an identical segment & q, ! number of distinct segments & scale, ! a value of majmin[] & stot(-7:7, 0:1), ! convert keysig to tonic pitch & tm12, ! a value of stot[], tonic - 12 & X(n-L+1, L), ! matrix of segments & xlim ! # bits, largest entry in X[] REAL & ave, ! average of log population & b, ! time of a beat & BIG, ! very large number & dif, ! time from an event to a beat & idzdz, ! external function inverse(Z'/Z) & rmachine, ! external function & near ! shortest value of dif DOUBLE PRECISION sum ! sum of log populations SAVE G, stot DATA G / 0, 1, 2, 3, 4, 6, 7, 8, 9, 10, 11, 12, ! major & 0, 1, 2, 4, 5, 6, 7, 8, 10, 11, 12, 13 / ! minor 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 *======================================================================= * begin *======================================================================= fault = 0 BIG = rmachine(3) IF (m < 10) THEN ! never mind the short stuff aff = 6. ! no affinity on scale 1...6 !!! fault = -1 RETURN END IF *======================================================================= * transpose *======================================================================= h = 1 scale = majmin(1) tm12 = stot(sharps(1),scale) * make invariant to tranpositions within the scale and between keys DO i = 1, m IF (h < k) THEN ! change key IF (key_time(h+1) .LE. note_time(i)) THEN h = h + 1 scale = majmin(h) tm12 = stot(sharps(h), scale) END IF END IF j = pitch(i) - tm12 ! Roman numeral form D(i) = 14 * (j / 12) + G( MOD(j, 12), scale ) END DO *======================================================================= * difference *======================================================================= xlim = 1 h = D(m) DO i = m, 2, -1 j = h h = D(i-1) o = j - h D(i) = o xlim = MAX(xlim, ABS(o)) END DO D(1) = 0 *======================================================================= * quantize, so inserting or deleting brief notes should have no effect *======================================================================= h = 1 i = 1 b = beat(1) sum = 0. DO j = 1, n ! for each beat near = BIG 10 IF (NINT(b) .LE. j) THEN ! is beat j the nearest to event i ? dif = ABS(b - j) IF (dif < near) THEN ! is event i the nearest to beat j ? near = dif h = i END IF IF (i < m) THEN ! more events ? i = i + 1 b = beat(i) GO TO 10 END IF END IF E(j) = D(h) END DO ! next j *======================================================================= * look for repeated patterns *======================================================================= * put segments into matrix DO i = 1, L ! for each pitch in segment DO h = 1, n-L+1 ! for each starting offset X(h,i) = E(i+h-1) END DO END DO * bits to represent largest entry in array i = BIT_SIZE(xlim) DO WHILE (.NOT. BTEST(xlim, i-1)) i = i - 1 END DO xlim = i * sort each column to bring matches together DO i = 1, L DO h = 1, n-L+1 ind(h) = h END DO CALL RSORTI(X(1,i), ind, n-L+1, -xlim, JSW) DO h = 1, L IF (h .EQ. i) CYCLE ! already sorted CALL IORDER(X(1,h), ind, n-L+1) END DO END DO * compare ajacent rows sum = 0. p = 1 q = 1 DO i = 2, n-L+1 DO h = 1, L IF (X(i,h) .NE. X(i-1,h)) GO TO 20 END DO p = p + 1 ! same CYCLE 20 CONTINUE ! different q = q + 1 IF (p .NE. 1) sum = sum + LOG(FLOAT(p)) p = 1 END DO IF (p .NE. 1) sum = sum + LOG(FLOAT(p)) * the Zipf exponent ave = SNGL(sum / DBLE(q)) IF (ave < 0.0223) THEN aff = 5. ELSE IF (ave > 10 000.) THEN aff = 1. ELSE aff = idzdz(-ave) END IF END ! of affine *ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ * inverse of the derivative of the Riemann zeta function divided * by the Riemann zeta function, for use in calculating the Zipf distribution * * an approximation of an approximation of an approximation!, by fitting a * spline to the CERN MATHLIB zeta function, differentiating the spline, * dividing, swapping X and Y, and fitting a quartic rational algebraic * function. * If you know a better way to invert Z'/Z, use it. * * Recommended domain: -2500 < X < -0.025 error < 5% ? * *ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ FUNCTION IDZDZ(X) IMPLICIT NONE REAL IDZDZ, X REAL C(10), numer, denom SAVE C DATA C / 0.2078623159E-03, 0.2160934794, 16.62735147, & -31.33392134, 5.937427891, 0.2078622858E-03, 0.2163012693, & 16.84372205, -14.37591128, 1.000000000 / * synthetic division numer = (((C(1) * X + C(2)) * X + C(3)) * X + C(4)) * X + C(5) denom = (((C(6) * X + C(7)) * X + C(8)) * X + C(9)) * X + C(10) idzdz = numer / denom END ! of IDZDZ *++++++++++++++++++++++++ End of file affine.f +++++++++++++++++++++++++