*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
* 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 +++++++++++++++++++++++++