Float precision when dealing with logarithms in musical pitch calculations-CodePudding

I'm writing a simple program to determine the difference between two musical pitches in cents; one cent is equal to 1/100th of a semitone. Dealing in cents is preferable for comparing musical pitches because the frequency scale is logarithmic, not linear. In theory, this is an easy calculation: the formula for determining the number of cents between two frequencies is:

1200 * log2(pitch_a / pitch_b)

I've written a small piece of code to automate this process:

import numpy as np
import math

def cent_difference(pitch_a, pitch_b)
     cents = 1200 * np.abs(math.log2(pitch_a / pitch_b))
     return cents

This works perfectly when I give the program octaves:

In [28]: cent_difference(880, 440)
Out[28]: 1200.0

...but misses the mark by about two cents on a perfect fifth:

In [29]: cent_difference(660, 440)
Out[29]: 701.9550008653875

...and keeps getting worse as I go, missing by about 14 cents on a major third:

In [30]: cent_difference(550, 440)
Out[30]: 386.31371386483477

Is this all float precision nonsense? Why does the perfect 5th example overestimate the cents, but the major third example underestimate the cents? What's going on here?

Much obliged for any help!

CodePudding user response：

The issue you're having isn't about the accuracy of Python's float type, but about the discrepancy between equal temperament and just intonation in music.

>>> cent_difference(660, 440)
701.9550008653874

This is assuming that a P5 interval represents a frequency ratio of 3/2. But in 12-ET, it doesn't: It has a ratio of 2^7/12 ≈ 1.4983070768766815. With the proper ET value for the higher note, you do get the expected 700.

>>> cent_difference(659.2551138257398, 440)
700.0

CodePudding user response：

What's going on here?

You're inputting frequency intervals in just intonation and expecting results in equal temperament..

If you feed the equal-tempered major third frequency ratio of 2^(4/12) into your formula, you indeed get a result of 400 cents (within floating point accuracy, as explained by the other answers and comments).

CodePudding user response：

The problem here is that floating point numbers use a set number of bits to represent any of the real numbers. Since there's infinitely many of those and only 2**32 values for a 32-bit float (at best), you can see how there will effectively be infinitely many reals that will have to be approximated. And if you keep computing with those approximations, errors occur.

You don't have to use large or long numbers to run into one either. My favourite:

>>> .1   .1   .1
0.30000000000000004

You can use more accurate types, which use better representations at the cost of some speed (and sometimes use operations that are slower but less likely to introduce errors).

For example Decimal, but make sure you use integers to define them:

>>> .1   .1   .1
0.30000000000000004
>>> from decimal import Decimal
>>> Decimal(.1)   Decimal(.1)   Decimal(.1)
Decimal('0.3000000000000000166533453694')
>>> Decimal (1)/Decimal(10)   Decimal(1)/Decimal(10)   Decimal(1)/Decimal(10)
Decimal('0.3')

The best solution, if one exists for your problem, is to avoid floating point math altogether.

By the way, here's your problem using Decimal:

from decimal import Decimal, Context


def cent_difference(pitch_a, pitch_b, ctx):
    ratio = ctx.divide(pitch_a, pitch_b)
    cents = Decimal(1200) * ctx.copy_abs(ratio.ln(ctx) / Decimal(2).ln(ctx))
    return cents


ctx = Context(prec=20)
print(cent_difference(Decimal(880), Decimal(440), ctx))
print(cent_difference(Decimal(660), Decimal(440), ctx))

Result:

1200
701.95500086538741774000

So, not that different. I'm not sure what result you expected on the second one there. If you pop over to Wolfram Alpha and task it with 1200 * log2(660 / 440), there seems to be no clean way to write this without the log still in there - precision will be lost for any representation in digits of an irrational number.