IEEE/ACM TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING, VOL. XX, NO. XX, XX 2020

1

Bayesian Melody Harmonization
Based on a Tree-Structured Generative Modelof Chord Sequences and Melodies
Hiroaki Tsushima, Eita Nakamura, and Kazuyoshi Yoshii Member, IEEE

Abstract—This paper describes a melody harmonization
method that generates a sequence of chords (symbols and onset
positions) for a given melody (a sequence of musical notes). A
typical approach to melody harmonization is to use a hidden
Markov model (HMM) that represents chords and notes as latentand observed variables, respectively. This approach, however,does not consider the syntactic functions (e.g., tonic, dominant,and subdominant) and hierarchical structure of chords that playvital roles in traditional harmony theories. In this paper, we
propose a unified hierarchical generative model consisting ofa probabilistic context-free grammar (PCFG) model generating
chord symbols associated with syntactic functions, a metrical
Markov model generating chord onset positions, and a Markov
model generating a melody conditioned by a chord sequence. To
estimate a musically natural tree structure, the PCFG is trainedin a semi-supervised manner by using chord sequences with tree
structure annotations. Given a melody, a sequence of a variable
number of chords can be estimated by using a Markov chain
Monte Carlo method that partially and iteratively updates the
symbols, onset positions, and tree structure of chords according
to the posterior distribution of chord sequences. Experimental
results show that the proposed method outperformed the HMMbased method and a conventional rule-based method in terms ofpredictive abilities.
Index Terms—Melody harmonization, probabilistic contextfree grammar, Markov model.

I. I NTRODUCTION
HORD arrangement is one of the most important tasks
in the composition of popular music because chord
sequences characterize musical moods and styles. To helpmusically untrained people compose original pieces, automatic
melody harmonization, i.e., automatic generation of a chord
sequence for a given melody (note sequence), has been studied[1]–[14]. For this, it is necessary to computationally characterize acceptable chord sequences and the relationships between
chords and melody notes. Statistical modeling is effective forC

Manuscript received XXXX XX, 2019; revised XXXX XX, 2019; accepted
XXXX XX, 2020. Date of publication XXXX XX, 2020; date of current
version XXXX XX, 2020. This work was partially supported by JSPS
KAKENHI Nos. 26700020, 16H01744, 16J05486, 16H02917, 16K00501, and19K20340, JST ACCEL No. JPMJAC1602, and Kayamori Foundation. The
associate editor coordinating the review of this manuscript and approving it
for publication is Prof. XXX YYY (Corresponding author: Eita Nakamura).
H. Tsushima, E. Nakamura, and K. Yoshii are with Graduate School ofInformatics, Kyoto University, Kyoto 606-8501, Japan, email: {tsushima,
enakamura, yoshii}@sap.ist.i.kyoto-u.ac.jp. E. Nakamura is also with the
Hakubi Center for Advanced Research, Kyoto University, Kyoto 606-8501,
Japan.
Color versions of one or more of the figures in this paper are available
online at http://ieeexplore.ieee.org.
Digital Object Identifier XXXX

Probabilistic
context-freegrammar (PCFG)
for chord symbols













C

Am

F

G

Metrical Markov model
for chord onset positions
Markov model
for musical notes




C

Am

F



G



Fig. 1. A tree-structured hierarchical generative model that stochasticallygenerates a chord sequence and a melody in this order. The labels Si represent
syntactic categories behind chord symbols explained in Section III.

data-driven induction of such musical grammar and can beapplied to various musical styles [7], [10], [11].
A typical approach for melody harmonization based on
statistical modeling is to integrate a probabilistic model representing a chord sequence with one representing the relationships between chords and a melody. If the former is a
Markov or semi-Markov model and the latter is a generativemodel of a melody given a chord sequence, the integrated
model is described as a hidden Markov model (HMM) [7]or hidden semi-Markov model (HSMM) [10]. In this model,
the sequential dependency of chord symbols is described withtransition probabilities, which can be learned from data.
This approach, however, does not consider that some chords
play a similar syntactic role and that music has a hierarchicalstructure (phrase, section, etc.) [15]–[17]. In harmony theories,
the syntactic roles of chord symbols are known as harmonic
functions such as tonic (T), dominant (D), and subdominant(SD)1 . The hierarchical structure of chords or harmonic functions is often represented as a tree [17]–[21]. For example, achord sequence (C, Dm, G, Am, C, F, G, C) in C major key
can be interpreted as a function sequence (T, SD, D, T, T, SD,
D, T), which is further parsed as a binary tree ((T, ((SD, D),T)), (T, ((SD, D), T))), where subtrees such as (SD, D) and (T,((SD, D), T)) appear repeatedly in different abstraction levels.
It has been shown that statistical models such as an HMM anda probabilistic context-free grammar (PCFG) model that canrepresent the syntactic functions and hierarchical structure ofchords outperform Markov models in predictive ability [22].
1 In the simple cases consisting of triads, these functions are associated withthe three main triads (C, G, and F in C major key) and their relative minor
chords.

2

IEEE/ACM TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING, VOL. XX, NO. XX, XX 2020

Another problem with the conventional Markov or semiMarkov models is that the metrical structure of chord onset
positions cannot be modeled. Chord transitions are considered
for a fixed time step (e.g., 16th note) and in semi-Markov
models chord durations are explicitly modeled. While chord
onsets tend to appear more frequently at strong beats, such
metrical structure cannot be described with these conventional
Markov or semi-Markov models.
In this paper we propose a statistical melody harmonization
method incorporating syntactic functions, tree structure, andthe metrical structure of chords, aiming to generate musicallymeaningful chord sequences. We formulate a unified hierarchical generative model that consists of (1) a PCFG model
representing the generative process of chord symbols, (2) ametrical Markov model [23], [24] representing that of chord
onset positions, and (3) a Markov model representing that of anote sequence from a chord sequence (Fig. 1). The harmonic
functions and hierarchical structures of chords are represented
by the leaf and inside nodes of a latent tree behind chords.
The metrical structure of chord rhythms is represented by the
metrical Markov model. Given a melody, chord sequenceswith a variable number of symbols can be sampled from theposterior distribution by using a Markov chain Monte Carlo(MCMC) method.
For application of statistical models, it is crucial to train
the parameters properly from data. In our previous work, wedeveloped an unsupervised learning method for the PCFG
model and experimentally showed the potential of this model
for melody harmonization [12], [13]. The estimated latent trees
behind chord sequences, however, often mismatch human’s
structural interpretations and the leaf nodes of the trees failto represent the harmonic functions of chords. This is mainly
because the PCFG learning tends to easily get stuck in bad
local optima due to its strong sensitivity to initialization. It
is thus hard to induce musically meaningful latent trees from
non-annotated chord sequences like (C, Dm, G, Am, C, F, G,
C) in an unsupervised manner.
To solve this problem, we propose two refined methods forlearning the PCFG model. First, we develop an initialization
method to make the model’s grammar more similar to human’sinterpretation. We initialize the PCFG model by using an
HMM that represents harmonic functions and chords as latentand observed variables, respectively (Section III-C3). Second,
we develop a learning method using additional data to makethe model predict latent trees similar to human’s structural
annotations. We train the PCFG model in a semi-supervisedmanner by using weakly annotated chord sequences like ((C,((Dm, G), Am)), (C, ((F, G), C))) as explained in Section
III-C2. We confirm that these methods contribute to estimating
latent trees behind chords that represent hierarchical musicalstructure (Section V).
The major contribution of this study is to present a Bayesianinference method for the unified probabilistic model of chords
and melodies, which is an extension of a previous study[12]. Another contribution is a comprehensive report of theexperimental evaluation results. We also show the effectivenessof the initialization and semi-supervised learning methods bysystematic comparisons. Sections III-C2 and III-C3 present

new technical improvements and Section V presents newexperimental results. The other sections describe the method
proposed in [12] in full detail.
In Section II of this paper, we review related work on
melody harmonization and music language modeling. Section III formulates the proposed model and Section IV presents
the melody harmonization method based on the model. Section V reports the experimental evaluations. We conclude witha brief summary and mention of future work in Section VI.
II. R ELATED W ORKThis section reviews related studies on melody harmonization and on music language models of chords and melodies.
A. Melody Harmonization
Melody harmonization systems can be roughly categorized
into two types in terms of their objectives. In the first type ofsystems, the aim is to generate a sequence of chord symbols
given a melody, and in the second type of systems, the aim isto generate multiple voices of musical notes.
As a study on the first type of systems, Chuan and Chew [6]proposed a method that selects musical notes from a given
melody that are likely to form chords by using a support
vector machine, constructs triads from the selected notes,
and makes a chord sequence by using hand-crafted rules. Simon et al. [7] developed a commercial system called MySong
based on an HMM representing sequential chord transitions,and Raczyński et al. [11] proposed a similar Markov model
representing chords conditioned by melodies and time-varying
keys. Tsushima et al. [12] proposed a harmonization method
based on a PCFG model representing the hierarchical structureof chords and a Markov model representing the transitions of
melody notes. De Prisco et al. [9] proposed a harmonization
method for bass-line inputs based on a distinctive network that
represents the relationships between bass notes, the previous
chord, and the current chord.
Recently, a deep neural network (DNN)-based method for
melody harmonization was proposed by Lim et al. and wasshown to outperform an HMM-based method in both objective
and subjective evaluation metrics [25]. An advantage of ourapproach over such DNN-based approaches is that the proposed model can learn and represent the grammatical structureof music in a similar way as humans do, which is consideredessential for enhancing the directability of interactive composition systems [13]. In contrast, it is generally difficult toanalyze or interpret the internal structure of DNNs.
As a study on the second type of systems, Ebcioğlu [1]proposed a rule-based method for generating four-part chorales
in Bach’s style. Several methods using variants of genetic algorithms based on music theories have also been proposed [2],[3], [8]. Allan and Williams [4] proposed an HMM-based
method that represents chords as latent variables and notes
as observed outputs. An HSMM has been used for explicitly
representing the durations of chords [10]. Paiement et al. [5]proposed a hierarchical tree-structured model that describes
chord movements from the viewpoint of hierarchical timeH. TSUSHIMA et al.: BAYESIAN MELODY HARMONIZATION

3scales by dividing the notations of chords. To generate fourpart chorales, a deep recurrent neural network has also beenused for capturing the long-term dependency characteristics of
melodies and harmonies [14].
B. Statistical Modeling of Melodies and Chords
Statistical models of melodies play an important role in thetasks of automatic harmonization and melody generation. The
standard approach for modeling sequences of musical notes(melodies) is to use Markov models [26]. Pachet et al. [27]proposed an efficient approach for generating a melody basedon Markov models with constraints. To describe the generative
process of melodies conditioned on chords, Simon et al. [7]proposed a generative model based on HMMs and bag-ofwords that considers chord symbols as latent variables andnotes as observed words.
In the generative theory of tonal music (GTTM) [17], a notesequence is assumed to form a tree that describes the degreesof the relative importance of individual notes. This theory consists of many hand-crafted rules used for recursively reducing anote sequence to a single note. Computational implementationof the GTTM and its application to music analysis have beenstudied [28], [29]. A probabilistic formulation of the GTTM
based on a PCFG was recently proposed for learning the
production rules of the PCFG from melodies [30].
Statistical models of chord sequences have been actively
investigated for automatic chord estimation [31]–[33], musicalgrammar analysis [18], [19], and automatic music arrangement[21], [34]. The most popular approach is to use a Markov (orn-gram) model [31], [34]. To avoid the sparseness problem
with a large value of n, various smoothing methods have been
proposed [35]. Yoshii et al. [32] proposed a vocabulary-free
infinity-gram model that represents the dependency of a chordon a variable-length history of chords.
To make generated chord sequences musically natural andinterpretable, it is considered effective to describe the musical properties of chords (e.g., syntactic categories, cadential
properties, and hierarchical structure) as in harmony theories. Tsushima et al. [22] attempted unsupervised learning of
HMMs for spontaneously discovering the syntactic categoriesof chords. As a result, the learned syntactic categories corresponded to the harmonic functions and the trained HMM
outperformed Markov models in predictive ability.
To capture the syntactic structure behind chord sequences,tree-structured models have been studied. Paiement et al. [33]formulated several hidden layers of state transitions for representing the tree structure behind chords. Some studies
attempted to explicitly describe the generative grammar torepresent the hierarchical structure of chords [18]–[21]. Steedman [18] and Rohrmeier [19] proposed CFG-based production
rules for chord sequences. A probabilistic extension was later
studied for music arrangement [21]. In these studies, a listof non-terminal symbols and that of production rules weremanually given based on music theories or musical intuition.
III. P ROBABILISTIC M ODELINGThis section explains the probabilistic model of chords andmelody notes. After introducing some mathematical notations








 

Leaf nodes









 

 













Fig. 2. Mathematical notation for tree structure t behind chord symbols z
based on a PCFG model. In this example, the number of chords N is 6 andS represents the start symbol.
TABLE IM ATHEMATICAL NOTATION IN THIS PAPER
SymbolMeaning

SectionN
T

Number of chord symbols
Number of 16th-note-level time units

III-A
III-A

z
y
x
o

Chord symbols
Chord onset positions
Pitches of melody notes
Onset positions of melody notes

III-A
III-A
III-A
III-A

t
t̂

Tree
Tree structure without node labels

III-B1
III-C2

Number of distinct non-terminal symbols
Parameters of PCFG
Transition probabilities over chord onset positions
Transition probabilities over melody pitches

III-B1
III-B1
III-B2
III-B3K
θ, η, λπτin Section III-A, the model is formulated in Section III-B.
Methods for training the model parameters are described in
Section III-C. For simplicity of notation, we consider thatdata consist of only one pair of a chord sequence and a
melody. Extension for the case with multiple sequences isstraightforward.
A. Mathematical Notation
A list of the mathematical symbols is provided in Table I.
In this paper, we assume that musical pieces have the time
signature of 4/4 and the onset positions of chords and melody
notes are on the 16th-note-level grid. Let L be the number ofmeasures of a musical piece (L = 8 in this paper) and T =16L be the total number of time units. A sequence of chordN
symbols and their onset positions are denoted by z = {zn }n=1
N
and y = {yn }n=1 , respectively, where N is the number ofchords and each yn takes an integer in [0, T ). Similarly, thepitches and onset positions of melody notes in the region ofnnchord zn are denoted by xn = {xn,i }Ii=1
and on = {on,i }Ii=1,
respectively, where In is the number of melody notes in theregion, xn,i is a MIDI note number from 33 (A1) to 92 (G#6),
and on,i takes an integer in [yn , yn+1 ). The whole melody
PN
N
is denoted by x = {xn }n=1 and y. Let I =n=1 In be
the total number of melody notes. We use the integer pitch
representation instead of the spelled pitch representation andtreat enharmonic notes equivalently. Since most notes in the4IEEE/ACM TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING, VOL. XX, NO. XX, XX 2020

Dirichlet Dirichlet

ξ

ζ

θ

η

Chord tree

t

Beta

ν
Dirichlet

PCFG
parametersβλDirichlet

Chord symbolszMetrical Markov
model parametersπChord onset
positionsyMelodypitchesxγp(yn |yn−1 ) = πyn−1 mod16,yn −yn−1 ,Markov model
parametersτFig. 3. Dependencies of variables and parameters of the generative model.

popular music data we use in this study are on the diatonic
scale, the information lost in the integer pitch representation
is little and is considered to have a small effect.
B. Model Formulation
We formulate a unified hierarchical generative model ofchord symbols z, chord onset positions y, and melody pitches
x (Fig. 3). This model consists of three sub-models: (1) a
PCFG model of z, (2) a metrical Markov model of y, and (3)a Markov model of x conditioned by z and y.
1) PCFG Model of Chord Symbols: Following the ideaof previous studies [18], [19], [21], we assume that chordsymbols z are derived from a latent tree t in a similar way thatwords are derived from a syntactic tree in the phrase structure
grammar of natural language. A derivation tree t and chordsymbols z are generated in this order according to a PCFG
G = (V, Σ, R, S) defined by a set of non-terminal symbols
V = {Sk }K
k=1 , a set of terminal symbols (chord symbols) Σ,a set of rule probabilities R, and a start symbol S (a nonterminal symbol at the root of t), where K is the numberof distinct non-terminal symbols. The non-terminal symbols
are expected to represent the syntactic roles and hierarchicalstructures of chord symbols. There are three types of ruleprobabilities. θ(A → BC) is the probability that a nonterminal symbol A ∈ V branches into non-terminal symbols
B ∈ V and C ∈ V . η(A → α) is the probability thatA ∈ V emits terminal symbol α ∈ Σ. λA ∈ [0, 1] is theprobability that a non-terminal symbol A ∈ V emits a terminal
symbol, and otherwise it branches. The first two probabilities
are normalized as follows:XX
θ(A → BC) = 1,
η(A → α) = 1.
(1)
B,C∈V

We also use the notation tm:n to indicate the root node of
the subtree for simplicity. In this paper, we refer to the nodes{tn:n }Nn=1 that emit chord symbols as leaf nodes.
2) Metrical Markov Model of Chord Onset Positions: The
metrical Markov model [23], [24] of chord onset positions y
on the 16th-note-level grid is defined as follows:α∈ΣA tree and a chord sequence are generated by a cascading
process of generating non-terminal symbols from the startsymbol and finally terminal symbols from the non-terminal
symbols. The probability of the tree and chord sequence isgiven as a product of relevant branching and emitting probabilities. We write θA = {θ(A → BC)}B,C∈V , θ = {θA }A∈V ,ηA = {η(A → α)}α∈Σ , η = {ηA }A∈V , and λ = {λA }A∈V .
Similar notations are used throughout this paper.
A subtree of t that derives zm:n = {zm , zm+1 , ..., zn } isdenoted by tm:n (Fig. 2). In particular, we have t = t1:N .

(2)

where πa,b indicates the probability that a chord has an onset
at the a-th position in a measure (0 ≤ a < 16) and a durationof b time units (0 < b ≤ T ). We write πa = {πa,b }Tb=03) Markov Model of Melody Pitches: The Markov model of
melody pitches x conditioned by chord symbols z and chord
onset positions y is defined as follows:p(xn,1 |xn−1,In−1 , zn ) = τ zn (xn−1,In−1 , xn,1 ),(3)p(xn,i |xn,i−1 , zn ) = τ zn (xn,i−1 , xn,i ) (2 ≤ i ≤ In ), (4)where τ c (a, b) is the transition probability from pitch a topitch b under chord symbol c, that is, the probability that pitch
a would be generated following pitch b in the time span ofchord symbol c. We write τac = {τ c (a, b)}92
b=33 (33 = A1 and
92 = G#6).
4) Hierarchical Bayesian Integration of Three Sub-Models:Let Ω = {t, z, y, x} be the set of the latent and observed
variables and Θ = {θ, η, λ, π, τ } be the set of model
parameters. Assuming that chord symbols and melody pitches
are generated independently, the whole model is given byp(Ω, Θ) = p(Ω|Θ)p(Θ)= p(t, z|θ, η, λ)p(y|π)p(x|z, y, τ )p(Θ),(5)

where p(t, z|θ, η, λ) is the probability of chords and a tree
given by the PCFG model, p(y|π) is the probability ofchord onset positions given by the metrical Markov model,p(x|z, y, τ ) is the probability of melody pitches given by
the Markov model, and p(Θ) = p(θ)p(η)p(λ)p(π)p(τ ) is aprior distribution over the model parameters. We consider priordistributions on the model parameters to make the learning
process more efficiently (avoiding getting stuck in bad local
optima) [22] and to induce the model’s probability parameters
sparser, which tends to represent syntactic structure similar
to human’s interpretation [36]. To make Bayesian inference
tractable, we use conjugate Dirichlet and beta priors:θA ∼ Dir(ξA ), ηA ∼ Dir(ζA ), λA ∼ Beta(νA ),(6)

τac ∼ Dir(γac ),(7)

πa ∼ Dir(βa ),γacwhere ξA , ζA , νA , βa , andare hyperparameters. It is
known that smaller values of these hyperparameters make theprobability parameters sparser.
C. Model Training
The three sub-models described in Section III-B can betrained separately using different types of training data. For
training the PCFG model, the basic method [12] of unsupervised learning is explained in Section III-C1. As refinements
to this method, we propose a semi-supervised learning method
using tree structure annotations on z in Section III-C2 and an

H. TSUSHIMA et al.: BAYESIAN MELODY HARMONIZATION

5

Semi-supervised learning of PCFG

Unsupervised learning of PCFG

Q

Q

Q

p(l, B, C)

Q
Q

Q

Q

= p(tn:n+l−1 = B, tn+l:n+k = C | tn:n+k = A, zn:n+k )

Q

Q

Q

Q

Q

Q

Q

Q

Q

Q

Q

Q

Q

Q

Q

Q

Q

C

Dm

G

Am

C

F

G

C

C

Dm

G

Am

C

F

G

C

Q

Q
Q

Q

Q

C
A
= (1 − λA )θ(A → BC) pB
n,n+l−1 pn+l,n+k /pn,n+k , (11)where l ∈ {1, . . . , k} indicates a split position.
Finally, we sample parameters θ, η, and λ according top(θ, η, λ|t, z) = p(θ|t, z)p(η|t, z)p(λ|t, z) as follows:Q

Q

C derived from tn:n+k are then sampled according toQθA | t, z ∼ Dir(ξA + uA ),

(12)

































ηA | t, z ∼ Dir(ζA + vA ),

(13)

C

Dm

G

Am

C

F

G

C

C

Dm

G

Am

C

F

G

C

λA | t, z ∼ Beta(νA + wA ),(14)

HMM-based initialization of leaf nodes

Fig. 4. Four possible learning methods for the PCFG model. The left andright columns show unsupervised and semi-supervised learning methods,
respectively. The top and bottom rows show the cases with and without
an HMM-based initialization method. Question marks indicate the variables
estimated in unsupervised learning.

initialization method based on an HMM of z in Section III-C3,
to make a tree t estimated by the PCFG model close to humaninterpretation. Combining these two refinements, we have four
possible learning methods shown in Fig. 4. In Sections III-C4
and III-C5, we explain how to train the metrical Markov modelof chord onset positions y and the Markov model of melody
pitches x.
1) Unsupervised Learning of PCFG: Our goal is to obtainthe maximum-a-posteriori (MAP) estimate of the model parameters Θ = {θ, η, λ, π, τ }. To estimate the parameters θ,η, and λ from chord symbols z in an unsupervised manner, we
draw samples from the posterior distribution p(θ, η, λ, t|z) byusing a Gibbs sampling method based on the inside filteringoutside sampling algorithm [12], [36]. More specifically, thelatent tree t and the parameters θ, η, and λ are alternately
updated according to the conditional posterior distributionsp(t|θ, η, λ, z) and p(θ, η, λ|t, z), respectively.
In the inside filtering step, the conditional probability (insidemessage) that a subsequence zn:m is derived from a subtree
whose root node is A is denoted by
pA
n,m = p(zn:m |tn:m = A).

(8)This probability can be calculated recursively from the leaf
nodes to the root node as follows:pA
n,n = λA η(A → zn ),

(9)

pA
n,n+k
=

X
B,C∈V

(1−λA )θ(A → BC)XC
pB
n,n+l−1 pn+l,n+k . (10)1≤l≤k

In the outside sampling step, a latent tree t is obtained byrecursively branching paths from the start symbol S to the leaf
nodes according to p(t|θ, η, λ, z). Suppose that we alreadyhave a node tn:n+k = A. Two non-terminal symbols B andwhere u(A → BC) is the number of times that a binary
production rule A → BC is used, v(A → α) is the numberof times that an emission rule A → α is used, wA,0 is thenumber of times that a non-terminal symbol A branches into
two non-terminal symbols, and wA,1 is the number of times
that a non-terminal symbol A emits a chord.
2) Semi-Supervised Learning of PCFG: The PCFG model
can be trained in a semi-supervised manner by using treestructure annotations. Such annotations can be represented as
S-expressions on z (e.g., ((C, ((Dm, G), Am)))) that specify
the shape of a tree t. The tree shapes of t and tm:n without
information about non-terminal symbols are denoted by t̂ andt̂m:n , respectively.
To estimate the parameters θ, η, and λ from chord symbols z with a tree structure t̂, we extend the inside-filtering
algorithm. More specifically, we calculate the posterior distribution p(t|t̂, z, θ, η, λ) of a tree t having a fixed shape t̂ bymodifying Eq. (10) as follows:
X
pA(1−λA )θ(A → BC)q(B, C, n, k)(15)
n,n+k =B,C∈V

q(B, C, n, k)BCpn,n+l−1 pn+l,n+k if ∃l ∈ [1, k]=s.t. t̂n:n+l−1 ⇒ zn:n+l−1 , t̂n+l:n+k ⇒ zn+l:n+k ,
0 otherwise,
where t̂m:n ⇒ zm:n means that a subtree t̂m:n derives chordsymbols zm:n . Because of the constraints of the probabilities in
the above expression, the trees sampled in the outside sampling
step match the tree shape t̂.
3) HMM-Based Initialization of PCFG: In the unsupervised learning of the PCFG model (Section III-C1), nonterminal symbols given to the leaf nodes of a tree t tend tomismatch the musically meaningful harmonic functions. To
solve this, we propose a technique for initializing the PCFG
model by using an HMM. This technique is motivated by thefact that the HMM is able to automatically learn syntactic
categories corresponding to the harmonic functions [22].
First, an HMM that represents chords and syntactic functions as observed and latent variables, respectively, is trainedby using Gibbs sampling (see [22] for details). The most likely
sequence of latent variables {yn }Nn=1 is then estimated byusing the Viterbi algorithm. These latent variables are used to
specify the non-terminal symbols of the leaf nodes {tn:n }Nn=1

6IEEE/ACM TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING, VOL. XX, NO. XX, XX 2020of a latent tree. Specifically, in the first inside-filtering step of
Gibbs sampling for the PCFG, we use the following formula
instead of Eq. (9) to initialize the leaf nodes:(
λA (if A = yn ),A(16)pn,n =
0 (otherwise).
The probability parameters are randomly initialized. This
technique can be used with the semi-supervised learning of
the PCFG in Section III-C2.
4) Bayesian Learning of Metrical Markov Model: Given
chord onset positions y, the posterior distribution of π can beanalytically calculated as follows:
πa | y ∼ Dir(βa + qa ),(17)

where qa,b is the number of times that a chord has an onset
at the a-th position in a measure and has a duration of b timeunits in the training data.
5) Bayesian Learning of Pitch Markov Model: Given chordsymbols z, onset positions y, and melody pitches x, theposterior distribution of τ can be analytically calculated asfollows:
τac | z, y, x ∼ Dir(γac + rca ),

(18)

c

where r (a, b) is the number of times that pitch a transits topitch b in the region of chord symbol c in the training data.

where p(s) = p(t, z, y, x, o|Θ) is the complete joint probability of s based on the proposed model and p̄(s∗ |s) is a
proposal distribution. If the proposal is rejected, t, z, andy are not updated. In the global update, a proposed chord
sequence is always accepted. This process is iterated until
the posterior probability reaches a plateau. In the historyof generated samples, we use a sample that maximizes theposterior distribution p(t, z, y|x, o, Θ) as the final estimate ofchords.
B. Global Update
We propose two methods for jointly updating chord symbols
z and a latent tree t while fixing the chord onset positions y.
One is to use Gibbs sampling for stochastically drawing z andt from the posterior distribution p(z, t|y, x, o, Θ). The otheris to use a Viterbi algorithm for deterministically estimating z
and t that maximize p(z, t|y, x, o, Θ).
The stochastic method is similar to the inside filteringoutside sampling algorithm described in Section III-C1 exceptthat both t and z are estimated. The inside messages arecalculated recursively from the terminal symbols z to the startsymbol S as follows:X
pA
η(A → zn ) p(xn |zn ),(20)n,n = λAzn ∈ΣX

IV. M ELODY H ARMONIZATIONC
pA
θ(A → BC)pBn,n+k = (1 − λA )n,n+l−1 pn+l,n+k , (21)
B,C∈V
1≤l≤k

This section explains how to generate chord sequences fora given melody based on the proposed model.

where p(xn |zn ) is the probability that pitches xn are generated
from chord zn and is given by

A. Problem Specification
Given a melody with pitches x and onset positions o andmodel parameters Θ = {θ, η, λ, π, τ }, our goal is to estimate
an appropriate number of chords with symbols z and onset
positions y, where a latent tree t should also be estimated.
In this paper, we calculate the expectations of Θ under theposterior distributions given by Eqs. (12), (13), (14), (17), and(18) and use them as the point estimates of Θ.
To draw samples of t, z, and y from the posterior distribution p(t, z, y|x, o, Θ), we propose a Metropolis-Hastings(MH) sampler with four types of proposals:• Global update: Update z and t by using a stochastic or
deterministic method while keeping the number of chords
and y unchanged.
• Chord split: Update t, z, and y by choosing one of thechords and split it into two adjacent chords.
• Chord merge: Update t, z, and y by choosing two
adjacent chords derived from a node of t and merge theminto a single chord.
• Rhythm update: Update y by choosing a chord n andmove its onset position yn back or forth.
One of these operations is randomly selected and a new
sample s∗ = (t∗ , z∗ , y∗ ) is proposed by referring to a current
sample s = (t, z, y). The acceptance ratio of s∗ is given by
p(s∗ )p̄(s|s∗ )∗,(19)g(s , s) = min 1,
p(s)p̄(s∗ |s)p(xn |zn ) =InYτ zn (xn,i−1 , xn,i ),

(22)

i=1

where xn,0 is interpreted as xn−1,In−1 . A latent tree t is
obtained by recursively sampling paths from the start symbol
S to the leaf nodes as in the outside-sampling algorithm. Each
chord symbol zn is then sampled as follows:p(zn ) ∝ η(tn:n → z) p(xn |zn ).

(23)

In the deterministic method, the sum operations in Eqs. (20)and (21) are replaced with max operations. Instead of using
the outside sampling, a latent tree t is obtained by recursively
back-tracking the most likely paths from the start symbol S
to the chord symbols.
C. Chord Split and Merge
A chord is split or adjacent chords are merged by considering both a latent tree t and a melody. Tree t is locally updated
by these operations (Fig. 5) and the split operation has aninverse relationship with the merge operation.
In the split operation, a sample s∗ is given by splitting arandomly selected chord zn into z L and z R , estimating a new
onset position y ∗ ∈ [yn + 1, yn+1 − 1], and splitting a nonterminal symbol tn:n into two non-terminal symbols tL andtR . The proposal distribution is thus given by(θ(tn:n →tL tR )η(tL →z L )η(tR →z R )yn+1 ≥ yn + 1,∗
N (yn+1 −yn −1)p̄(s |s) =0otherwise.
(24)

H. TSUSHIMA et al.: BAYESIAN MELODY HARMONIZATION

 

 

 

 



Split a chord











7

V. E VALUATIONThis section reports three objective experiments for evaluating the unified model of chords and melody notes and anotherexperiment for evaluating the melody harmonization method
based on the unified model.

Merge chords




 

+1:+1









 





















Fig. 5. Chord split and merge operations.

In the merge operation, a sample s is given by merging arandomly selected pair of adjacent chords z L and z R into z
according to the probability η(tn:n → z). Thus we havep̄(s|s∗ ) =η(tn:n → z),
#MergeableNodes(s∗ )(25)where #MergeableNodes(s∗ ) is the number of pairs of adjacentchords that can be merged into s∗ , i.e., those chords forminga subtree with two leaves. The inverse proposal distribution
for the split operation is the same as the proposal distribution
for the merge operation and vice versa. In addition, we havep(s∗ )= (1 − λtn:n )θ(tn:n → tL tR )p(s)λ L η(tL → z L )λtR η(tR → z R )· t
λtn:n η(tn:n → zn )p(xL |z L )p(xR |z R )p(y ∗ |yn )p(yn+1 |y ∗ )·,p(xn |zn )p(yn+1 |yn )(26)

R
where xL
n and xn are sequences of pitches obtained by
splitting xn at a position y ∗ .
In the split operation, the acceptance ratio of s∗ given by
Eq. (19) is calculated by using Eqs. (24), (25), and (26). In themerge operation, the acceptance ratio of s∗ given by Eq. (19)is calculated by exchanging s and s∗ in Eqs. (24), (25), and(26). Through the split and merge operations, the number ofchords N is optimized stochastically.

D. Rhythm Update
Chord onset positions y are sampled from the conditional
posterior distribution p(y|t, z, x, o, Θ). A new sample s∗ isgiven by moving a randomly selected onset position yn to a
position yn∗ ∈ [yn−1 + 1, yn+1 − 1]. The proposal distributionp̄(s∗ |s) and its inverse version p̄(s|s∗ ) are given byp̄(s∗ |s) = p̄(s|s∗ ) =1
1.
N − 1 yn+1 − yn−1 − 1(27)

A. Experimental Conditions
To train the PCFG model, 399 chord sequences corresponding to musical sections of seven or eight measures (e.g., verse,
bridge, and chorus sections) were extracted from the Billboard
data [37] using its section annotations. Due to the limited
amount of available data, we limit the vocabulary of chordsymbols used for the experiments to the set of major andminor triads. The vocabulary of chord symbols consists ofthe combinations of the 12 root notes {C, C#, ..., B} andthe two chord types {major, minor}, and other chord types
are reduced to triads according to the root, third, and fifthnotes (diminished and augmented triads are discarded in theanalysis). We manually made tree structure annotations data2 .
In Eqs. (6) and (7), the elements of hyperparameter ν were
all set to 1.0 and those of ξ, ζ, β, and γ were all set to 0.1.
To train the parameters of an HMM used for initializing
the leaf nodes of a tree t, a set (called J-pop data) of chord
sequences of 4,872 Japanese popular songs obtained froma public web page3 was used. More specifically, an HMM
with four latent states was trained such that these statescorrespond to the three harmonic functions (tonic, dominant,subdominant) and an extra function as in [22]. Fig. 6 showsthe emission probabilities of the HMM. To train the metrical
Markov model described in Section III-B2 and the Markov
model described in Section III-B3, 9,902 pairs of melodies andchord sequences were extracted from 194 pieces of popular
music included in the Rock Corpus [38].
As shown in Fig. 4, we tested the following four learning
configurations of the PCFG model:(i) US-PCFG: Unsupervised learning (Section III-C1).
(ii) SS-PCFG: Semi-supervised learning based on tree
structure annotations (Section III-C2).
(iii) US-HMM-PCFG: Unsupervised learning with HMMbased initialization (Sections III-C1 and III-C3).
(iv) SS-HMM-PCFG: Semi-supervised learning based ontree structure annotations with HMM-based initialization(Sections III-C2 and III-C3).
All data (Billboard data, J-pop data, and Rock Corpus) weretransposed to the C major or C minor key (we used the
local key information when it is given, to deal with possible
modulations within each piece). To put emphasis on theMarkov model of a melody, each element of the trained τwas squared and then normalized.

The likelihood ratio is given byp(x∗n−1 |zn−1 )p(x∗n |zn )p(yn∗ |yn−1 )p(yn+1 |yn∗ )p(s∗ )=, (28)p(s)p(xn−1 |zn−1 )p(xn |zn )p(yn |yn−1 )p(yn+1 |yn )where x∗n−1 and x∗n are the sequences of pitches in the regionsof chords n−1 and n with the new onset position yn∗ . The
acceptance ratio of s∗ given by Eq. (19) is calculated by using
Eqs. (27) and (28).

B. Statistical Chord ModelingWe investigated the performance of the PCFG model inpredicting chord sequences, for varying numbers of nonterminal symbols K. The number K was changed from 1
2 The data is available at http://sap.ist.i.kyoto-u.ac.jp/members/tsushima/chord tree annotation/
3 J-Total Music: http://music.j-total.net

8IEEE/ACM TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING, VOL. XX, NO. XX, XX 2020

4
(Others)prob

3(Dominant)prob

2(Subdominant)

0.6
0.5
0.4
0.3
0.2
0.1
0

prob

1
(Tonic)

0.6
0.5
0.4
0.3
0.2
0.1
0

0.6
0.5
0.4
0.3
0.2
0.1
0

prob

State ID
(Interpretation)

0.6
0.5
0.4
0.3
0.2
0.1
0

C

Am

F

A

E♭

B♭

Em

F#m Cm

F#

F

Dm

D

B♭

A♭

Gm

Fm

Bm

B

C

G

C

Em

F

E

Am

Dm

B♭

Fm

D

A♭

Cm

C#

B♭

E♭

F#

Fm B♭m Gm E♭mFig. 6. The emission probabilities of an HMM. Ten chord symbols withhighest emission probabilities are shown for each state.

to 25 in the framework of five-fold cross validation. As anevaluation measure, we calculated the perplexity for an unseen
chord sequence z as follows:!
p(z|λ, θ, η)
1Pchord = exp − ln P
,
(29)
0
N
z0 :|z0 |=N p(z |λ, θ, η)
where p(z|λ, θ, η) = pS1,N is the marginal probability obtainedby using the inside algorithm given by Eqs. (9) and (10).
As shown in Fig. 7, US-HMM-PCFG achieved the lowest(best) perplexities and outperformed US-PCFG for almost all
K. This indicates that the HMM-based initialization method
improves the generalization capability of the PCFG model for
unseen chords. Interestingly, SS-HMM-PCFG underperformed
US-HMM-PCFG. This means that the semi-supervised learning method using tree structure annotations does not contributeto chord prediction. This is presumably because there are many
possible structural interpretations for chord sequences and theunsupervised learning method is better at dealing with such
ambiguity in a probabilistic manner. Nonetheless, SS-HMMPCFG can be considered useful for uniquely estimating themost likely tree behind chords (Section V-D).
The improvement of the test-set perplexity was saturated
before K reached 20. Considering the computational cost, theoptimal values of K were estimated as 12, 17, 16, and 18 forlearning methods (i)–(iv), respectively, which were used forthe following evaluations.

C. Statistical Melody ModelingWe investigated the performance of the unified model inpredicting melodies by using 172 melodies extracted from the
RWC music database [39]. These melodies were disjoint withthe datasets (Billboard data, J-pop data, and Rock Corpus)
used for training the unified model. In this experiment, we
assumed that the chord onset positions y were fixed to barlines, to enable analytically marginalizing out chord symbols
and the tree in the calculation of the melody perplexity defined

below. As an evaluation measure, we calculated the perplexity
for an unseen melody (x, o) (per note) as follows:!
p(x, o|y, λ, θ, η)
1,(30)
Pmelody = exp − ln P
I
z:|z|=N p(z|λ, θ, η)
where p(x, o|y, λ, θ, η) = pS1,N is the marginal probability
obtained by using the inside algorithm given by Eqs. (20)and (21). Note that chord symbols z and a tree t are bothmarginalized out.
For comparison, we tested a measure-wise HMM that
represents chords and melody notes as latent and observed
variables, respectively. To conduct a fair comparison with the
PCFG model, chords were allowed to change only at barlines in this model. As a generative process, one chord foreach measure is generated according to a Markov model, andmelody notes in each measure is generated by a Markov model
conditionally dependent on the generated chord as in Eq. (4).
The transition probabilities of chords were learned from theBillboard data. Given a melody, the most likely chord sequence
was estimated using the Viterbi algorithm.
As shown in Fig. 8, the results of melody modeling wereconsistent with those of chord modeling shown in Section V-B.
US-HMM-PCFG attained the lowest perplexity and US-PCFG,
US-HMM-PCFG, and SS-HMM-PCFG outperformed the conventional HMM (6.09). The HMM-based initialization of the
PCFG model improved the ability of melody prediction.
D. Tree-Structured Parsing
We investigated the performance of the unified model intree-structured parsing of chord sequences. 5-fold cross validation on the Billboard dataset [37] was conducted and theunsupervised learning method (Section III-C1) was comparedwith the semi-supervised learning method (Section III-C2). Asan evaluation measure, we calculated the tree edit distance(TED) [40] between an estimated tree and the ground-truth
tree. Since the non-terminal symbols of the ground-truth tree
are not given, only the estimated tree shape was evaluated.
As shown in Fig. 9, the semi-supervised learning method
turned out to be effective for improving the performance oftree-structured parsing. For the unsupervised learning methods, the HMM-based initialization also led to better treestructured parsing results. Interestingly, applying both thesemi-supervised learning and the HMM-based initialization of
the PCFG model slightly degraded the performance. Whenthe leaf nodes of a tree were guided to represent harmonic
functions, the production rules used in a higher level were
often disjointed from those used in a lower level. This might
prevent the model from inducing a globally consistent tree
from a chord sequence. Overall, we confirmed the efficacy ofthe proposed learning techniques for improving the similarity
between human structural annotations of chord sequences andthose generated by the PCFG model.
E. Melody Harmonization
We investigated the performance of the proposed melody
harmonization method based on the unified probabilistic modelH. TSUSHIMA et al.: BAYESIAN MELODY HARMONIZATION

US-PCFG
US-HMM-PCFG

10

Test-data perplexity per chord

9

10

9

9

8

8

7

7

6

6

5
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25

Number of distinct non-terminal symbols K(a)5

SS-PCFG
SS-HMM-PCFG

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25

Number of distinct non-terminal symbols K(b)Fig. 7. Test-set perplexities per chord symbol obtained by the four different learning configurations of the PCFG model: (a) US-PCFG and US-HMM-PCFG,(b) SS-PCFG and SS-HMM-PCFG. The band and boxes indicate the median and interquartile range (IQR), and the whiskers indicate the 1.5 times wider
range as IQR. Points outside this range are identified as outliers.

Test-data perplexity per note

6.14
6.12
6.1
6.093

Measure-wise HMM

6.08
6.06
6.04
6.02

20 21 22 23 24 25

6

US-PCFG

SS-PCFG

US-HMM-PCFG SS-HMM-PCFG

Fig. 8. Test-set perplexities per melody note. Indicators in the box plots are
the same as those in Fig. 7.
9Tree edit distance

8
7
6
5
4
3

20 21 22 23 24 25

2

US-PCFG

SS-PCFG

US-HMM-PCFG SS-HMM-PCFG

Fig. 9. Tree edit distances between estimated and ground-truth trees.
Indicators in the box plots are the same as those in Fig. 7.

of a chord sequence and a melody. For evaluation, we extracted
172 pairs of human-composed melodies and chord sequencesfrom the RWC music database [39]. We tested eight configurations by combining the four learning configurations (i)–(iv)with the sampling method and the Viterbi algorithm describedin Section IV-B. For comparison, we tested a standard HMMbased method, similar to the model of [7], that represents chordtransitions at the 16th-note level. We also tested a rule-based
method implemented in the Melisma Music Analyzer [41] thatgenerates a sequence of chord root notes for a given melody.
For the proposed method, we ran the MCMC sampler togenerate a sufficient number (= 5000) of chord sequences
for a given melody. From the generated chord sequences,
one with the highest posterior probability was chosen as the
final chord sequence. We also used a set of five samples ofchord sequences with the highest posterior probabilities in
the calculation of the mean reciprocal rank (MRR) explainedlater. Since the learning methods of the unified model (i)–(iv) and the proposed melody harmonization method are allbased on MCMC sampling, we took the average score of eachmeasure over 500 trials, i.e., we ran the harmonization method
five times using 100 different parameter sets of the PCFG
model. For the HMM-based method, we similarly sampled
a sufficient number (= 1000) of chord sequences using theforward filtering-backward sampling algorithm [43] and chose
five samples with the highest posterior probabilities for the
calculation of the MRR. For the other metrics, we used
the most probable chord sequence obtained by the Viterbi
algorithm.
The methods were evaluated by comparing the generated
chord sequences with the original human-composed chord
sequence. To evaluate the local correctness of generated
chord sequences, we calculated the accuracy by comparing
the generated sequence with the highest posterior probability
with the original sequence at the 16th-note level. We also
calculated the tonal pitch step distance [42], which measuresthe dissimilarity between two chord sequences based on a
cognitive model of tonality (see [42] for details). We also
calculated the mean reciprocal rank (MRR) by comparing thefive chord sequences with the original sequence at the 16thnote level. To evaluate the global correctness of generated
chord sequences, we calculated the Levenshtein edit distance
between the generated sequence and the original sequence. To
compare the proposed method with the rule-based method, we
also calculated the Levenshtein distance between the generated
and original sequences of root notes by ignoring the chordtypes.
As shown in Fig. 10, the Viterbi algorithm consistently
outperformed the sampling method and US-HMM-PCFG per-

10

IEEE/ACM TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING, VOL. XX, NO. XX, XX 2020(%)(a) Accuracy(b) Tonal pitch step distance
4

0.27

23

3.8

0.26

22

3.6

0.25

21

3.4

0.24

20

3.2

0.23

19

3

0.22

24

(c) MRR(a) Accuracy

18
USPCFGSSPCFG

USHMMPCFG

SSHMMPCFG

0.21

2.8
HMMbased

Rulebased

USPCFG(d) Levenshtein distance (chord root note)

SSPCFG

USHMMPCFG

SSHMMPCFG

HMMbased

Rulebased

USPCFGSSPCFG

USHMMPCFG

SSHMMPCFG

HMMbased(e) Levenshtein distance (chord symbol)
99

96

98
94

97
96

92

95
94

90

93
92

88

91
86
USPCFGSSPCFG

USHMMPCFG

SSHMMPCFG

based

90
HMMbased

Rulebased

USPCFGSSPCFG

USHMMPCFG: Viterbi

SSHMMPCFG

HMMbased: Sampling

Fig. 10. Evaluation results of melody harmonization: (a) accuracy, (b) tonal pitch step distance [42], (c) MRR, (d) Levenshtein distance for chord root notes,and (e) Levenshtein distance for chord symbols. For accuracy and MRR, the higher, the better. For Levenshtein distance, the lower, the better. We plot thepiece-wise mean and standard error corresponding to a 68% confidence interval.

formed best in terms of the MRR and the Levenshtein distance.
In terms of accuracy, the proposed method outperformed the
HMM-based method (19.6%), but underperformed the rulebased method (22.1%). In terms of the MRR, the proposed
method outperformed the HMM-based method (0.233). As
for the tonal pitch step distance and Levenshtein distance forroot notes, the proposed method with the Viterbi algorithm
outperformed both the HMM-based method and the rule-based
method. In the case of the Levenshtein distance for chordsymbols, the proposed method significantly outperformed the
HMM-based method (97.2). These results show that the proposed method was inferior to the rule-based method in termsof the accuracy, but outperformed it in terms of the othercompared metrics. A possible reason is that while the PCFG
is able to capture the global structure with a derivation tree, itcannot explicitly model the sequential dependency of chords,
unlike Markov models.
F. Example Results
Fig. 11 shows examples of chord sequences generated by the
proposed and conventional methods. The numbers of distinct
non-terminal symbols of US-PCFG and SS-HMM-PCFG were
set to 12 and 17, respectively. The HMM-based method often
generated a sequence consisting of a very few chord symbols.
The rule-based method generated a chord sequence consistingof a reasonable number of symbols but could not generatemusically natural rhythms.
In contrast, the proposed method (US-PCFG and SS-HMMPCFG) were able to successfully generate a chord sequencewith a reasonable number of symbols with natural rhythms.
This demonstrates the effectiveness of the metrical Markov

model of chord onset positions. In US-PCFG, leaf nodes above
different chord symbols were given different labels, whereas
in SS-HMM-PCFG, the leaf nodes above C-major chord andA-minor chord, which are usually categorized as tonic chords,
were given the same label S1 . In addition, in the latter case,the generated chord sequence had cadences and chord rhythms
that were similar to those of the ground truth. We observedthat the two models had a similar tendency for other melodies
as well (see online supplemental material4 ). This indicatesthat the SS-HMM-PCFG can successfully reflect syntactic
functions and hierarchical structure of chords in the generated
results.
G. Open Problems
The proposed melody harmonization method tends to generate conservative chord symbols (e.g., C major and G major)
that frequently appear in the training data because they tend
to maximize the posterior probabilities of chords for given
melodies. To increase the diversity and musical attractivenessof the generated chord sequences, it would be effective toincorporate prior knowledge about chord patterns (idioms andphrases) that frequently appear in actual musical pieces into
our MCMC sampler. Another solution is to replace the PCFG
model that splits a non-terminal node into two child nodes
with a tree-substitution grammar (TSG) model that splitsa non-terminal node into two variable-depth subtrees. This
enables us to automatically discover tree-structured patterns
frequently found in chord sequences. We also plan to replace
the generative model of a melody (Markov model) with a deep
4 https://pcfgharmonization.github.io/

H. TSUSHIMA et al.: BAYESIAN MELODY HARMONIZATION

Ground truth :

C

Am

G

HMM-based :11

CRule-based :

G

F

C

C

G

G

G

G

F

G

Dm

C

(only root notes)

Dm

G

S
S7
S5

US-PCFG :

S6

S1

S4

S2

S3

S1

C

Am

F

G

C

S
S6

S6

S4
SS-HMMPCFG :

S4

S5

S5

S1

S3

S1

S1

S2

S3

S2

S3

C

G

Am

C

Dm

G

Dm

G

Fig. 11. Examples of generated chord sequences. From top to bottom, the input melody, the human-composed chord sequence, and chord sequences generated
by the HMM-based method, the rule-based method, US-PCFG, and SS-HMM-PCFG are shown.

generative model (e.g., recurrent neural network), which iseffective for learning complex dependency in time series.
While our method can capture syntactic functions and the
local structure behind chords, it is still weak in capturing
global structure such as phrase structure. A possible solution
is to use the combinatory categorial grammar (CCG) [44] formodeling chord sequences. A CCG can explicitly describe theglobal syntactic structure (e.g., phrase, clause) by using symbols constructed from the combination of simple categories
of words (e.g., “noun”, “verb”). By inferring the CCG fromthe data of chord sequences, it would be possible to directly
analyze the phrase structure in chord sequences. It is also
important to improve the naturalness of chord rhythms. Since
chord rhythms depend on rhythms of melody notes, we plan
to formulate a model of melody note rhythms and integrate itwith our unified model.

accuracy of generated chord sequences. We also confirmed thatthe semi-supervised learning of the PCFG model is superior tothe completely unsupervised learning, in terms of the ability to
predict tree structures behind chords, and that the HMM-based
initialization of the PCFG model improved the performance of
melody harmonization.
Our next direction is to construct an interactive music composition/arrangement system. Interactive systems are usefulfor those people who want to incorporate their preference inthe generated music and can go beyond the reach of fully
automated systems. The ability of the proposed model to learn
and represent the grammatical structure of music in a similar
way as humans do is essential for enhancing the directabilityof the system, as partly shown in [13]. Through large-scale
user studies we plan to reveal human’s music creation processfrom the computational perspective.

ACKNOWLEDGEMENT
We are grateful to Ryo Nishikimi for helping the computer
program implementation of our method.
VI. C ONCLUSION
This paper has presented a Bayesian melody harmonization
method that aims to generate a chord sequence with natural
rhythms for a given melody5 . This method is based on aunified generative model of a chord sequence and a melody
involving a metrical Markov model describing chord rhythms
and a PCFG model describing the tree structure behind chords.
Experimental results showed that the proposed method outperformed the conventional HMM-based method in terms of
5 The source code of our melody harmonization method is available in the
online supplementary page (https://pcfgharmonization.github.io/).

R EFERENCES[1] K. Ebcioğlu, “An expert system for harmonizing four-part chorales,”Computer Music Journal, vol. 12, no. 3, pp. 43–51, 1988.
[2] G. Papadopoulos and G. Wiggins, “AI methods for algorithmic composition: A survey, a critical view and future prospects,” in AISB Symposium
on Musical Creativity, 1999, pp. 110–117.
[3] M. Towsey, A. Brown, S. Wright, and J. Diederich, “Towards melodic
extension using genetic algorithms,” Educational Technology & Society,vol. 4, no. 2, pp. 54–65, 2001.
[4] M. Allan and C. Williams, “Harmonising chorales by probabilistic
inference,” in NIPS, 2005, pp. 25–32.
[5] J. F. Paiement, D. Eck, and S. Bengio, “Probabilistic melodic harmonization,” in CSCSI, 2006, pp. 218–229.
[6] C. H. Chuan and E. Chew, “A hybrid system for automatic generationof style-specific accompaniment,” in IJWCC, 2007, pp. 57–64.
[7] I. Simon, D. Morris, and S. Basu, “Mysong: automatic accompaniment
generation for vocal melodies,” in SIGCHI Conference on Human
Factors in Computing Systems. ACM, 2008, pp. 725–734.

12IEEE/ACM TRANSACTIONS ON AUDIO, SPEECH, AND LANGUAGE PROCESSING, VOL. XX, NO. XX, XX 2020

[8] R. D. Prisco and R. Zaccagnino, “An evolutionary music composer
algorithm for bass harmonization,” in Applications of Evolutionary
Computing. Springer, 2009, pp. 567–572.
[9] R. D. Prisco, A. Eletto, A. Torre, and R. Zaccagnino, “A neural network
for bass functional harmonization,” in European Conference on the
Applications of Evolutionary Computation. Springer, 2010, pp. 351–
360.
[10] R. Groves, “Automatic harmonization using a hidden semi-Markov
model,” in AIIDE, 2013, pp. 48–54.
[11] S. A. Raczyński, S. Fukayama, and E. Vincent, “Melody harmonization
with interpolated probabilistic models,” Journal of New Music Research,
vol. 42, no. 3, pp. 223–235, 2013.
[12] H. Tsushima, E. Nakamura, K. Itoyama, and K. Yoshii, “Function- andrhythm-aware melody harmonization based on tree-structured parsing
and split-merge sampling of chord sequences,” in ISMIR, 2017, pp. 502–
508.
[13] ——, “Interative arrangement of chords and melodies based on a treestructured generative model,” in ISMIR, 2017, pp. 502–508.
[14] G. Hadjeres and F. Pachet, “DeepBach: A steerable model for Bachchorales generation,” in ICML, 2017, pp. 1362–1371.
[15] A. Cadwallader and D. Gagné, Analysis of Tonal Music: A SchenkerianApproach (3rd ed.). Oxford University Press, 2011.
[16] H. Riemann, Harmony Simplified: Or the Theory of the Tonal Functionsof Chords (2nd ed.). Augener, 1986.
[17] F. Lerdahl and R. Jackendoff, A Generative Theory of Tonal Music.
MIT press, 1985.
[18] M. J. Steedman, “A generative grammar for jazz chord sequence,” Music
Perception, vol. 2, no. 1, pp. 52–77, 1984.
[19] M. Rohrmeier, “Mathematical and computational approaches to music
theory, analysis, composition and performance,” Journal of Mathematics
and Music, vol. 5, no. 1, pp. 35–53, 2011.
[20] W. B. De Haas, J. P. Magalhães, F. Wiering, and R. C. Veltkamp,“Automatic functional harmonic analysis,” Computer Music Journal,
vol. 37, no. 4, pp. 37–53, 2013.
[21] D. Quick, “Learning production probabilities for musical grammars,”
Journal of New Music Research, vol. 45, no. 4, pp. 295–313, 2016.
[22] H. Tsushima, E. Nakamura, K. Itoyama, and K. Yoshii, “Generative
statistical models with self-emergent grammar of chord sequences,”
Journal of New Music Research, 2018.
[23] C. Raphael, “A hybrid graphical model for rhythmic parsing,” Artificial
Intelligence, vol. 137, no. 1, pp. 217–238, 2002.
[24] M. Hamanaka, M. Goto, H. Asoh, and N. Otsu, “A learning-based
quantization: Unsupervised estimation of the model parameters,” in
ICMC, 2003, pp. 369–372.
[25] H. Lim, S. Rhyu, and K. Lee, “Chord generation from symbolic melody
using BLSTM networks,” in Proc. ISMIR, 2017, pp. 621–627.
[26] C. Ames, “The Markov process as a compositional model: a survey andtutorial,” Leonardo, vol. 22, no. 2, pp. 175–187, 1989.
[27] F. Pachet, R. Roy, and G. Barbieri, “Finite-length markov processes
with constraints,” in Twenty-Second International Joint Conference on
Artificial Intelligence, 2011.
[28] M. Hamanaka, K. Hirata, and S. Tojo, “Implementing ‘A generative
theory of tonal music’,” Journal of New Music Research, vol. 35, no. 4,
pp. 249–277, 2006.
[29] ——, “Musical structural analysis database based on GTTM,” in ISMIR,
2014, pp. 325–330.
[30] E. Nakamura, M. Hamanaka, K. Hirata, and K. Yoshii, “Tree-structured
probabilistic model of monophonic written music based on the generative theory of tonal music,” in IEEE ICASSP, 2016, pp. 276–280.
[31] R. Scholz, E. Vincent, and F. Bimbot, “Robust modeling of musicalchord sequences using probablistic N-grams,” in IEEE ICASSP, 2009,
pp. 53–56.
[32] K. Yoshii and M. Goto, “A vocabulary-free infinity-gram model for
nonparametric bayesian chord progression analysis.” in ISMIR, 2011,
pp. 645–650.
[33] J. F. Paiement, D. Eck, and S. Bengio, “A probabilistic model for chord
progressions,” in ISMIR, 2005, pp. 312–319.
[34] S. Fukayama, K. Yoshii, and M. Goto, “Chord-sequence-factory: A
chord arrangement system modifying factorized chord sequence probabilities,” in ISMIR, 2013, pp. 457–462.
[35] S. F. Chen and J. Goodman, “An empirical study of smoothing techniques for language modeling,” in ACL, 1996, pp. 310–318.
[36] M. Johnson, T. L. Griffiths, and S. Goldwater, “Bayesian inference for
PCFGs via Markov chain Monte Carlo,” in NAACL-HLT, 2007, pp. 139–
146.

[37] J. B. L. Smith, J. A. Burgoyne, I. Fujinaga, D. D. Roure, and J. S.
Downie, “Design and creation of a large-scale database of structural
annotations,” in ISMIR, 2011, pp. 555–560.
[38] T. D. Clercq and D. Temperley, “A corpus analysis of rock harmony,”
Popular Music, vol. 30, no. 01, pp. 47–70, 2011.
[39] M. Goto, H. Hashiguchi, T. Nishimura, and R. Oka, “RWC music
database: Popular, classical and jazz music databases.” in ISMIR, 2002,
pp. 287–288.
[40] K. Zhang and D. Shasha, “Simple fast algorithms for the editing distance
between trees and related problems,” SIAM journal on computing,vol. 18, no. 6, pp. 1245–1262, 1989.
[41] D. Temperley and D. Sleator, “The Melisma Music Analyzer,” 2001,
http://www.link.cs.cmu.edu/melisma/.
[42] W. B. de Haas, F. Wiering, and R. C. Veltkamp, “A geometricaldistance measure for determining the similarity of musical harmony,”International Journal of Multimedia Information Retrieval, vol. 2, no. 3,
pp. 189–202, 2013.
[43] J. Van Gael, Y. Saatci, Y. W. Teh, and Z. Ghahramani, “Beam sampling
for the infinite hidden Markov model,” in Proc. ICMC, 2008, pp. 1088–1095.
[44] M. Steedman and J. Baldridge, “Combinatory categorial grammar,” NonTransformational Syntax: Formal and Explicit Models of Grammar.
Wiley-Blackwell, 2011.

Hiroaki Tsushima He received a masters degree in Informatics from Kyoto
University in 2019. He was a member of Speech and Audio Processing Group,
Graduate School of Informatics, Kyoto University. His research interests
include automatic music composition and machine learning.

Eita Nakamura Eita Nakamura He received a Ph.D. degree in physics fromthe University of Tokyo in 2012. He has been a post-doctoral researcher at
the National Institute of Informatics, Meiji University, and Kyoto University.
He is currently an Assistant Professor at the Hakubi Center for Advanced
Research and Graduate School of Informatics, Kyoto University. His researchinterests include music modelling and analysis, music information processing
and statistical machine learning.

Kazuyoshi Yoshii He received the Ph.D. degree in
informatics from Kyoto University, Japan, in 2008.
He is currently a Senior Lecturer at Kyoto University. His research interests include music signal
processing and machine learning. He is a Memberof the Information Processing Society of Japan andInstitute of Electronics, Information, and Communication Engineers.