# 1 Overview

Do the movements musicians’ makes during their performance relate the music?

The short answer is yes they do! How is complicated and requires the use of recurrence quantification analysis (see Marwan, 2008). Briefly, it’s a way of seeing when a time-series is self-similar (recurs).

To show how this analysis works, we will examine a pilot study were I recorded my own postural sway using a Wii BalanceBoard (& Matlab) when I played Bach’s Prelude BWV. 1007

First, we examine if there are any recurrent patterns of movement while I play. Since body movement during performance conveys information to those who watch (see Davidson, 1993, 1994), would expect my movements to have structure and not be just random noise. My movement patterns also should line up with important elements of the musical structure. Let’s examine: 1) a recurrence plot my movements and 2) compare them to recurrence plot of the musical score.

This re-analysis my movement data was created for a two day workshop that Moreno Coco, PhD and I developed to teach recurrence quantification analysis and dynamical systems. Moreno developed the crqa package and helped with this re-analysis.

You can find a conference paper on this data set and also see the studies conducted on professional musicians here.

library(crqa)

# 2 Data

## 2.1 Movement Data

For this section, we use the postural sway movements of me playing Bach’s Cello Suite 1, Prelude (transposed up a 5th for violin), twice. The movements were re-sampled at 12 HZ (which is approximately 83.3 millisecond per unit of time), for a total of 1200 data points. This corresponds roughly to 100 seconds per performance, which is time-Warped so that interval between each 2 bars was made to be equal.

## 2.2 Musical Score Data

BachPrelude.txt = Midi note numbers of all notes. So D and D# are different numbers. Every note in the piece was 16th note, except the last note.

# 3 Analysis of Movement Data

P1 = read.csv("crqa/P1.csv", header = FALSE)
colnames(P1) = colnames(P2) = c("ML", "AP")
P1$ML = as.numeric(as.matrix(P1$ML)); P1$AP = as.numeric(as.matrix(P1$AP))
P2$ML = as.numeric(as.matrix(P2$ML)); P2$AP = as.numeric(as.matrix(P2$AP))
# P1 it generates an extra NA row, which we need to remove
P1 = P1[-1, ]

## 3.2 Different analyses, different hypotheses, the same data.

### 3.2.1 Postural similarity between AP and ML - within performances

This data could be analysed in several different ways, each of which pursuing a particular theoretical direction. We could, for example, be interested on the relation between x (mediolaterial) and y (anterior-posterior) within each performance. There is evidence that AP and ML sway are controlled separately under some conditions (Winter, Prince, Frank, Powell, & Zabjek, 1996), and are coupled when the task demands it (Balasubramaniam, Riley, & Turvey, 2000; Mochizuki, Duarte, Amadio, Zatsiorsky, & Latash, 2006). On the violin, we can assume that AP and ML are weakly coupled, but the degree of coupling depends on the placement of the performer’s feet and either he locks his knees! These initial conditions were constant for my performances (feet at about 30 degree angle from each other and locked knees). If wanted to find when the ML and AP shared patterns, we need to examine the windowed CRQA and see if there a musical driver for the ML an AP system to align.

### 3.2.2 Postural similarity between ML VS ML and AP vs AP - between performances

This question asks how similar are the movements of same performer, playing the same piece. There are two ways to ask this question with recurrence (cross and joint).

#### 3.2.2.1 Cross-Recurrence

• Cross-recurrence will ask the weaker of the two possible questions, “did the same movement pattern repeat across the two performances, regardless of time in the performance.”

• Joint-recurrence asks the more strict question, “did the same movement pattern repeat across the two performances, at exactly the same in each the performance.”

• ML and AP movements might differ this in their relationship the overlap to each other and between performers. AP sway might be more closely tied to the sound-producing movements of the bow, and that ML sway will be more loosely coupled to it as ML sway will provide a clear measure of musician’s expressive and stylistic intentions.

## 3.3 Cross-recurrence between ML and AP across performances.

In order to do RQA on continuous data, we need to figure out what could the best parameters for delay, embedding dimension and radius be. In fact, for each new data set, we start pretty agnostic about the values that should be used to efficiently run crqa(). We examine how the parameters can be unstable (when tested by hand) but we cautiously use optimizeParam to do this.

## 3.4 Finding optimal parameters for each Performance

We calculate AMI and FNN by hand for ML and AP, and also run the optimizeParam on the first performance.

# select a bin size using Sturges' rule (1926)
bin = ceil(1 + log2(nrow(P1)))

# Lets get ML for AP for P1
mutual(P1$ML, bin, 30, plot = TRUE) delay.ML1 = 18 mutual(P1$AP, bin, 30, plot = TRUE)

delay.AP1 = 15

# Select Dimention
plot(false.nearest(P1$ML, m = 6, d = delay.ML1, t = 0, rt = 10, eps = sd(P1$ML)/10))

M.ML1 = 5

plot(false.nearest(P1$AP, m=6, d = delay.AP1, t = 0, rt = 10, eps = sd(P1$AP)/10))

M.AP1 = 5

## Automatic selection
## we choose the parameters of the optimization

par = list(lgM =  20, fnnpercent = 100,
typeami = "mindip",
min.rec = 9, max.rec = 11,
normalize = 0, rescale = 1, mindiagline = 2,
minvertline = 2, tw = 1, whiteline = FALSE,
recpt = FALSE)

Parm.ML1 = optimizeParam(P1$ML, P1$ML, par)
Parm.AP1 = optimizeParam(P1$AP, P1$AP, par)

print(unlist(Parm.ML1))
##  radius  emddim   delay
## 53.4031  4.0000 18.0000
print(unlist(Parm.AP1))
##   radius   emddim    delay
## 57.04129  5.00000 12.00000

A question that that immediately arise is: what do we do with the other performance? Do we need to optimize that one too? Do we need to tailor each analysis to individual trial? If we cap recurrence within a specific range, than we will likely show the exact same recurrence rate.

• Solution to apply ONE set of parameters to all time-series
• So we optimize performance 2 as well and see how much they differ
# Lets get ML for AP for P1
mutual(P2$ML, bin, 30, plot = TRUE) delay.ML2 = 18 mutual(P2$AP, bin, 30, plot=TRUE)

delay.AP2 = 15

# Select Dimention
plot(false.nearest(P2$ML, m = 6, d = delay.ML2, t = 0, rt = 10, eps = sd(P2$ML)/10))

M.ML2 = 5

plot(false.nearest(P1$AP, m = 6, d = delay.AP2, t = 0, rt = 10, eps = sd(P2$AP)/10))

M.AP2 = 5

## Automatic selection

par = list(lgM =  20, fnnpercent = 100,
typeami = "mindip", normalize = 0, rescale = 1, mindiagline = 2,
minvertline = 2, tw = 0, whiteline = FALSE,
min.rec = 9, max.rec = 11,
recpt = FALSE)

Parm.ML2 = optimizeParam(P2$ML, P2$ML, par)
Parm.AP2 = optimizeParam(P2$AP, P2$AP, par)

print(Parm.ML2)
## $radius ## [1] 50.77814 ## ##$emddim
## [1] 4
##
## $delay ## [1] 13 print(Parm.AP2) ##$radius
## [1] 54.4415
##
## $emddim ## [1] 4 ## ##$delay
## [1] 11
• Safe bet for this dataset is to the take the largest lag, and the mode or median of dimensions, and mean of the radius.
• This will allow us to compare between recurrence rates as all the parameters were set the same
delay = max(Parm.ML1$delay, Parm.ML2$delay, Parm.AP1$delay, Parm.AP2$delay)
embed =  median(c(Parm.ML1$emddim, Parm.ML2$emddim, Parm.AP1$emddim, Parm.AP2$emddim))
radius = mean(Parm.ML1$radius, Parm.ML2$radius, Parm.AP1$radius, Parm.AP2$radius)

rescale =  1; normalize = 0; minvertline = 2; mindiagline = 2; whiteline = FALSE;
recpt = FALSE; tw = 0

CRP.ML = crqa(P1$ML, P2$ML, delay, embed, rescale, radius,
normalize, minvertline, mindiagline, tw,  whiteline, recpt)

CRP.AP = crqa(P1$AP, P2$AP, delay, embed, rescale, radius,
normalize, minvertline, mindiagline, tw,  whiteline, recpt)

Take out the results and print them side-by-side. Save the two CRPs for plotting, and convert it into a numeric matrix.

round(cbind(unlist(CRP.ML[1:9]), unlist(CRP.AP[1:9])), digit = 2)
##            [,1]     [,2]
## RR         9.63     9.24
## DET       99.52    99.02
## NRLINE 12457.00 17774.00
## maxL     206.00   146.00
## L         10.11     6.76
## ENTR       2.99     2.56
## rENTR      0.62     0.57
## LAM       99.80    99.68
## TT        11.05     9.35

### 3.4.1 Findings

• RR was higher in ML, that this mean there is more sharing of patterns in ML over AP?
• DET was higher in ML, might expressive movements be more stable?
• L was also higher for ML, again suggesting expressive movements were more stable.

## 3.5 Visualize the two CRP

Save the two CRPs for plotting, and convert it into a numeric matrix.

CRP1 = CRP.ML$RP; CRP2 = CRP.AP$RP

CRP1 = matrix(as.numeric(CRP1), nrow = nrow(CRP1))
CRP2 = matrix(as.numeric(CRP2), nrow = nrow(CRP2)) 

We re-use code shown yesterday to plot the CRQ results of P1 and P2, side by side. We want to reconstruct the time-course in seconds to get a better sense of the scale.

tstamp =  1:nrow(CRP1) # generate the time indeces
realtime = round((tstamp*83.3)/1000) # this is the real time of the the recording.
indnorm = seq(1, max(tstamp), 120) # just take less points to put in the x-axis and y-axis
indnorm = c(indnorm, max(tstamp))

cols = c("white","blue4") # set up the colours

par(mfrow = c(1,2), mar = c(3.8, 3.8, 0.2,2), font.axis = 2, cex.axis = 1,
font.lab = 2, cex.lab = 1.2)

plot(tstamp, tstamp,
type = "n", xaxt = "n", yaxt = "n",
xlab = "", ylab = "")

# l = 1
for (l in 1:ncol(CRP1)){
ind = which(CRP1[,l] == 1)
points(rep(l,length(ind)), ind, cex = 1.2,
col = "blue4", pch = 16)
}

# add the name of the axes
mtext("P1 ML", at = mean(tstamp), side = 1, line = 2, cex = 1.2, font = 2)
mtext("P2 ML", at = mean(tstamp), side = 2, line = 2, cex = 1.2, font = 2)
# add the real time axis
mtext(realtime[indnorm], at = indnorm, cex = 1, font = 2, line = 1, side = 1)
mtext(realtime[indnorm], at = indnorm, cex = 1, font = 2, line = 1, side = 2)

plot(tstamp, tstamp,
type = "n", xaxt = "n", yaxt = "n",
xlab = "", ylab = "")

# l = 1
for (l in 1:ncol(CRP2)){
# print(l)
ind = which(CRP2[,l] == 1)
points(rep(l,length(ind)), ind, cex = 1.2,
col = "blue4", pch = 16)

}

mtext("P1 AP", at = mean(tstamp), side = 1, line = 2, cex = 1.2, font = 2)
mtext("P2 AP", at = mean(tstamp), side = 2, line = 2, cex = 1.2, font = 2)
# add the real time axis
mtext(realtime[indnorm], at = indnorm, cex = 1, font = 2, line = 1, side = 1)
mtext(realtime[indnorm], at = indnorm, cex = 1, font = 2, line = 1, side = 2)

# 4 Generate JRQA

A JRPQ is the RQA of each plot overlapped. So we can simply multiply our RQA plot in R. We must make sure to have Theiler window of at least 1 to ensure the major diagonal has meaning.

There is debate over the size of the Theiler window for data like we have.

## 4.1 JRQA for ML only

# Run RQA for ML
RQA.ML1 = crqa(P1$ML, P1$ML, delay, embed, rescale, radius,
normalize, minvertline, mindiagline, tw = 1,  whiteline, recpt)

RQA.ML2 = crqa(P2$ML, P2$ML, delay, embed, rescale, radius,
normalize, minvertline, mindiagline, tw = 1,  whiteline, recpt)

round(cbind(unlist(RQA.ML1[1:9]), unlist(RQA.ML2[1:9])), digit = 2)
##            [,1]     [,2]
## RR        10.68    10.04
## DET       99.61    99.45
## NRLINE 11590.00 13496.00
## maxL    1145.00  1145.00
## L         12.05     9.71
## ENTR       3.07     2.88
## rENTR      0.63     0.62
## LAM       99.86    99.76
## TT        10.58     9.32
RP.ML1 = RQA.ML1$RP RP.ML2 = RQA.ML2$RP

RP.ML1 = matrix(as.numeric(RP.ML1), nrow = nrow(RP.ML1))
RP.ML2 = matrix(as.numeric(RP.ML2), nrow = nrow(RP.ML2))

# Make JRP for ML
JRP.ML = RP.ML1*RP.ML2

JRP.ML.Plot = matrix(as.numeric(JRP.ML), nrow = nrow(JRP.ML)) 

## 4.2 Visualize JRP for ML

# generate the plot of the JRP
tstamp =  1:nrow(JRP.ML) # generate the time indeces
realtime = round((tstamp*83.3)/1000) # this is the real time of the the recording.
indnorm = seq(1, max(tstamp), 120) # just take less points to put in the x-axis and y-axis
indnorm = c(indnorm, max(tstamp))

cols = c("white","blue4") # set up the colours

par(mar = c(3.8, 3.8, 0.2,2), font.axis = 2, cex.axis = 1,
font.lab = 2, cex.lab = 1.2)

plot(tstamp, tstamp,
type = "n", xaxt = "n", yaxt = "n",
xlab = "", ylab = "")

# l = 1
for (l in 1:ncol(JRP.ML.Plot)){
ind = which(JRP.ML.Plot[,l] == 1)
points(rep(l,length(ind)), ind, cex = 1.2,
col = "blue4", pch = 16)
}

# add the name of the axes
mtext("Time", at = mean(tstamp), side = 1, line = 2, cex = 1.2, font = 2)
mtext("Time", at = mean(tstamp), side = 2, line = 2, cex = 1.2, font = 2)
# add the real time axis
mtext(realtime[indnorm], at = indnorm, cex = 1, font = 2, line = 1, side = 1)
mtext(realtime[indnorm], at = indnorm, cex = 1, font = 2, line = 1, side = 2)

## 4.3 Quantify JRP

We can also quantify the recurrence measures of the JRP plot by inputting the RP obtain again in crqa(). To do so, we need to set the argument recpt = TRUE, and input the RP as the argument ts1, and leave ts2 = NA.

rescale =  1; normalize = 0; minvertline = 2;
mindiagline = 2; whiteline = FALSE;  tw = 0;
ts2 = NA; recpt = TRUE;

resML = crqa(JRP.ML, ts2, delay, embed, rescale, radius,
normalize, minvertline, mindiagline, tw, whiteline, recpt)

print(round(unlist(resML[1:9]), digit = 2))
##      RR     DET  NRLINE    maxL       L    ENTR   rENTR     LAM      TT
##    2.34   98.40 2962.00 1145.00   10.19    2.63    0.64   99.00    6.44

### 4.3.1 Findings

We dont see a lot of recurrence. The time warping on these performances were conducted at the 2 bar level, not the note-by-note level. Had we done note-by-note the results might have looked different. However, to make sense of these we would have to manually examine the musical score and examine JRQA. This analysis would address when does the performer do exactly the same thing at the same time in each performance and how that result might be driven by the underlying music.

# 5 Musical Score

Load the score of midi note numbers.

Score = read.table("crqa/BachPrelude.txt", header = FALSE)
Score = as.numeric(ScoreV1) The problem is going to be that we do not know how long each note was held for, but we do know how long each 2 bars were held for (5 seconds per 2 bars x 12 HZ = 30 samples per bar). The score has 16 notes per bar that means we cannot overlap the two signals as 30/16 does not create an interger ratio. However, we can simply resample the data to make it work. [Note this is a terrible solution, but the idea is that you need to align he music and movement. A better method is code how long each note was played and ensure you match the notes to fit that length] #install.packages("signal") require(signal) # a way to force R to install the package if missing # we need a new dataframe as our orginal will not match the new size P1R <- NULL P2R <- NULL P1RML <- resample(P1$ML,13.46,12) P2R$ML <- resample(P2ML,13.46,12) So we need to upconvert the time on the notes to match the time length of the time-series [again this not a proper alignment, but an exmaple of the process] length(Score) ## [1] 673 Score.long<-rep(Score, each = 2)  # 6 Overlap RQA of ML with RQA of Score First we need to regenerate the RQA for larger time-series par = list(lgM = 20, fnnpercent = 100, radiusspan = 100, radiussample = 20, typeami = "mindip", min.rec = 9, max.rec = 11, normalize = 0, rescale = 1, mindiagline = 2, minvertline = 2, tw = 1, whiteline = FALSE, recpt = FALSE) Parm.ML1.RS = optimizeParam(P1RML, P1R$ML, par) Parm.ML2.RS = optimizeParam(P2R$ML, P2R$ML, par) delay.2 = max(Parm.ML1.RS$delay, Parm.ML2.RS$delay) embed.2 = median(c(Parm.ML1.RS$emddim, Parm.ML2.RS$emddim)) radius.2 = mean(Parm.ML1.RS$radius, Parm.ML2.RS$radius) rescale = 1; normalize = 0; minvertline = 2; mindiagline = 2; whiteline = FALSE; recpt = FALSE; tw = 1 RQA.ML1.RS = crqa(P1R$ML, P1R$ML, delay.2, embed.2, rescale, radius.2, normalize, minvertline, mindiagline, tw=1, whiteline, recpt) RQA.ML2.RS = crqa(P2R$ML, P2R$ML, delay.2, embed, rescale, radius.2, normalize, minvertline, mindiagline, tw=1, whiteline, recpt) RQA.ML1.RS = RQA.ML1.RS$RP
RQA.ML2.RS = RQA.ML2.RS$RP RQA.ML1.RS = matrix(as.numeric(RQA.ML1.RS), nrow = nrow(RQA.ML1.RS)) RQA.ML2.RS = matrix(as.numeric(RQA.ML2.RS), nrow = nrow(RQA.ML2.RS))  ## 6.1 Lets RQA the score Remember RQA loses data: • When we time-lag and embedded we lost 57 number data points (delay X embed-1) • Thus be before we multipy the music RQA X Movement RQA we need to trim the music RQA matrix trim = delay.2*(embed.2 - 1) # the difference due to having different delays and embedding dimensions used to run RQA on the posture and on the music score recpt = FALSE; ## reset this back to FALSE Score = Score.long[1:(length(Score.long)-trim)] RQA.Score = crqa(Score, Score, d = 1, e = 1, rescale = 1, radius = .00001, normalize, minvertline, mindiagline, tw = 1, whiteline, recpt) round(unlist(RQA.Score[1:9]), digit = 1) ## RR DET NRLINE maxL L ENTR rENTR LAM TT ## 8.8 53.7 33590.0 34.0 2.3 0.5 0.2 99.1 2.1 RQA.Score = RQA.Score$RP
RQA.Score = matrix(as.numeric(RQA.Score), nrow = nrow(RQA.Score))

## plot the RQA of the Score

tstamp =  1:nrow(RQA.Score) # generate the time indeces
realtime = round((tstamp*83.3)/1000) # this is the real time in seconds of the the recording.
indnorm = seq(1, max(tstamp), 120) # just take less points to put in the x-axis and y-axis

cols = c("white","blue4") # set up the colours

par(mar = c(3.8, 3.8, 0.2, 2), font.axis = 2, cex.axis = 1,
font.lab = 2, cex.lab = 1.2)

plot(tstamp, tstamp,
type = "n", xaxt = "n", yaxt = "n",
xlab = "", ylab = "")

# l = 1
for (l in 1:ncol(RQA.Score)){
ind = which(RQA.Score[,l] == 1)
points(rep(l,length(ind)), ind, cex = .2,
col = "blue4", pch = 16)
}

# add the name of the axes
mtext("Time (Score)", at = mean(tstamp), side = 1, line = 2, cex = 1.2, font = 2)
mtext("Time (Score)", at = mean(tstamp), side = 2, line = 2, cex = 1.2, font = 2)
# add the real time axis
mtext(realtime[indnorm], at = indnorm, cex = 1, font = 2, line = 1, side = 1)
mtext(realtime[indnorm], at = indnorm, cex = 1, font = 2, line = 1, side = 2)

### 6.1.1 JRP between ML and the Score

# generate the JRP plot which integrates the RPs of postures (ML1, ML2) and music played
JRP.Score.ML1 <- RQA.Score*RQA.ML1.RS
JRP.Score.ML2 <- RQA.Score*RQA.ML2.RS

### 6.1.2 Quantify and examine JRP

We quantify the recurrence measures of the JRP plot between the performances (as expressed in the ML sway) and the score (as expressed in midi notes) by inputting the RP obtain again in crqa(). To do so, we need to set the argument recpt = TRUE, and input the RP as the argument ts1, and leave ts2 = NA.

rescale =  1; normalize = 0; minvertline = 2; mindiagline = 2;
whiteline = FALSE; recpt = TRUE; tw = 0; ts2 = NA

resML1 = crqa(JRP.Score.ML1, ts2, delay, embed, rescale, radius,
normalize, minvertline, mindiagline, tw, whiteline, recpt)
resML2 = crqa(JRP.Score.ML2, ts2, delay, embed, rescale, radius,
normalize, minvertline, mindiagline, tw, whiteline, recpt)

cbind(round(print(unlist(resML1[1:9])), digit = 2),
round(print(unlist(resML2[1:9])), digit = 2))
##           RR          DET       NRLINE         maxL            L
##    1.0916504   47.4252950 3578.0000000   20.0000000    2.4041364
##         ENTR        rENTR          LAM           TT
##    0.6383032    0.2568721   86.6578454    2.0714286
##           RR          DET       NRLINE         maxL            L
##    0.9675472   45.4590694 3134.0000000   11.0000000    2.3318443
##         ENTR        rENTR          LAM           TT
##    0.5720987    0.2484593   84.6230406    2.0559166
##           [,1]    [,2]
## RR        1.09    0.97
## DET      47.43   45.46
## NRLINE 3578.00 3134.00
## maxL     20.00   11.00
## L         2.40    2.33
## ENTR      0.64    0.57
## rENTR     0.26    0.25
## LAM      86.66   84.62
## TT        2.07    2.06
## plot the JRP plots
tstamp =  1:nrow(JRP.Score.ML1) # generate the time indeces
realtime = round((tstamp*83.3)/1000) # this is the real time in seconds of the the recording.
indnorm = seq(1, max(tstamp), 120) # just take less points to put in the x-axis and y-axis

cols = c("white","blue4") # set up the colours

par(mfrow = c(1,2), mar = c(3.8, 3.8, 0.2, 2), font.axis = 2, cex.axis = 1,
font.lab = 2, cex.lab = 1.2)

plot(tstamp, tstamp,
type = "n", xaxt = "n", yaxt = "n",
xlab = "", ylab = "")

# l = 1
for (l in 1:ncol(JRP.Score.ML1)){
ind = which(JRP.Score.ML1[,l] == 1)
points(rep(l,length(ind)), ind, cex = 1.2,
col = "blue4", pch = 16)
}

# add the name of the axes
mtext("Time (ML1*Score)", at = mean(tstamp), side = 1, line = 2, cex = 1.2, font = 2)
mtext("Time (ML1*Score)", at = mean(tstamp), side = 2, line = 2, cex = 1.2, font = 2)
# add the real time axis
mtext(realtime[indnorm], at = indnorm, cex = 1, font = 2, line = 1, side = 1)
mtext(realtime[indnorm], at = indnorm, cex = 1, font = 2, line = 1, side = 2)

plot(tstamp, tstamp,
type = "n", xaxt = "n", yaxt = "n",
xlab = "", ylab = "")

# l = 1
for (l in 1:ncol(JRP.Score.ML2)){
# print(l)
ind = which(JRP.Score.ML2[,l] == 1)
points(rep(l,length(ind)), ind, cex = 1.2,
col = "blue4", pch = 16)

}

mtext("Time (ML2*Score)", at = mean(tstamp), side = 1, line = 2, cex = 1.2, font = 2)
mtext("Time (ML2*Score)", at = mean(tstamp), side = 2, line = 2, cex = 1.2, font = 2)
# add the real time axis
mtext(realtime[indnorm], at = indnorm, cex = 1, font = 2, line = 1, side = 1)
mtext(realtime[indnorm], at = indnorm, cex = 1, font = 2, line = 1, side = 2)

### 6.1.3 Findings

We can see recurrent patterns between the ML movement and the score played. Even with this rough approximation, there are interesting places which we would want to follow up. For example, the diagonal lines you see in the plot towards the end of the score reflect the location in the music where the musician is running up the scale with pedal notes in between.