*----------------------------------------------------------------------- * puntal.f - contrapuntal analysis * by Andy Allinger, 2012, released to the public domain * This program may be used by any person for any purpose. *----------------------------------------------------------------------- *----------------------------------------------------------------------- * edist - distance between notes, such that they belong to different voices * * refer to : * Wai Man Szeto and Man Hon Wong * "Stream Segregation Algorithm for Pattern Matching * in Polyphonic Music Databases," * draft report, 2003. * *___Name_______Type______In/Out___Description_________________ * TON(N) Real In Note ontimes * TOFF(N) Real In Note offtimes * PITCH(N) Integer In MIDI note # * N Integer In Length of arrays * I Integer In Index of first point * J Integer In Index of second point * EDIST Real Out Distance *----------------------------------------------------------------------- FUNCTION EDIST (TON, TOFF, PITCH, N, I, J) IMPLICIT NONE INTEGER PITCH, N, I, J REAL TON, TOFF, EDIST DIMENSION TON(N), TOFF(N), PITCH(N) REAL BETA PARAMETER (BETA = 400.) ! dif from A=57 to B=59 is worth 0.08 seconds REAL TOT, F, FMELA, FMELB REAL DIF1, DIF2, DIF3 *----------------------------------------------------------------------- * Time separation distance *----------------------------------------------------------------------- DIF1 = MIN(ABS(TON(I) - TOFF(J)), ABS(TON(J) - TOFF(I))) * Overlap penalty DIF2 = DIM(MIN(TOFF(I), TOFF(J)), MAX(TON(I), TON(J))) * Pitch difference (use the Mel scale?) * Convert MIDI note number to Hertz * f = 440 @ 2^((n-69) / 12) F = 440. * 2. ** ((PITCH(I) - 69) / 12.) * Convert Hertz to Mel FMELA = 2595. * LOG10(1. + F / 700.) F = 440. * 2. ** ((PITCH(J) - 69) / 12.) FMELB = 2595. * LOG10(1. + F / 700.) DIF3 = ABS(FMELA - FMELB) * Make units comparable DIF3 = DIF3 / BETA * Combine TOT = DIF1 **2 + DIF2 **2 + DIF3 **2 EDIST = SQRT(TOT) RETURN END ! of EDIST *----------------------------------------------------------------------- * nextcomb - lists the next combination of n choose k in lexical order *----------------------------------------------------------------------- SUBROUTINE NEXTCOMB (N, K, R, DONE) IMPLICIT NONE INTEGER N, K, R ! I/O list LOGICAL DONE DIMENSION R(K) INTEGER C, I, J, P ! local variables * Begin IF (K .EQ. 0) THEN DONE = .TRUE. RETURN END IF I = K C = N - K * Check incremented value p at position i 10 P = R(I) + 1 IF (P > C + I) THEN ! carry I = I - 1 IF (I < 1) THEN DONE = .TRUE. RETURN END IF GO TO 10 END IF * Apply increment DO J = I, K ! roll over digits R(J) = P P = P + 1 END DO RETURN END ! of NEXTCOMB *----------------------------------------------------------------------- * puntal - assign stream ID's to a set of note events * Refer to: * Elaine Chew and Xiaodan Wu * "Separating Voices in Polyphonic Music: A Contig Mapping Approach" * Computer Music Modeling and Retrieval, pp. 1-20, 2005 * * see further references below * *___Name_______________Type__________I/O_______Description______________ * NOTE_TIME(N) Real In Note on times * PITCH(N) Integer In Note numbers * OFFON(N) Integer In NoteOn or NoteOff * N Integer In Number of note on/off events * V(N) Integer Out Voice membership of note events * MSV Integer Out Number of voices made * B(0:VLIM,N) Integer Neither Backpointers, reused for voice * NIV(N+1) Integer Neither Implicit voices each segment * NOV(N+1) Integer Neither Observed voices each segment * NOTIND(PLIM,N+1) Integer Neither What notes are in each sement * TOFF(N) Real Neither Note-off time each note-on * PLIST(N+1) Logical Neither Whether a partition is kept * IFAULT Integer Out Error code 0 = OK *----------------------------------------------------------------------- SUBROUTINE PUNTAL (NOTE_TIME, PITCH, OFFON, N, V, MSV, B, $ NIV, NOV, NOTIND, TOFF, PLIST, IFAULT) IMPLICIT NONE INTEGER PLIM, VLIM PARAMETER (PLIM = 24, ! polyphony requirement in General MIDI Level 1 $ VLIM = 10) ! fingers INTEGER PITCH, OFFON, N, V, MSV, B, NIV, NOV, NOTIND, IFAULT REAL NOTE_TIME, TOFF LOGICAL PLIST DIMENSION NOTE_TIME(N), PITCH(N), OFFON(N), V(N), B(0:VLIM,N), $ NIV(N+1), NOV(N+1), NOTIND(PLIM,N+1), TOFF(N), PLIST(N+1) * Local variables INTEGER $ BACK, ! a back pointer to best analysis $ CBEST, C, ! integer-valued connection scores $ DIF, $ H, I, J, K, L, ! counters $ INEW, IOLD, ! whether data is in even or odd column $ M, O, ! temp storage for niv(L), nov(L) $ N1, N2, ! note indices $ NSV, ! number of sounding voices $ NOTDEL ! index of a note event to delete INTEGER $ P, ! # partitions needing combining $ PALL, ! # of minimal partitions $ PSI(0:127), ! pointer to sound inception $ PVECT(VLIM), ! temporary storage $ R(VLIM), ! combination vector, lexical order $ RTRY(VLIM), ! trial combination $ REPTAB(VLIM), ! replacement table $ XI, XJ ! indices into note data REAL $ BIG, ! very large number $ COBSF, COBSM, ! penalty for mismatch of observed vs expected $ CTRANS, ! penalty for changing number of voices $ D(PLIM,PLIM), ! distances between simultaneous notes $ EDIST, ! external distance function $ RMACHINE, ! external function $ PIP, ! minimal time interval $ TNEW, ! time of current event $ TPREV, ! time of previous event $ U(0:VLIM,0:1), ! dynamic programming table $ UBEST, UTRY ! log-probabilities (arbitrary units) LOGICAL ALLDEC, ! all contigs have declining voice count $ DONE, $ USEDEC ! merge partitions with declining voice count SAVE COBSF, COBSM, CTRANS, PIP DATA COBSF, COBSM, CTRANS / -0.2, -1.0, -7.0 / DATA PIP / 0.035 / *----------------------------------------------------------------------- * Begin *----------------------------------------------------------------------- IFAULT = 0 BIG = RMACHINE(3) IF (N < 3) THEN ! trivial cases DO J = 1, N V(J) = 0 END DO MSV = 0 RETURN END IF * Limit the number of voices to make IF (VLIM > PLIM) THEN PRINT *, 'voice: Too many voices requested. Maximum ', PLIM IFAULT = 2 RETURN END IF *----------------------------------------------------------------------- * Segment *----------------------------------------------------------------------- DO K = 0, 127 PSI(K) = 0 END DO MSV = 0 NSV = 0 TPREV = 0 I = 0 DO J = 1, N TNEW = NOTE_TIME(J) IF (TNEW - TPREV > PIP) THEN I = I + 1 K = 0 DO L = 0, 127 IF (PSI(L) > 0) THEN K = K + 1 IF (K .LE. PLIM) NOTIND(K,I) = PSI(L) END IF END DO NOV(I) = MIN(NSV, PLIM) TPREV = TNEW END IF IF (OFFON(J) .EQ. 1) THEN IF (PSI(PITCH(J)) .EQ. 0) THEN PSI(PITCH(J)) = J NSV = NSV + 1 ELSE ! PRINT *, 'event ', J, ' attempt to turn on sounding note!' CONTINUE END IF ELSE IF (OFFON(J) .EQ. 0) THEN IF (PSI(PITCH(J)) > 0) THEN TOFF(PSI(PITCH(J))) = TNEW PSI(PITCH(J)) = 0 NSV = NSV - 1 ELSE ! PRINT *, 'event ', j, ' attempt to hush silent note!' CONTINUE END IF ELSE PRINT *, 'event ', J, ' bad offon value: ', OFFON(J) END IF END DO PALL = I * All notes briefer than a pip. Mark them stray. IF (PALL .EQ. 1) THEN DO I = 1, N V(I) = 0 END DO MSV = 0 RETURN END IF *----------------------------------------------------------------------- * Viterbi smoothing pass *----------------------------------------------------------------------- DO J = 0, VLIM ! uninformative distribution U(J,0) = 0. END DO * Forwards phase, for each time step DO I = 1, PALL * Alternate storage in array U(o,0:1) on even and odd j IF (BTEST(I, 0)) THEN INEW = 1 IOLD = 0 ELSE INEW = 0 IOLD = 1 END IF * For each # voices at present time step DO J = 0, VLIM UBEST = -BIG BACK = 0 ! avoid compiler warning * For each # voices at previous time step DO 10 K = 0, VLIM UTRY = U(K,IOLD) ! prior probability IF (UTRY < UBEST) GO TO 10 ! bail out * Penalize change IF (J .NE. K) UTRY = UTRY + CTRANS IF (UTRY > UBEST) THEN UBEST = UTRY BACK = K END IF 10 CONTINUE ! next k * Observation likelihood DIF = DIM(J, NOV(I)) ! too few observed notes UBEST = UBEST + COBSF * DIF DIF = DIM(NOV(I), J) ! too many observed notes U(J,INEW) = UBEST + COBSM * DIF * save pointer to most likely prior state B(J,I) = BACK END DO ! next j END DO ! next i * Find the highest scoring path UBEST = -BIG BACK = 0 ! avoid compiler warning IF (BTEST(PALL, 0)) THEN ! avoid compiler warning INEW = 1 ELSE INEW = 0 END IF DO J = 0, VLIM IF (U(J,INEW) > UBEST) THEN BACK = J UBEST = U(J,INEW) END IF END DO * Follow back pointers NIV(PALL) = 0 NOV(PALL) = 0 DO I = PALL, 2, -1 NIV(I) = MIN(BACK, NOV(I)) MSV = MAX(MSV, NIV(I)) BACK = B(BACK,I) END DO NIV(1) = MIN(BACK, NOV(1)) MSV = MAX(MSV, NIV(1)) ! End of Viterbi *----------------------------------------------------------------------- * Reduce polyphony *----------------------------------------------------------------------- DO L = 1, PALL M = NIV(L) O = NOV(L) IF (O > M) THEN ! remove some notes from this partition DO J = 1, O-1 XJ = NOTIND(J,L) IF (XJ < 1) PRINT*, 'notind=0 at j=', J, ' partit:', L DO I = J+1, O ! find all similarities XI = NOTIND(I,L) IF (XI < 1) PRINT*, 'notind=0 at i=', I, ' partit:', L D(I,J) = EDIST(NOTE_TIME, TOFF, PITCH, N, XI, XJ) END DO ! next i END DO ! next j * Try all combinations to find the maximum distance DO K = 1, M RTRY(K) = K END DO UBEST = 0. DONE = .FALSE. DO WHILE (.NOT. DONE) UTRY = 0. DO J = 1, M-1 DO I = J+1, M UTRY = UTRY + D(RTRY(I),RTRY(J)) END DO END DO IF (UTRY > UBEST) THEN DO K = 1, M R(K) = RTRY(K) END DO UBEST = UTRY END IF CALL NEXTCOMB(O, M, RTRY, DONE) END DO * Shorten note events not marked by combination vector R J = 1 DO 30 K = 1, O IF (J .LE. M) THEN IF (K .EQ. R(J)) THEN ! skip this note J = J + 1 GO TO 30 END IF END IF NOTDEL = NOTIND(K,L) DO I = L, PALL DONE = .TRUE. DO H = 1, NOV(I) IF (NOTIND(H,I) .EQ. NOTDEL) THEN ! swap to end NOTIND(H,I) = NOTIND(NOV(I),I) NOV(I) = NOV(I) - 1 ! forget existence NIV(I) = MIN(NIV(I), NOV(I)) ! keep consistent DONE = .FALSE. GO TO 20 END IF END DO ! next h 20 IF (DONE) GO TO 30 ! no trace of this event END DO ! next i 30 CONTINUE ! next k END IF ! too many events?! * Sort remaining events by frequency in ascending order DO K = 1, M PVECT(K) = PITCH(NOTIND(K,L)) END DO CALL QSORTI(PVECT, R, M) CALL IORDER(NOTIND(1,L), R, M) * Label local voice numbers DO K = 1, M B(K,L) = M - K + 1 END DO PLIST(L) = .TRUE. ! mark partition as not joined END DO ! next L PLIST(PALL) = .TRUE. ! DEBUG: DO I = 1, PALL IF (NIV(I) .NE. NOV(I)) THEN PRINT *, 'at ', I, 'niv nov mismatch: ', NIV(I), NOV(I) END IF END DO *----------------------------------------------------------------------- * Prioritized contig mapping * refer to: * Asako Ishigaki, Masaki Matsubara, Hiroaki Saito * "Prioritized Contig Combining to Segregate Voices in Polyphonic Music" * draft report, 2011 *----------------------------------------------------------------------- P = PALL - 2 ! remove all partitions except the beginning and end ALLDEC = .FALSE. 41 IF (P .EQ. 0) GO TO 45 IF (ALLDEC) THEN USEDEC = .TRUE. ELSE USEDEC = .FALSE. END IF ALLDEC = .TRUE. L = 1 42 IF (L .EQ. PALL-1) GO TO 41 ! end of list? IF (PLIST(L+1)) THEN I = L + 1 ! find voice count at last partit of contig DO WHILE (.NOT. PLIST(I+1)) I = I + 1 END DO IF (NOV(L) .LE. NOV(I)) THEN ! merge increasing ALLDEC = .FALSE. GO TO 44 ! and merge ELSE IF (USEDEC) THEN USEDEC = .FALSE. ALLDEC = .FALSE. GO TO 44 END IF END IF 43 L = L + 1 ! next partition GO TO 42 44 IF (NOV(L) > NOV(I)) THEN ! merge declining contigs H = I ! partition with fewer notes I = L ! partition with more notes ELSE ! merge increasing contigs H = L END IF M = NOV(H) O = NOV(I) DO K = 1, M RTRY(K) = K END DO CBEST = 32767 DONE = .FALSE. DO WHILE (.NOT. DONE) ! find connection costs for all combinations C = 0 DO J = 1, M N1 = NOTIND(J,H) N2 = NOTIND(RTRY(J),I) IF (N1 .EQ. N2) THEN C = C - 1 ! subtract for same note ELSE DIF = ABS(PITCH(N1) - PITCH(N2)) C = C + DIF ! penalty is pitch difference END IF END DO ! next j IF (C < CBEST) THEN ! lowest connection cost DO K = 1, M R(K) = RTRY(K) END DO CBEST = C END IF CALL NEXTCOMB(O, M, RTRY, DONE) END DO * Replace local note numbers with global note numbers DO J = 1, M REPTAB(B(J,H)) = B(R(J),I) ! fill replacement table END DO * Renumber voices DO K = H, 1, -1 DO J = 1, NOV(K) B(J,K) = REPTAB(B(J,K)) END DO IF (PLIST(K)) GO TO 50 END DO 50 PLIST(L+1) = .FALSE. P = P - 1 GO TO 43 * Final results 45 DO J = 1, N V(J) = 0 END DO DO J = 1, PALL DO I = 1, NOV(J) K = B(I,J) L = NOTIND(I,J) IF (V(L) .EQ. 0) V(L) = K END DO END DO RETURN END ! of puntal *+++++++++++++++++++++ End of file puntal.f ++++++++++++++++++++++++++++