*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% * rhythm.f - subroutines for music analysis * by Andy Allinger, 2011-2017, released to the public domain * This program may be used by any person for any purpose. *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% * OSCILL - a frequency analyzer for pulse trains * A bank of oscillatory filters responds to a sequence of events in time. * Each filter has the undriven form of a decaying logarithmic spiral. * Filters respond to events by adjusting their amplitude, phase, and frequency. * For general reference, see: * * Theodosius Pavlidis * Biological Oscillators: their Mathematical Analysis * Academic Press, 1973 * * Frank C. Hoppensteadt * "Signal Processing by Model Neural Networks" * SIAM Review v.34 n.3 pp.426-444 * Society for Industrial and Applied Mathematics, September 1992 * * Edward W. Large and John F. Kolen * "Resonance and the Perception of Musical Meter" * Connection Science v.6 n.1 pp.177-208, 1994 * *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% * I/O LIST *___Name_________Type______In/Out___Description_________________________ * T(M) Real In Onset times of events * S(M) Real In Saliences of events * M Integer In Number of events * N Integer In Number of evaluation timepoints * O Integer In Number of oscillators * STEP Real In Step size in time * FDEEP Real In Center frequency of deepest oscillator * FSPAC Real In Frequency ratio of adjacent oscillators * FRANG Real In Ratio of min to max frequency of a single * oscillator * A(0:O,0:N) Real Out Amplitude * P(0:O,0:N) Real Out Phase * F(0:O,0:N) Real Out Frequency * IFAULT Integer Out nonzero on error * *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SUBROUTINE OSCILL (T,S,M,N,O,STEP,FDEEP,FSPAC,FRANG,A,P,F,IFAULT) IMPLICIT NONE INTEGER M, N, O, IFAULT REAL T, S, STEP, FDEEP, FSPAC, FRANG, A, P, F DIMENSION T(M), S(M), A(0:O,0:N), P(0:O,0:N), F(0:O,0:N) * constants INTEGER Q REAL ETA, L, PI, RINIT, RMIN, RMAX, SM, TWOPI, V, ZETA PARAMETER ( $ ETA = 1., ! target salience per event * $ g = 7., ! sharpening coefficient $ L = 3., ! half-life, seconds $ PI = 3.141592654, ! $ Q = 24, ! number of initial phase offsets $ RINIT = 0.1, ! initial amplitude $ RMIN = 0.0001, ! distance to stay away from origin $ RMAX = 0.9999, ! closest approach to unit circle $ SM = 0.05, ! minimal salience $ TWOPI = 2 * PI, ! $ V = 1.0, ! volume of the beaker $ ZETA = 2.618034 ) ! target events per cycle * local variables INTEGER H, I, J, K ! counters REAL $ BETA, ! decay rate $ B1, ! decay for whole step $ BEST, ! best amplitude result $ BIG, $ DS, ! adjustment to each salience $ DT, DX, DY, ! changes in t, x and y $ FMIN, F0, FMAX, ! min, central, max freq. each oscillator $ FM, ! M as Real $ GAMMA, MU, ! deviation, average of salience $ NU, ! frequency, in Hertz $ PHI ! phase, measured in beats elapsed REAL $ RHO, R2, ! amplitude, amplitude squared $ S1, ! adjusted salience $ SMAX, ! greatest entry in S $ THETA, ! phase % 1, in radians $ T0, ! time of prior event $ TNEXT, ! time of next timepoint $ TPREV, ! time last processed event or step $ TT, ! same as T(i) $ U, ! sharpening coefficient given rho, theta $ W, XI, ! weight factor of salience, crit. fraction $ X, Y ! x = rho @ cos(theta) ; y = rho @ sin(theta) DOUBLE PRECISION $ DEPS, ! machine double precision $ C2, ST0, ST1, ST2, SP, STP, DET, ! least-squares sums $ SS, SS2 ! salience sum, sum of squares LOGICAL E ! whether exponentiation is required !!! stats INTEGER processed, ignored * external functions REAL RMACHINE, $ SNAP, SLOG, SSIN, SCOS, SARC ! swift EXP, LOG, SIN, COS, ATAN2 DOUBLE PRECISION DMACHINE *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% * begin IFAULT = 0 ! clear error codes BIG = RMACHINE(3) DEPS = DMACHINE(1) * validate IF (M < 1 .OR. N < 1 .OR. O < 1) THEN IFAULT = 1 RETURN END IF * check that enough timepoints have been allocated for the given events IF (T(M) > N * STEP) THEN PRINT *, 'oscill: t(m) is ', T(M), ' n@step is ', N * STEP IFAULT = 2 RETURN END IF * get f0, take root of frang F0 = FDEEP / FSPAC ! OK to multiply by fspac on iteration k = 1 FMIN = F0 / SQRT(FRANG) FMAX = F0 * SQRT(FRANG) BETA = LOG(0.5) / L B1 = EXP(BETA * STEP) * Total salience SS = 0.D0 SS2 = 0.D0 SMAX = 0. DO I = 1, M SS = SS + S(I) SS2 = SS2 + S(I) **2 IF (S(I) > SMAX) SMAX = S(I) END DO FM = FLOAT(M) MU = SNGL(SS / FM) GAMMA = 0.6745 * SNGL(SQRT(SS2 / FM - (SS / FM) **2)) *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% * for each oscillator *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DO K = 1, O F0 = F0 * FSPAC ! central frequency FMIN = FMIN * FSPAC ! minimum frequency FMAX = FMAX * FSPAC ! maximum frequency BEST = -BIG XI = (FM - ZETA * T(M) * F0) / FM ! fraction of events to ignore XI = MIN(MAX(RMIN, XI), RMAX) DS = MU + GAMMA * TAN(PI * (XI - 0.5)) ! Cauchy cumulative dist. DS = MIN(MAX(0., DS), SMAX - SM) SS = 0.D0 J = 0 ! count nonzero events DO I = 1, M IF (S(I) > DS) THEN SS = SS + (S(I) - DS) J = J + 1 END IF END DO W = ETA * J / MAX(SM, SNGL(SS)) ignored = 0 processed = 0 DO H = 0, Q-1 ! for each initial phase E = .FALSE. ! no stimulus event intervening RHO = RINIT ! initial amplitude PHI = FLOAT(H) / FLOAT(Q) ! initial phase NU = F0 ! default frequency A(0,0) = RHO P(0,0) = PHI F(0,0) = NU I = 1 T0 = 0. TNEXT = 0. TPREV = 0. ST0 = 0.D0 ST1 = 0.D0 ST2 = 0.D0 SP = 0.D0 STP = 0.D0 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% * for each evaluation point *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DO J = 1, N TNEXT = TNEXT + STEP 10 IF (I > M) GO TO 50 ! no more events ? TT = T(I) IF (TT > TNEXT) GO TO 50 ! no event before tnext ? *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% * apply the stimulus - first advance r, phi to present * then make an adjustment in Cartesian coordinates *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% * phi' = phi + nu @ dt DT = TT - TPREV PHI = PHI + NU * DT * rho' = rho @ e^(b @ dt) RHO = RHO * SNAP(BETA * DT) ! approximate exponential decay RHO = MIN(MAX(RMIN, RHO), RMAX) *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% * Ignore stimuli below critical value. *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% S1 = S(I) - DS IF (S1 .LE. 0.) THEN ignored = ignored + 1 I = I + 1 TPREV = TT E = .TRUE. GO TO 10 ELSE processed = processed + 1 S1 = S1 * W END IF THETA = TWOPI * MOD(PHI, 1.0) ! 0 .LE. theta < 2 @ pi X = RHO * SCOS(THETA) ! cosine Y = RHO * SSIN(THETA) ! sine *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% * Use activation sharpening. Refer to: * J. Devin McAuley * Perception of Time as Phase: Toward an Adaptive Oscillator Model of * Rhythmic Pattern Processing * University of Indiana Ph.D. Thesis, July 1995, p. 58 * u = ((1 + COS(theta)) / 2) ^ (g @ r) * * also refer to: * J. Devin McAuley * Learning to Perceive and Produce Rhythmic Patterns * in an Artificial Neural Network * Indiana University Department of Computer Science * Technical Report 371, 1 February 1993 * which cites: * Carme Torras * Temporal-Pattern Learning in Neural Models * Berlin: Springer-Verlag, 1985 * * Another choice might be the wrapped Cauchy distribution. Refer to: * Claudio Agostinelli * R CRAN package circular * 22 May 2006 * u = (1 - r^2) / ( 2 @ pi @ (1 + r^2 - 2 @ r @ COS(theta)) ) *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% * u = ((1 + COS(theta)) / 2) ** (g * rho) R2 = RHO **2 U = (1. - R2) / (TWOPI * (1. + R2 - 2. * X)) *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% * Consider a chemical oscillator with reactant concentrations x and y, * subject to a stimulus of adding a substance with y = 0, and x = v, * the maximum concentration of x. Then the phase transition will fulfill: * dy/dx = -y / (v - x) * Refer to: * Arthur T. Winfree * The Geometry of Biological Time * Springer-Verlag, 1980, p. 139 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% U = (U * S1) / SQRT((V-X) **2 + Y **2) DX = (V-X) * U DY = (-Y) * U X = MIN(X + DX, V) IF (ABS(DY) .GE. ABS(Y)) THEN ! don't overcorrect IF (Y < 0.) PHI = PHI + 1. ! increment phase Y = 0. ELSE Y = Y + DY END IF RHO = SQRT(X **2 + Y **2) IF (RHO < RMIN) THEN ! don't fall on zero (the phase singularity) X = RMIN RHO = RMIN END IF IF (RHO > RMAX) RHO = RMAX THETA = SARC(Y, X) ! arctangent PHI = AINT(PHI) + THETA / TWOPI *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% * Measure frequency as the slope of a line of phase as a function of time * nu = dphi / dt *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% * decay the running sums for the least squares line C2 = DBLE(SNAP(BETA * (TT-T0))) ! decay factor ST0 = ST0 * C2 + 1.D0 ! SUM t^0 ST1 = ST1 * C2 + DBLE(TT) ! SUM t^1 ST2 = ST2 * C2 + DPROD(TT, TT) ! SUM t^2 SP = SP * C2 + DBLE(PHI) ! SUM phi STP = STP * C2 + DPROD(TT, PHI) ! SUM t @ phi DET = ST0 * ST2 - ST1 * ST1 ! determinant IF (DET > DEPS) NU = SNGL( (ST0 * STP - SP * ST1) / DET ) ! NU = MIN(MAX(FMIN, NU), FMAX) ! impose limits IF (NU < FMIN .OR. NU > FMAX) THEN ! reset NU = F0 RHO = RMIN ST0 = 0. ST1 = 0. ST2 = 0. SP = 0. STP = 0. END IF I = I + 1 T0 = TT TPREV = TT E = .TRUE. GO TO 10 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 50 CONTINUE ! advance time to an evaluation point for A, P, F *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% IF (E) THEN DT = TNEXT - TPREV RHO = RHO * SNAP(BETA * DT) ! approximate exponential function PHI = PHI + NU * DT ELSE RHO = RHO * B1 ! avoid an exponentiation when possible PHI = PHI + NU * STEP END IF IF (RHO < RMIN) RHO = RMIN A(0,J) = RHO P(0,J) = PHI F(0,J) = NU TPREV = TNEXT E = .FALSE. END DO ! next j * find sum of logarithm of amplitude U = 0. DO J = 1, N U = U + SLOG(A(0,J)) END DO IF (U > BEST) THEN ! new best answer BEST = U DO J = 0, N A(K,J) = A(0,J) P(K,J) = P(0,J) F(K,J) = F(0,J) END DO END IF END DO ! next h ! PRINT *, 'osc# ', k, ' processed ', processed, ' events ', ! $ 'and ignored ', ignored, ! $ ' using ds: ', DS, ! $ ' using w: ', W END DO ! next k RETURN END ! of OSCILL *----------------------------------------------------------------------- * orlp - logical OR of logarithms of probabilities *----------------------------------------------------------------------- SUBROUTINE ORLP (X, N, RESULT) IMPLICIT NONE INTEGER N REAL X, RESULT DIMENSION X(N) INTEGER J REAL XMAX *----------------------------------------------------------------------- * Begin. *----------------------------------------------------------------------- IF (N < 1) RETURN XMAX = X(1) DO J = 2, N IF (X(J) > XMAX) XMAX = X(J) END DO RESULT = 0. DO J = 1, N RESULT = RESULT + EXP(X(J) - XMAX) END DO RESULT = LOG(RESULT) + XMAX RETURN END ! of ORLP *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% * strata - The best rhythic patterns are chosen from the output of the * oscillators with the Viterbi algorithm. Dropped beats are filled in. * The phase and time estimates are approximated with a smoothing a spline * yielding a function for phase given time, for each of several rhythmic * strata. Refer to : * * Anssi P. Klapuri, Antti J. Eronen, and Jaakko T. Astola, * "Analysis of the Meter of Acoustic Musical Signals" * IEEE Transactions on Speech and Audio Processing, 2004 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% * TRANSFER LIST *__Name___________Type______In/Out___Description_________________________ * T(m) Real in times of events * S(m) Real in saliences of events * L Integer in number of levels ("strata") to make * TL Integer in index of the tactus level * m Integer in number of events * n Integer in number of evaluation timepoints * NODES Integer in number of spline interpolation points * step Real in duration between oscillator evaluations * urpar(8) Real in parameters for probability functions * opt_S Logical in option to calculate an average tempo * opt_T Logical in option to get tempo segments * full Logical in calculate all the strata * prob Real out total log probability * CX(NODES+1) Real out spline coefficients * CY(NODES,L) Real out " " * CA(NODES,L) Real out " " * CB(NODES,L) Real out " " * CC(NODES,L) Real out " " * TX(0:n) Real out tempo line segment coefficients * TA(0:n) Real out " " * nt Integer out number of tempo line segments * tempo Real out average tempo * A(0:o,0:n,2) Real neither oscillator amplitude * P(0:o,0:n,2) Real neither oscillator phase * F(0:o,0:n,2) Real neither oscillator output frequency * B(o,n) Real neither back pointers * C(0:n,L) Real neither best paths thru the Viterbi lattice * W(0:n) Real neither weights for least-squares spline * X(0:n) Real neither time input for spline * Y(0:n) Real neither phase input for spline * Z(7*NODES+25) Real neither work array for spline * fault Integer out error code *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SUBROUTINE strata (T, S, L, TL, M, N, NODES, step, urpar, opt_S, $ opt_T, usecfg, full, prob, CX, CY, CA, CB, CC, TX, TA, nt, $ tempo, A, P, F, B, C, W, X, Y, Z, fault) IMPLICIT NONE INTEGER o PARAMETER (o = 27) ! number of oscillators INTEGER L, TL, m, n, NODES, nt, B, C, fault REAL T, S, step, urpar, prob, CX, CY, CA, CB, CC, TX, TA, tempo, $ A, P, F, W, X, Y, Z LOGICAL opt_S, opt_T, usecfg, full DIMENSION T(m), S(m), urpar(8), CX(NODES+1), CY(NODES,L), $ CA(NODES,L), CB(NODES,L), CC(NODES,L), TX(0:n), TA(0:n), $ A(0:o,0:n,2), P(0:o,0:n,2), F(0:o,0:n,2), B(o,n), C(0:n,L), $ W(0:n), X(0:n), Y(0:n), Z(7*NODES+25) * local variables INTEGER GN PARAMETER (GN = 1024) ! evaluation points of spline INTEGER $ back, ! back pointer $ h, i, ii, j, k, k1, ! counters $ jnew, jold ! even or odd time step REAL $ afact, ! scale down amplitude importance $ BIG, ! maximum number $ CP(4), ! cubic polynomial coefficients $ dp, dt, ! change of phase, change of time $ g, ! the Large-Klapuri g function $ gfact, ! scale down resonance importance $ GX(gn), GY(gn), GA(gn), GB(gn), GC(gn), ! g function parameters $ mu, nu, ! average frequency, present frequency $ pobs, ppre, ! observed, predicted phase $ r, ! ratio of frequencies $ rpar(8), ! rhythm parameters $ sigi, ! parameters for frequency change, initial ! frequency distribution, and phase error $ tcum, ! cumulative time $ twovf, twovi, twovp, ! double variance $ U(o,0:1), ! cumulative log probabilities at even/odd steps $ ubest, utry, ! log-prob. of likliest, other transitions $ V(o) ! log-prob. of each prior state LOGICAL first, $ flip(o) SAVE first, GX, GY, GA, GB, GC, rpar * Constants * AMPIMP parameter will be scaled by STEP. Thus, it is per * unit time. Figure 4 eighth notes per second, at 5% phase error * gives a penalty of -0.5. A resonator tracking at 10% amplitude * has a log of about -2.3 INTEGER minstp ! minimum timesteps for SPFIT REAL ampimp, gfnimp, ! importance of amplitude, resonance $ fdeep, frang, fspac ! deepest natural frequency, range, spacing PARAMETER (ampimp = 0.2, ! because amplitude is not probability $ gfnimp = 0.3, ! not a real probability either $ fdeep = 0.1, ! center frequency of slowest oscillator $ frang = 4. / 3., ! range of each oscillator $ fspac = 38. / 31., ! chosen to avoid spurious resonances $ minstp = 7) ! don't break spline interpolator * External functions INTEGER CEIL ! round up REAL $ rmachine, ! machine constants, single precision $ SLOG ! Swift logarithm !! timing info REAL TARRAY(2), RESULT * The default rhythm parameters are from Klapuri, Eronen, & Astola DATA rpar / 0.18, ! tatum period ! $ 0.39, ! standard deviation of log ratio $ 0.27, ! standard deviation of log ratio $ 0.55, ! tactus period $ 0.28, ! standard deviation of log ratio $ 2.10, ! measure period $ 0.26, ! standard deviation of log ratio ! $ 0.2, ! std deviation of tempo changes $ 0.1, ! std deviation of tempo changes $ 0.1 / ! std deviation of phase error DATA first / .TRUE. / ! calculate the g function on first run *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% * begin *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% * define constants BIG = rmachine(3) afact = ampimp * step gfact = gfnimp IF (first) THEN CALL gfunc(GX, GY, GA, GB, GC, gn) first = .FALSE. END IF jnew = 0 ! avoid meaningless compiler warnings * custom config data IF (usecfg) THEN DO i = 1, 8 rpar(i) = urpar(i) END DO END IF * compare number of strata against number of oscillators IF (L > o) THEN PRINT *, 'Too many rhythmic strata for oscillators allocated' fault = 1 RETURN END IF CALL DTIME (tarray, result) * set up evaluation points dt = (n + 0.5) * step / NODES tcum = 0. DO j = 1, NODES+1 CX(j) = tcum tcum = tcum + dt END DO * run the data thru the filterbank CALL OSCILL (T,S,m,n,o,step,fdeep,fspac,frang,A,P,F,fault) IF (fault .NE. 0) PRINT *, 'oscill returns error #', fault !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! PRINT *, 'forward oscillators: ' ! 41 FORMAT ('step: ', I5) ! 42 FORMAT (F9.3, $) ! 43 FORMAT () ! DO j = 1, n ! print 41, j ! DO h = 1, o ! PRINT 42, A(h,j,1) ! END DO ! PRINT 43 ! DO h = 1, o ! print 42, P(h,j,1) ! END DO ! PRINT 43 ! DO h = 1, o ! PRINT 42, F(h,j,1) ! END DO ! PRINT 43 ! END DO !!!!!!!!!!!!!!!!!!!! * Reverse order of the data tcum = n * step ! PRINT *, 'using tcum of ', tcum DO i = 1, m ! negate time T(i) = -T(i) + tcum T(i) = MIN(MAX(0., T(i)), tcum - .001) END DO DO i = 1, m/2 ! swap order r = T(i) T(i) = T(m-i+1) T(m-i+1) = r r = S(i) S(i) = S(m-i+1) S(m-i+1) = r END DO * Filter the reversed sequence. CALL OSCILL (T, S, m, n, o, step, fdeep, fspac, frang, $ A(0,0,2), P(0,0,2), F(0,0,2), fault) IF (fault .NE. 0) THEN PRINT *, 'rev. ocsill returns #', fault PRINT *, 'T(m): ', T(m) PRINT *, 'n*step: ', N * STEP PRINT *, 'tcum: ', tcum END IF * Reflect phase DO h = 1, o dp = FLOAT(CEIL(P(h,n,2))) DO j = 0, n P(h,j,2) = -P(h,j,2) + dp END DO END DO * Keep backwards result until forwards pass has higher amplitude DO h = 1, o flip(h) = .TRUE. END DO DO j = 0, n DO h = 1, o IF (flip(h)) THEN IF (A(h,j,1) < A(h,n-j,2)) THEN A(h,j,1) = A(h,n-j,2) P(h,j,1) = P(h,n-j,2) F(h,j,1) = F(h,n-j,2) ELSE flip(h) = .FALSE. END IF END IF END DO END DO CALL DTIME (tarray, result) ! PRINT *, 'rhythm: time filtering: ', RESULT prob = 0. twovf = 2 * rpar(7) **2 ! same for all strata twovp = 2 * rpar(8) **2 DO h = 1, L ! for each stratum CALL DTIME (tarray, result) *||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| * Viterbi algorithm: * The basic equation for conditional probability is: * P(A&B) = P(A|B) @ P(B) * Suppose each event e in a sequence of n events is conditional only on * the event before it. This yields the equation for a Markov chain: * P(E) = P(e[1]) @ PRODUCT P(e[i+1]|e[i]) where i = 1...n-1 * where P(E) is the probability of the entire sequence. * Suppose also that knowledge of each event is supplemented by an observation * at each step. Applying Bayes's rule: * P(E&O) = P(E|O) @ P(O) = P(O|E) @ P(E) * disregarding P(O) which is the same for all events, * P(E|O) = P(O|E) @ P(E) * Bayesian reasoning is less elegant when the quantities involved must * be measured. For the present problem, transition probability P(e[i+1]|e[i]) * is based on the assumption that frequency remains steady and therefore phase * advances steadily. The observation likelihood, P(O|E), is taken to be * proportional to the oscillator amplitude, and to be higher if a frequency * is proportional to the tactus (first detected, (strongest?)) frequency. * The Viterbi algorithm asks at each step in the sequence, for each possible * event, what event was most likely at the prior step? The answers for each * step are kept, and at the end of the sequence the most likely chain of * events may be traced back. * *||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| * choose parameters for initial frequency distribution IF (h .EQ. 1) THEN ! tatum mu = rpar(1) sigi = rpar(2) ELSE IF (h .EQ. TL) THEN ! tactus mu = rpar(3) sigi = rpar(4) ELSE IF (h .EQ. TL+1) THEN ! measure mu = rpar(5) sigi = rpar(6) ELSE IF (h < TL) THEN ! subtatum? mu = 0.09 sigi = 0.40 ELSE ! strophe? mu = 4.2 sigi = 0.40 END IF twovi = 2. * sigi **2 IF (h > 1) gfact = gfnimp / FLOAT(h - 1) *||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| * Log-normal initial frequency distribution: * * p(T) = (1 / T sigma (2@pi)^(1/2)) @ e^((-1/2@sigma^2) @ (ln(T/m))^2) * * log probability, disregarding constants is then: * * ln p(T) = (-1/2@sigma^2) @ ln(T/m)^2 - ln(T) *||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| DO k = 1, o nu = F(k,0,1) r = mu * nu ! frequency @ time: non-dimensional IF (r > 1.) r = 1. / r r = SLOG(r) r = -(r **2) / twovi IF (nu > 0.) r = r + LOG(nu) U(k,0) = r END DO * forwards phase, for each time step DO 30 j = 1, n IF (BTEST(j, 0)) THEN ! alternate storage in array U jnew = 1 jold = 0 ELSE jnew = 0 jold = 1 END IF * for each oscillator at present time step DO 20 k = 1, o IF (A(k,j,1) < 0.) THEN ! this cell has been marked U(k,jnew) = -BIG GO TO 20 END IF ubest = -BIG nu = F(k,j,1) ppre = P(k,j,1) - nu * step ! predicted value for prior phase *||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| * transition probabilities *||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| * for each oscillator at previous time step ii = 0 DO 10 i = 1, o IF (A(i,j-1,1) < 0.) GO TO 10 ! this cell has been marked ii = ii + 1 utry = U(i,jold) ! prior probability * penalize change in frequency r = nu / F(i,j-1,1) IF (r > 1.) r = 1. / r ! restrict to 0 < r .LE. 1 r = SLOG(r) utry = utry - (r **2) / twovf * penalize phase mismatch dp = ppre - P(i,j-1,1) r = dp - ANINT(dp) utry = utry - (r **2) / twovp V(ii) = utry IF (utry > ubest) THEN ubest = utry back = i END IF 10 CONTINUE ! next i !! DEBUG: IF (ubest .LE. -BIG) THEN PRINT *, 'no way out from j=', j, ' k=', k, 'ii=', II PRINT *, 'unew=', U(k,jnew) END IF * save pointer to most likely prior state B(k,j) = back CALL ORLP (V, ii, utry) ! add probabilities *||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| * observation likelihood *||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| * oscillator amplitude utry = utry + afact * SLOG(A(k,j,1)) * is frequency a Farey relation to previous strata? DO i = 1, h-1 r = F(C(j,i),j,1) / nu IF (r > 1.) r = 1. / r ! evaluate g function CALL SCOMP(GX, GY, GA, GB, GC, gn, r, g, 1, fault) IF (fault .NE. 0) PRINT *, 'scomp: error#', fault utry = utry + gfact * g END DO * does frequency match prior expected frequency? r = mu * nu IF (r > 1.) r = 1. / r r = SLOG(r) r = -(r **2) / twovi IF (nu > 0.) r = r + LOG(nu) utry = utry + r U(k,jnew) = utry 20 CONTINUE ! next k 30 CONTINUE ! next j * find the highest scoring path ubest = -BIG DO k = 1, o IF (U(k,jnew) > ubest) THEN back = k ubest = U(k,jnew) END IF END DO * cumulative probability prob = prob + ubest * trace back pointers DO j = n, 1, -1 IF (back < 1 .OR. back > o) THEN PRINT *, 'bad back pointer! :', back, ' j=', j, ' n=', n PRINT *, 'B= ', B END IF C(j,h) = back back = B(back,j) END DO C(0,h) = back ! End of Viterbi !!!!!!!!!!!!!! debugging dump ! PRINT *, 'stratum ', h, ' C[]=' ! 45 FORMAT (I3, $) ! DO j = 0, n ! PRINT 45, C(j,h) ! END DO !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! CALL DTIME (tarray, result) ! PRINT *, 'rhythm: level ', h, ' time in DP lattice: ', result *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% * replace dropped beats *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% W(0) = A(back,0,1) X(0) = 0. Y(0) = P(back,0,1) tcum = 0. k1 = back DO j = 1, n k = k1 ! index of previous k1 = C(j,h) ! index of new ppre = P(k,j-1,1) + step * F(k1,j,1) ! backwards integral pobs = P(k1,j,1) pobs = pobs + ANINT(ppre - pobs) ! adjust by an integer P(k1,j,1) = pobs W(j) = A(k1,j,1) A(k1,j,1) = -1. ! prevent this cell from being used again tcum = tcum + step X(j) = tcum Y(j) = pobs END DO ! next j *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% * fit a least-squares spline *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% IF (n < minstp) THEN ! fit a cubic polynomial CALL cfit(W, X, Y, n+1, CP, fault) IF (fault .NE. 0) PRINT *, 'cfit: error#', fault CY(1,h) = CP(1) CA(1,h) = CP(2) CB(1,h) = CP(3) CC(1,h) = CP(4) * restore X X(0) = 0. DO J = 1, N X(J) = X(J-1) + step END DO ELSE * least squares spline fitting by deBoor and Morris (NSWC library) CALL SPFIT (X, Y, W, n+1, CX, NODES+1, $ CY(1,h), CA(1,h), CB(1,h), CC(1,h), Z, fault) IF (fault .NE. 0) PRINT *, 'spfit: error #', fault END IF * tactus level only IF (h .EQ. TL) THEN * average tempo (for dance mode) IF (opt_S) THEN CALL wlsl (W, X, Y, n+1, r, tempo, fault) ! re-use r IF (fault .NE. 0) PRINT *, 'wlsl: error#', fault END IF * line-segment approximation of tempo curve IF (opt_T) THEN CALL mullin (W, X, Y, TX, TA, n+1, nt, fault) IF (fault .NE. 0) PRINT *, 'mullin: error #', fault END IF ! opt_T ? CALL DTIME (tarray, result) ! PRINT *, 'rhythm: level ', h, ' time with splines: ', result IF (.NOT. full) RETURN ! don't need all the splines END IF ! tactus level ? END DO ! next h RETURN END ! of strata *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% * add data to rhythm parameters * *__Name________Type______________In/Out___Description___________________ * T(M) Real In event times * S(M) Real In saliences of events * M Integer In # of events * N Integer In # evaluation points * STEP Real In time between evaluation points * A(0:o,0:n) Real Neither amplitude for whole filterbank * P(0:o,0:n) Real Neither phase for whole filterbank * F(0:o,0:n) Real Neither frequency for whole filterbank * rsum(9) Double Precision Both sums & sums of squares (see next f'n) *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SUBROUTINE adrp (T, S, m, n, step, A, P, F, rsum) IMPLICIT NONE INTEGER o PARAMETER (o = 3) ! how many oscillators INTEGER M, N REAL T, S, STEP, A, P, F DOUBLE PRECISION rsum DIMENSION T(m), S(m), A(0:o,0:n), P(0:o,0:n), F(0:o,0:n), rsum(9) * local variables INTEGER fault, $ i, j ! counters REAL $ dp, ! change in phase $ e, ! remainder phase error $ fdeep, fspac, frang, ! input to OSCILL $ flog, ! logarithm of ratio of consecutive frequencies $ lnt, ! logarithm of period (to match Klapuri) $ nu ! frequency SAVE fdeep, fspac, frang DATA fdeep, fspac, frang / 0.5, 4., 5. / * begin fault = 0 * one oscillator for tatum, tactus, measure CALL oscill (T,S,m,n,o,step,fdeep,fspac,frang,A,P,F,fault) IF (fault .NE. 0) PRINT *, 'adrp: error in oscillators' DO i = 1, 3 ! for each stratum DO j = 1, n nu = F(i,j) lnt = -LOG(nu) IF (i .EQ. 1) THEN ! measure rsum(5) = rsum(5) + lnt rsum(6) = rsum(6) + DPROD(lnt, lnt) ELSE IF (i .EQ. 2) THEN ! tactus rsum(3) = rsum(3) + lnt rsum(4) = rsum(4) + DPROD(lnt, lnt) ELSE ! tatum rsum(1) = rsum(1) + lnt rsum(2) = rsum(2) + DPROD(lnt, lnt) END IF * frequency change flog = LOG( nu / F(i,j-1) ) rsum(7) = rsum(7) + DPROD(flog, flog) * phase mismatch dp = P(i,j) - nu * step ! what prior phase should be dp = P(i,j-1) - dp ! difference e = dp - ANINT(dp) ! remainder phase error rsum(8) = rsum(8) + DPROD(e, e) END DO ! next j END DO ! next i * count total events rsum(9) = rsum(9) + n RETURN END ! of adrp *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% * rsum DOUBLE PRECISION rsum(1) sum of tatum log period * rsum(2) sum of square tatum log period * rsum(3) sum of tactus log period * rsum(4) sum of square tactus log period * rsum(5) sum of measure log period * rsum(6) sum of square measure log period * rsum(7) sum of square log frequency ratio * rsum(8) sum of square phase error * rsum(9) total evaluation points * * rpar REAL the estimated parameters: * rpar(1) mean tatum period * rpar(2) tatum standard deviation * rpar(3) tactus mean period * rpar(4) tactus standard deviation * rpar(5) measure mean period * rpar(6) measure standard deviation * rpar(7) std dev. of logarithm of frequency * ratio from one step to the next * This probably depends on the * value chosen for step * rpar(8) std. dev. of remainder phase error * ie- how far observed phase differs * from what would be predicted based * on a steady tempo *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SUBROUTINE MKRPAR (rsum, rpar) IMPLICIT NONE REAL rpar DOUBLE PRECISION rsum DIMENSION rsum(9), rpar(8) DOUBLE PRECISION ave, dev, n n = rsum(9) ave = rsum(1) / n dev = SQRT(rsum(2) / n - ave **2) rpar(1) = SNGL(EXP(ave)) ! units of time rpar(2) = SNGL(dev) ave = rsum(3) / n dev = SQRT(rsum(4) / n - ave **2) rpar(3) = SNGL(EXP(ave)) rpar(4) = SNGL(dev) ave = rsum(5) / n dev = SQRT(rsum(6) / n - ave **2) rpar(5) = SNGL(EXP(ave)) rpar(6) = SNGL(dev) * multiply by 3 for the three levels counted n = n * 3.D0 dev = SQRT(rsum(7) / n) ! assume mean of zero rpar(7) = SNGL(dev) dev = SQRT(rsum(8) / n) rpar(8) = SNGL(dev) RETURN END ! of MKRPAR ************************************************************************ * interpolate the g function ************************************************************************ SUBROUTINE gfunc (X, Y, A, B, C, n) IMPLICIT NONE INTEGER i, j, n, ibeg, iend, ifault REAL X(n), Y(n), A(n), B(n), C(n), sr2pi REAL farey(23), width(23), height(23), wf, hf, tot REAL alpha, beta ************************************************************************ * These constants represent the strength at which an oscillator latches on * to a signal, depending on the ratio of the signal frequency to the * oscillator frequency. The array farey is made up of fractions from the * Stern-Brocott-Farey sequence. Width and height refer to an "empirical regime * diagram" by Large, in which width represents a tolerance around an ideal * integer ratio, and height is a parameter that controls how strongly * the oscillator entrains to its input. These variables can be reinterpreted * as intensity (height), and standard deviation (width) around a mean (ratio). * Refer to: * Edward Wilson Large * Dynamic Representation of Musical Structure * Ph.D. Thesis, Ohio State University, 1994, Fig. 39 ************************************************************************ DATA farey / 0.0, .125, .1428571429, 0.1666666667, 0.2, 0.25, $ 0.2857142858, 0.3333333333, 0.375, 0.400, 0.4285714286, 0.5, $ 0.5714285714, 0.6, 0.625, 0.6666666667, 0.7142857143, 0.75, $ 0.8, 0.8333333333, 0.8571428571, 0.875, 1.000 / DATA width / 18, 2, 3, 4, 7, 10, 3, 12, 3, 7, 3, 19, 4, 7, 2, 13, $ 3, 10, 7, 5, 3, 2, 21 / DATA height / 5, 2, 2, 2.2, 2.2, 3.3, 1.9, 3.3, 1.7, 2, 1.7, 4, $ 1.6, 1.8, 1.6, 3.5, 1.4, 3, 1.6, 1.5, 1.2, 1.2, 5 / * scale factors DATA wf, hf / 2500., 3.16 / DATA alpha, beta / 0., 0. / ! not used ************************************************************************ * evaluate the function - sum of bell curves * with different means and deviations ************************************************************************ sr2pi = SQRT(2. * ACOS(-1.0)) DO j = 1, 23 width(j) = width(j) / wf height(j) = height(j) / hf END DO tot = 0. DO i = 1, n X(i) = FLOAT(i-1) / (n-1) Y(i) = 23. / FLOAT(n) DO j = 1, 23 Y(i) = Y(i) + (height(j)**2 / sr2pi) * $ EXP( - (X(i)-farey(j))**2 / (2 * width(j))**2 ) END DO tot = tot + Y(i) END DO DO i = 1, n Y(i) = LOG(Y(i) / tot) END DO *----------------------------------------------------------------------- * call the spline interpolator *----------------------------------------------------------------------- ibeg = 0 iend = 0 ifault = 0 CALL CBSPL(X, Y, A, B, C, N, ibeg, iend, alpha, beta, ifault) IF (ifault .NE. 0) PRINT *, 'cbspl: error #', ifault RETURN END ! of gfunc *++++++++++++++++++++++++ End of file rhythm.f +++++++++++++++++++++++++