*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
* OSCILL - a frequency analyzer for pulse trains
* by Andy Allinger, 2011, released to the public domain.
* This program may be used by any person for any purpose.
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
* 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 Ratio of frequencies of adjacent oscillators
* frang REAL In Ratio of min to max frequency of a single
* oscillator
* A[o,0:n] REAL Out Amplitude
* P[o,0:n] REAL Out Phase
* F[o,0:n] REAL Out Frequency
* fault INTEGER Out nonzero on error
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
SUBROUTINE oscill (T,S,m,n,o,step,fdeep,fspac,frang,A,P,F,fault)
IMPLICIT NONE
INTEGER m, n, o, fault
REAL T, S, step, fdeep, fspac, frang, A, P, F
DIMENSION T(m), S(m), A(o,0:n), P(o,0:n), F(o,0:n)
* constants
REAL L, LS, PS, rmin, rmax, twopi, v
PARAMETER (
* & g = 7.0, ! sharpening coefficient
& L = 3.0, ! half-life of amplitude decay, cycles
& LS = 0.25, ! half-life of satiation decay, cycles
& PS = 0.618034, ! max duration of ignoring, cycles
& rmin = .0001, ! distance to stay away from origin
& rmax = .9999, ! closest approach to unit circle
& twopi = 6.283185307179586,
& v = 1.0 ) ! volume of the beaker
* local variables
INTEGER i, j, k ! counters
REAL
& b, ! amplitude decay rate
& b1, ! amplitude decay for whole step
& c, ! salience decay rate
& dt, dx, dy, ! changes in t, x, and y
& fmin, f0, fmax, ! min, central, max freq. each oscillator
& nu ! frequency, in Hertz
REAL
& phi, ! phase, measured in beats elapsed
& pprev, ! phase at previous event
& rho, r2, ! amplitude, amplitude squared
& s1, sprev, ss, ! salience: present, previous, accumulation
& SNAP, SSIN, SCOS, SARC, ! replacements for EXP, SIN, COS, ATAN2
& 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, ! weight coefficient of an event
& x, y ! x = rho @ cos(theta) ; y = rho @ sin(theta)
DOUBLE PRECISION
& dmachine, DEPS, ! very small number
& b2, st0, st1, st2, sp, stp, det ! decay factor, least-squares sums
LOGICAL e ! whether exponentiation is required
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
* begin
fault = 0 ! clear error codes
DEPS = dmachine(1)
* validate
IF (m < 1 .OR. n < 1 .OR. o < 1) THEN
fault = 1
RETURN
END IF
* check that enough timepoints have been allocated for the given events
IF (n < T(m) / step) THEN
fault = 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)
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
* for each oscillator
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
DO k = 1, o
f0 = f0 * fspac ! central frequency
fmin = fmin * fspac ! minimum frequency
fmax = fmax * fspac ! maximum frequency
e = .FALSE. ! no stimulus event intervening between steps
b = LOG(0.5) * f0 / L ! decay rate where half-life is L cycles
b1 = EXP(b * step)
c = LOG(0.5) * f0 / LS ! longer memory for slow oscillators
rho = rmin ! minimal amplitude
phi = 0. ! initial phase
nu = f0 ! default frequency
A(k,0) = rho
P(k,0) = phi ! pass thru origin
F(k,0) = f0
i = 1
t0 = 0.
tnext = 0.
tprev = 0.
ss = 0.
sprev = 0.
pprev = 0.
st0 = 0.
st1 = 0.
st2 = 0.
sp = 0.
stp = 0.
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
* 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(b * dt) ! approximate exponential decay
rho = MIN(MAX(rmin, rho), rmax)
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
* ignore stimuli that occur too close together
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
s1 = S(i)
ss = ss * SNAP(c * (tt-t0))
IF (s1 < ss .AND. phi-pprev < PS .AND. s1 < sprev) THEN
i = i + 1
tprev = tt
e = .TRUE.
GO TO 10
ELSE
ss = ss + s1
sprev = s1
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
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
w = (s1 * u) / SQRT((v-x) **2 + y **2)
dx = w * (v-x)
dy = w * (-y)
x = MIN(x + dx, v)
IF (ABS(dy) .GE. ABS(y)) THEN ! don't overcorrect
IF (y < 0.) phi = phi + 1. ! increment phase if x-axis is reached
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
b2 = SNAP(b * (tt-t0)) ! decay factor
st0 = st0 * b2 + 1 ! SUM t^0
st1 = st1 * b2 + tt ! SUM t^1
st2 = st2 * b2 + DPROD(tt, tt) ! SUM t^2
sp = sp * b2 + phi ! SUM phi
stp = stp * b2 + DPROD(tt, phi) ! SUM t @ phi
det = st0 * st2 - st1 * st1 ! determinant of normal equations
IF (det > DEPS) nu = SNGL( (st0 * stp - sp * st1) / det )
nu = MIN(MAX(fmin, nu), fmax) ! impose limits
i = i + 1
t0 = tt
tprev = tt
pprev = phi
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(b * dt) ! approximate exponential function
phi = phi + nu * dt
ELSE
rho = rho * b1 ! avoid an exponentiation when possible
phi = phi + nu * step
END IF
A(k,j) = rho
P(k,j) = phi
F(k,j) = nu
tprev = tnext
e = .FALSE.
END DO ! next j
END DO ! next k
RETURN
END ! of oscill
*++++++++++++++++++++++++ end of file oscill.f +++++++++++++++++++++++++