#!/usr/bin/env python
# -*- encoding: utf-8 -*-
"""Functions for interval construction"""
from importlib import resources
from typing import Collection, Dict, List, Union, overload, Iterable
from typing_extensions import Literal
import msgpack
import numpy as np
from numpy.typing import ArrayLike
from .._cache import cache
from .._typing import _FloatLike_co
with resources.path("librosa.core", "intervals.msgpack") as imsgpack:
with imsgpack.open("rb") as _fdesc:
# We use floats for dictionary keys, so strict mapping is disabled
INTERVALS = msgpack.load(_fdesc, strict_map_key=False)
[docs]@cache(level=10)
def interval_frequencies(
n_bins: int,
*,
fmin: _FloatLike_co,
intervals: Union[str, Collection[float]],
bins_per_octave: int = 12,
tuning: float = 0.0,
sort: bool = True
) -> np.ndarray:
"""Construct a set of frequencies from an interval set
Parameters
----------
n_bins : int
The number of frequencies to generate
fmin : float > 0
The minimum frequency
intervals : str or array of floats in [1, 2)
If `str`, must be one of the following:
- `'equal'` - equal temperament
- `'pythagorean'` - Pythagorean intervals
- `'ji3'` - 3-limit just intonation
- `'ji5'` - 5-limit just intonation
- `'ji7'` - 7-limit just intonation
Otherwise, an array of intervals in the range [1, 2) can be provided.
bins_per_octave : int > 0
If `intervals` is a string specification, how many bins to
generate per octave.
If `intervals` is an array, then this parameter is ignored.
tuning : float
Deviation from A440 tuning in fractional bins.
This is only used when `intervals == 'equal'`
sort : bool
Sort the intervals in ascending order.
Returns
-------
frequencies : array of float
The frequencies
Examples
--------
Generate two octaves of Pythagorean intervals starting at 55Hz
>>> librosa.interval_frequencies(24, fmin=55, intervals="pythagorean", bins_per_octave=12)
array([ 55. , 58.733, 61.875, 66.075, 69.609, 74.334, 78.311,
82.5 , 88.099, 92.812, 99.112, 104.414, 110. , 117.466,
123.75 , 132.149, 139.219, 148.668, 156.621, 165. , 176.199,
185.625, 198.224, 208.828])
Generate two octaves of 5-limit intervals starting at 55Hz
>>> librosa.interval_frequencies(24, fmin=55, intervals="ji5", bins_per_octave=12)
array([ 55. , 58.667, 61.875, 66. , 68.75 , 73.333, 77.344,
82.5 , 88. , 91.667, 99. , 103.125, 110. , 117.333,
123.75 , 132. , 137.5 , 146.667, 154.687, 165. , 176. ,
183.333, 198. , 206.25 ])
Generate three octaves using only three intervals
>>> intervals = [1, 4/3, 3/2]
>>> librosa.interval_frequencies(9, fmin=55, intervals=intervals)
array([ 55. , 73.333, 82.5 , 110. , 146.667, 165. , 220. ,
293.333, 330. ])
"""
if isinstance(intervals, str):
if intervals == "equal":
# Maybe include tuning here?
ratios = 2.0 ** (
(tuning + np.arange(0, bins_per_octave, dtype=float)) / bins_per_octave
)
elif intervals == "pythagorean":
ratios = pythagorean_intervals(bins_per_octave=bins_per_octave, sort=sort)
elif intervals == "ji3":
ratios = plimit_intervals(
primes=[3], bins_per_octave=bins_per_octave, sort=sort
)
elif intervals == "ji5":
ratios = plimit_intervals(
primes=[3, 5], bins_per_octave=bins_per_octave, sort=sort
)
elif intervals == "ji7":
ratios = plimit_intervals(
primes=[3, 5, 7], bins_per_octave=bins_per_octave, sort=sort
)
else:
ratios = np.array(intervals)
bins_per_octave = len(ratios)
# We have one octave of ratios, tile it up to however many we need
# and trim back to the right number of bins
n_octaves = np.ceil(n_bins / bins_per_octave)
all_ratios = np.multiply.outer(2.0 ** np.arange(n_octaves), ratios).flatten()[
:n_bins
]
if sort:
all_ratios = np.sort(all_ratios)
return all_ratios * fmin
@overload
def pythagorean_intervals(
*,
bins_per_octave: int = ...,
sort: bool = ...,
return_factors: Literal[False] = ...
) -> np.ndarray:
...
@overload
def pythagorean_intervals(
*, bins_per_octave: int = ..., sort: bool = ..., return_factors: Literal[True]
) -> List[Dict[int, int]]:
...
@overload
def pythagorean_intervals(
*, bins_per_octave: int = ..., sort: bool = ..., return_factors: bool = ...
) -> Union[np.ndarray, List[Dict[int, int]]]:
...
[docs]@cache(level=10)
def pythagorean_intervals(
*, bins_per_octave: int = 12, sort: bool = True, return_factors: bool = False
) -> Union[np.ndarray, List[Dict[int, int]]]:
"""Pythagorean intervals
Intervals are constructed by stacking ratios of 3/2 (i.e.,
just perfect fifths) and folding down to a single octave::
1, 3/2, 9/8, 27/16, 81/64, ...
Note that this differs from 3-limit just intonation intervals
in that Pythagorean intervals only use positive powers of 3
(ascending fifths) while 3-limit intervals use both positive
and negative powers (descending fifths).
Parameters
----------
bins_per_octave : int
The number of intervals to generate
sort : bool
If `True` then intervals are returned in ascending order.
If `False`, then intervals are returned in circle-of-fifths order.
return_factors : bool
If `True` then return a list of dictionaries encoding the prime factorization
of each interval as `{2: p2, 3: p3}` (meaning `3**p3 * 2**p2`).
If `False` (default), return intervals as an array of floating point numbers.
Returns
-------
intervals : np.ndarray or list of dictionaries
The constructed interval set. All intervals are mapped
to the range [1, 2).
See Also
--------
plimit_intervals
Examples
--------
Generate the first 12 intervals
>>> librosa.pythagorean_intervals(bins_per_octave=12)
array([1. , 1.067871, 1.125 , 1.201355, 1.265625, 1.351524,
1.423828, 1.5 , 1.601807, 1.6875 , 1.802032, 1.898437])
>>> # Compare to the 12-tone equal temperament intervals:
>>> 2**(np.arange(12)/12)
array([1. , 1.059463, 1.122462, 1.189207, 1.259921, 1.33484 ,
1.414214, 1.498307, 1.587401, 1.681793, 1.781797, 1.887749])
Or the first 7, in circle-of-fifths order
>>> librosa.pythagorean_intervals(bins_per_octave=7, sort=False)
array([1. , 1.5 , 1.125 , 1.6875 , 1.265625, 1.898437,
1.423828])
Generate the first 7, in circle-of-fifths other and factored form
>>> librosa.pythagorean_intervals(bins_per_octave=7, sort=False, return_factors=True)
[
{2: 0, 3: 0},
{2: -1, 3: 1},
{2: -3, 3: 2},
{2: -4, 3: 3},
{2: -6, 3: 4},
{2: -7, 3: 5},
{2: -9, 3: 6}
]
"""
# Generate all powers of 3 in log space
pow3 = np.arange(bins_per_octave)
# Using modf here to quickly get the fractional part of the log,
# accounting for whatever power of 2 is necessary to get 3**k
# within the octave.
log_ratios: np.ndarray
pow2: np.ndarray
log_ratios, pow2 = np.modf(pow3 * np.log2(3))
# If the fractional part is negative, add
# one more power of two to get it into the range [0, 1).
too_small = log_ratios < 0
log_ratios[too_small] += 1
pow2[too_small] += 1
# Convert powers of 2 to integer
pow2 = pow2.astype(int)
idx: Iterable[int]
if sort:
# Order the intervals
idx = np.argsort(log_ratios)
log_ratios = log_ratios[idx]
else:
# If not sorting, we'll take powers in order
idx = range(bins_per_octave)
if return_factors:
return list({2: -pow2[i], 3: pow3[i]} for i in idx)
return np.power(2, log_ratios)
def __harmonic_distance(logs, a, b):
"""Compute the harmonic distance between ratios a and b.
Harmonic distance is defined as `log2(a * b) - 2*log2(gcd(a, b))` [#]_.
Here we are expressing a and b as prime factorization exponents,
and the prime basis are provided in their log2 form.
.. [#] Tenney, James.
"On ‘Crystal Growth’ in harmonic space (1993–1998)."
Contemporary Music Review 27.1 (2008): 47-56.
"""
a = np.array(a)
b = np.array(b)
# numerator = positive exponents
a_num = np.maximum(a, 0)
# denominator = negative exponents
a_den = a_num - a
b_num = np.maximum(b, 0)
b_den = b_num - b
# log2(ab / gcd(a,b)**2) = log(a) + log(b) - 2 * log(gcd(a,b))
# gcd(a,b) for rationals: gcd(a_num, b_num) / lcm(a_den, b_den)
# gcd = minimum(a_num, b_num) and lcm = maximum(a_den, b_den)
gcd = np.minimum(a_num, b_num) - np.maximum(a_den, b_den)
# Rounding this to 6 decimals to avoid floating point weirdness
return np.around(logs.dot(a + b - 2 * gcd), 6)
def _crystal_tie_break(a, b, logs):
"""Given two tuples of prime powers, break ties."""
return logs.dot(np.abs(a)) < logs.dot(np.abs(b))
@overload
def plimit_intervals(
*,
primes: ArrayLike,
bins_per_octave: int = ...,
sort: bool = ...,
return_factors: Literal[False] = ...
) -> np.ndarray:
...
@overload
def plimit_intervals(
*,
primes: ArrayLike,
bins_per_octave: int = ...,
sort: bool = ...,
return_factors: Literal[True]
) -> List[Dict[int, int]]:
...
@overload
def plimit_intervals(
*,
primes: ArrayLike,
bins_per_octave: int = ...,
sort: bool = ...,
return_factors: bool = ...
) -> Union[np.ndarray, List[Dict[int, int]]]:
...
[docs]@cache(level=10)
def plimit_intervals(
*,
primes: ArrayLike,
bins_per_octave: int = 12,
sort: bool = True,
return_factors: bool = False
) -> Union[np.ndarray, List[Dict[int, int]]]:
"""Construct p-limit intervals for a given set of prime factors.
This function is based on the "harmonic crystal growth" algorithm
of [#1]_ [#2]_.
.. [#1] Tenney, James.
"On ‘Crystal Growth’ in harmonic space (1993–1998)."
Contemporary Music Review 27.1 (2008): 47-56.
.. [#2] Sabat, Marc, and James Tenney.
"Three crystal growth algorithms in 23-limit constrained harmonic space."
Contemporary Music Review 27, no. 1 (2008): 57-78.
Parameters
----------
primes : array of odd primes
Which prime factors are to be used
bins_per_octave : int
The number of intervals to construct
sort : bool
If `True` then intervals are returned in ascending order.
If `False`, then intervals are returned in crystal growth order.
return_factors : bool
If `True` then return a list of dictionaries encoding the prime factorization
of each interval as `{2: p2, 3: p3, ...}` (meaning `3**p3 * 2**p2`).
If `False` (default), return intervals as an array of floating point numbers.
Returns
-------
intervals : np.ndarray or list of dictionaries
The constructed interval set. All intervals are mapped
to the range [1, 2).
See Also
--------
pythagorean_intervals
Examples
--------
Compare 3-limit tuning to Pythagorean tuning and 12-TET
>>> librosa.plimit_intervals(primes=[3], bins_per_octave=12)
array([1. , 1.05349794, 1.125 , 1.18518519, 1.265625 ,
1.33333333, 1.40466392, 1.5 , 1.58024691, 1.6875 ,
1.77777778, 1.8984375 ])
>>> # Pythagorean intervals:
>>> librosa.pythagorean_intervals(bins_per_octave=12)
array([1. , 1.06787109, 1.125 , 1.20135498, 1.265625 ,
1.35152435, 1.42382812, 1.5 , 1.60180664, 1.6875 ,
1.80203247, 1.8984375 ])
>>> # 12-TET intervals:
>>> 2**(np.arange(12)/12)
array([1. , 1.05946309, 1.12246205, 1.18920712, 1.25992105,
1.33483985, 1.41421356, 1.49830708, 1.58740105, 1.68179283,
1.78179744, 1.88774863])
Create a 7-bin, 5-limit interval set
>>> librosa.plimit_intervals(primes=[3, 5], bins_per_octave=7)
array([1. , 1.125 , 1.25 , 1.33333333, 1.5 ,
1.66666667, 1.875 ])
The same example, but now in factored form
>>> librosa.plimit_intervals(primes=[3, 5], bins_per_octave=7,
... return_factors=True)
[
{},
{2: -3, 3: 2},
{2: -2, 5: 1},
{2: 2, 3: -1},
{2: -1, 3: 1},
{3: -1, 5: 1},
{2: -3, 3: 1, 5: 1}
]
"""
primes = np.atleast_1d(primes)
logs = np.log2(primes, dtype=np.float64)
# The seed set are primes and their reciprocals
# These are the values that we can use to expand our
# interval set. These are expressed in terms of the
# prime factorization exponents
seeds = []
for i in range(len(primes)):
# Add the prime
seed = [0] * len(primes)
seed[i] = 1
seeds.append(tuple(seed))
# Add the inverse
seed[i] = -1
seeds.append(tuple(seed))
# The frontier is the set of candidate intervals for inclusion
frontier = seeds.copy()
# The distances table will let us keep track of the harmonic
# distances between all selected intervals
distances = dict()
# Initialize the interval set with the root (1)
intervals = list()
root = tuple([0] * len(primes))
intervals.append(root)
while len(intervals) < bins_per_octave:
# Find the element on the frontier that minimizes the total
# harmonic distance to the existing set
score = np.inf
best_f = 0
for f, point in enumerate(frontier):
# Compute harmonic distance (HD) to each selected interval
HD = 0.0
for s in intervals:
if (s, point) not in distances:
distances[s, point] = __harmonic_distance(logs, point, s)
distances[point, s] = distances[s, point]
HD += distances[s, point]
if HD < score or (
np.isclose(HD, score)
and _crystal_tie_break(point, frontier[best_f], logs)
):
score = HD
best_f = f
new_point = frontier.pop(best_f)
intervals.append(new_point)
for _ in seeds:
new_seed = tuple(np.array(new_point) + np.array(_))
if new_seed not in intervals and new_seed not in frontier:
frontier.append(new_seed)
pows = np.array(list(intervals), dtype=float)
log_ratios: np.ndarray
pow2: np.ndarray
log_ratios, pow2 = np.modf(pows.dot(logs))
# If the fractional part is negative, add
# one more power of two to get it into the range [0, 1).
too_small = log_ratios < 0
log_ratios[too_small] += 1
pow2[too_small] -= 1
# Convert powers of 2 to integer
pow2 = pow2.astype(int)
idx: Iterable[int]
if sort:
# Order the intervals
idx = np.argsort(log_ratios)
log_ratios = log_ratios[idx]
else:
# If not sorting, we'll take powers in order
idx = range(bins_per_octave)
if return_factors:
# Collect the factorized intervals into a list
factors = []
for i in idx:
v = dict()
if pow2[i] != 0:
v[2] = -pow2[i]
v.update({p: int(power) for p, power in zip(primes, pows[i]) if power != 0})
factors.append(v)
return factors
# Otherwise, just return intervals as floats
return np.power(2, log_ratios)