Author Archives: Robert Schwarz

Consonant Triads

We use a physiological interpretation of consonance and dissonant musical sound for human listeners and use a dissonance score to rate different possible triads (combinations of three notes in a scale). This is first done on the conventional scale of 12 equidistant tones in an octave, as well as the Bohlen-Pierce scale.

In [1]:
using Plots
gr()
Out[1]:
Plots.GRBackend()

Plomp-Levelt dissonance curve

As a basis of our investigation, we use the Plomp-Levelt dissonance curve as explained in this article by William Sethares.

The curve models the fact that similar frequencies appear dissonant because of “friction” between the adjacent sensors in our ears.

In [2]:
"Dissonance of sound with two sinusoidal sounds with given frequencies and amplitudes."
function diss(freq1, ampl1, freq2, ampl2)
    if freq2 < freq1 # swap args
        return diss(freq2, ampl2, freq1, ampl1)
    end
    s = 0.24 * (freq2 - freq1) / (0.0207 * freq1 + 18.96)
    5.0 * ampl1 * ampl2 * (exp(-3.51 * s) - exp(-5.75 * s))
end
Out[2]:
diss

Let’s draw the curve for a pair of sine waves, the first fixed at given base frequency and the other going through a range of relatively higher frequencies.

In [3]:
basefreq = 440.0
ratrange = range(1.0, stop=3.0, length=200)
otherfreq = basefreq * ratrange
disss = diss.(basefreq, 1.0, otherfreq, 1.0);
In [4]:
plot(otherfreq, disss, xlabel="Frequency", ylabel="Dissonance", legend=false)
Out[4]:

600 800 1000 1200 0.0 0.2 0.4 0.6 0.8 Frequency Dissonance

Dissonance with harmonics

Rather than analyzing dissonance between sinusoidal sounds, we also want to consider sounds with rich harmonics. These have sounds at multiple frequencies with declining amplitudes. When computing the dissonance between two pitches with harmonics, we have to add the dissonance values for all combinations of frequencies.

In [5]:
# Let's first give formulas for the amplitudes of common waveshapes.

enumrange(number::Int) = collect(range(1.0, stop=number, length=number))

abstract type Wave end

struct Sine <: Wave end
harm(::Sine, number::Int) = vcat([1.0], zeros(number - 1))

struct Sawtooth <: Wave end
harm(::Sawtooth, number::Int) = 1.0 ./ enumrange(number)

struct Triangle <: Wave end
harm(::Triangle, number::Int) = [(((i+2)%4)-2)*(i%2) for i=1:number] ./ enumrange(number).^2

struct Square <: Wave end
harm(::Square, number::Int) = [i%2 for i in 1:number] ./ enumrange(number)
Out[5]:
harm (generic function with 4 methods)
In [6]:
@show harm(Sine(), 9)
@show harm(Sawtooth(), 9)
@show harm(Square(), 9)
@show harm(Triangle(), 9);

harm(Sine(), 9) = [1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
harm(Sawtooth(), 9) = [1.0, 0.5, 0.333333, 0.25, 0.2, 0.166667, 0.142857, 0.125, 0.111111]
harm(Square(), 9) = [1.0, 0.0, 0.333333, 0.0, 0.2, 0.0, 0.142857, 0.0, 0.111111]
harm(Triangle(), 9) = [1.0, 0.0, -0.111111, 0.0, 0.04, 0.0, -0.0204082, 0.0, 0.0123457]

Let’s also quickly check how our approximated waveforms look:

In [7]:
xx = range(0, stop=4π, length=100)

function yy(wave, n)
    coefs = harm(wave, n)
    sum(coefs[f] * sin.(f*xx) for f=1:n)
end

plot(xx,[yy(w, 20) for w in [Sawtooth(), Square(), Triangle()]])
Out[7]:

0 2 4 6 8 10 12 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 y1 y2 y3

In [8]:
# We add a new type to couple frequency and amplitute of a harmonic overtone

struct Tone
    freq::Float64
    ampl::Float64
end

diss(t1::Tone, t2::Tone) = diss(t1.freq, t1.ampl, t2.freq, t2.ampl)
Out[8]:
diss (generic function with 2 methods)
In [9]:
# More functions to generate all overtones

function harmonics(wave::Wave, basefreq::Float64, number::Int)
    # frequencies with integer multiples
    freqs = basefreq .* enumrange(number)
    ampls = abs.(harm(wave, number))
    # keep tones with nonzero amplitude
    [Tone(tup...) for tup in zip(freqs, ampls) if tup[2] > 0.0]
end
Out[9]:
harmonics (generic function with 1 method)
In [10]:
@show harmonics(Sine(), 440.0, 4)
@show harmonics(Triangle(), 440.0, 4);

harmonics(Sine(), 440.0, 4) = Tone[Tone(440.0, 1.0)]
harmonics(Triangle(), 440.0, 4) = Tone[Tone(440.0, 1.0), Tone(1320.0, 0.111111)]
In [11]:
function diss(tones::AbstractArray{Tone})
    res = 0.0
    # iterate over all pairs
    n = length(tones)
    for i = 1:n
        for j=i+1:n
            res += diss(tones[i], tones[j])
        end
    end
    res
end
Out[11]:
diss (generic function with 3 methods)
In [12]:
function together(tones1::AbstractArray{Tone}, tones2::AbstractArray{Tone})
    # TODO: should do proper merge and match those tones with same frequency?
    vcat(tones1, tones2)
end

function diss(tones1::AbstractArray{Tone}, tones2::AbstractArray{Tone})
    diss(together(tones1, tones2))
end
Out[12]:
diss (generic function with 4 methods)

As a first example, we look at the square waveform approximated with 10 harmonices and draw the dissonance curve for two pitches played simultaneously.

Note the local maxima near frequency rations 1 and 2 and the minima at integer (or “simple”) ratios.

In [13]:
numharm = 10
basetone = harmonics(Square(), basefreq, numharm)
othertones = [harmonics(Square(), of, numharm) for of in otherfreq]
disss = [diss(basetone, ot) for ot in othertones]

plot(ratrange, disss, legend=false, xlabel="Freq ratio", ylabel="Dissonance", title="Square $basefreq")
Out[13]:

1.0 1.5 2.0 2.5 3.0 0.00 0.25 0.50 0.75 1.00 Square 440.0 Freq ratio Dissonance

Next, we repeat the experiment, now based on an approximation of the sawtooth waveform.
Because of the richer harmonices, there is more interaction and more dissonance (in absolute terms).

As a consequence, the dissonance curve shows more local minima, again corresponding to “simple” frequency ratios, that is, ratios of small integers.

In [14]:
numharm = 10
basetone = harmonics(Sawtooth(), basefreq, numharm)
othertones = [harmonics(Sawtooth(), of, numharm) for of in otherfreq]
disss = [diss(basetone, ot) for ot in othertones]

plot(ratrange, disss, legend=false, xlabel="Freq ratio", ylabel="Dissonance", title="Sawtooth $basefreq")
Out[14]:

1.0 1.5 2.0 2.5 3.0 0.00 0.25 0.50 0.75 1.00 1.25 Sawtooth 440.0 Freq ratio Dissonance

Triads

We will now compute the dissonaces for triads of tones in an analoguous fashion.
But this time, we will only evaluate a discrete set of pitches, based on a 12-tone scale, dividing the octave in equal steps.

In [15]:
# dissonance for triads
function diss(tones1::AbstractArray{Tone}, tones2::AbstractArray{Tone}, tones3::AbstractArray{Tone})
    diss(together(together(tones1, tones2), tones3))
end
Out[15]:
diss (generic function with 5 methods)
In [16]:
# range of one octave
halfstep = 2.0^(1.0/12.0)
@show halfstep
ratrange = [halfstep^i for i=0:12]

basefreq = 440.0
otherfreq = basefreq * ratrange;

halfstep = 1.0594630943592953

We start using the square waveform. In the heatmap below, we can see that dissonance is high whenever a single halfstep is present in the triad. This is true for the horizontal and vertical lines that are off by one from the origin, as well as the two diagonals off by one from the center. To a lesser extent, this pattern repeats for 11 halfsteps (one octave up, one halfstep down).

In [17]:
# start with square wave again
numharm = 10
basetone = harmonics(Square(), basefreq, numharm)
othertones = [harmonics(Square(), of, numharm) for of in otherfreq]

# compute dissonances of all triads
disss = [diss(basetone, ot1, ot2) for ot1 in othertones, ot2 in othertones]

heatmap(0:12, 0:12, disss, aspect_ratio=1.0, title="Square wave on 12-tone octave")
Out[17]:

0 2 4 6 8 10 12 0 2 4 6 8 10 12 Square wave on 12-tone octave 0.5 1.0 1.5 2.0 2.5

If we use the sawtooth wave instead, we see a richer structure in the dissonance heatmap.
Here, we should also be able to recognize familiar concepts such as the major or minor triads.

In [18]:
# now use richer sawtooth wave
numharm = 10
basetone = harmonics(Sawtooth(), basefreq, numharm)
othertones = [harmonics(Sawtooth(), of, numharm) for of in otherfreq]

# compute dissonances of all triads
disss = [diss(basetone, ot1, ot2) for ot1 in othertones, ot2 in othertones]

heatmap(0:12, 0:12, disss, aspect_ratio=1.0, title="Sawtooth wave on 12-tone octave")
Out[18]:

0 2 4 6 8 10 12 0 2 4 6 8 10 12 Sawtooth wave on 12-tone octave 0.5 1.0 1.5 2.0 2.5 3.0 3.5

Bohlen-Pierce Scale

One way to view the conventional western scale is as a division of an octave (that is, a 2:1 ratio of frequency) in 12 (more or less) equal steps. Alternatively, one can start at the ratios 4:5:6 (appearing naturally as harmonics) and collect pitches that are reached from these recursively.

In the Bohlen-Pierce scale, we range over a tritave (a 3:1 ratio) in 13 steps. This can also be derived using ratios 3:5:7. The scale is well-suited for instruments that only feature odd harmonics, such as the clarinet or the pan flute.

In [19]:
# range of one tritave
bpstep = 3.0^(1.0/13.0)
@show bpstep
ratrange = [bpstep^i for i=0:13]

basefreq = 440.0
otherfreq = basefreq * ratrange;

bpstep = 1.0881822434633168
In [20]:
# start with square wave again
numharm = 10
basetone = harmonics(Square(), basefreq, numharm)
othertones = [harmonics(Square(), of, numharm) for of in otherfreq]

# compute dissonances of all triads
disss = [diss(basetone, ot1, ot2) for ot1 in othertones, ot2 in othertones]

heatmap(0:13, 0:13, disss, aspect_ratio=1.0, title="Square wave on BP tritave")
Out[20]:

0 2 4 6 8 10 12 0 2 4 6 8 10 12 Square wave on BP tritave 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00

The heatmap above looks fairly boring, with not much dissonance apart from the obvious lines. Maybe this is related to the wider spacing of the steps in this scale.

We will also look at the sawtooth wave, even though we should expect high dissonance through-out.

In [21]:
# now use richer sawtooth wave
numharm = 10
basetone = harmonics(Sawtooth(), basefreq, numharm)
othertones = [harmonics(Sawtooth(), of, numharm) for of in otherfreq]

# compute dissonances of all triads
disss = [diss(basetone, ot1, ot2) for ot1 in othertones, ot2 in othertones]

heatmap(0:13, 0:13, disss, aspect_ratio=1.0, title="Sawtooth wave on BP tritave")
Out[21]:

0 2 4 6 8 10 12 0 2 4 6 8 10 12 Sawtooth wave on BP tritave 0.5 1.0 1.5 2.0 2.5

Larger Pitch Ranges

We will generate some more heatmaps covering two octaves or tritaves, each time based on sawtooth waves.

First the 12-tone octaves:

In [22]:
# range of two octaves
ratrange = [halfstep^i for i=0:24]
basefreq = 440.0
otherfreq = basefreq * ratrange

# sawtooth wave
numharm = 10
basetone = harmonics(Sawtooth(), basefreq, numharm)
othertones = [harmonics(Sawtooth(), of, numharm) for of in otherfreq]

# compute dissonances of all triads
disss = [diss(basetone, ot1, ot2) for ot1 in othertones, ot2 in othertones]

heatmap(0:24, 0:24, disss, aspect_ratio=1.0, title="Sawtooth wave on 12-tone octaves")
Out[22]:

0 5 10 15 20 0 5 10 15 20 Sawtooth wave on 12-tone octaves 0.5 1.0 1.5 2.0 2.5 3.0 3.5

Then the BP tritaves:

In [23]:
# range of two tritaves
ratrange = [bpstep^i for i=0:26]
basefreq = 440.0
otherfreq = basefreq * ratrange

# sawtooth wave
numharm = 10
basetone = harmonics(Sawtooth(), basefreq, numharm)
othertones = [harmonics(Sawtooth(), of, numharm) for of in otherfreq]

# compute dissonances of all triads
disss = [diss(basetone, ot1, ot2) for ot1 in othertones, ot2 in othertones]

heatmap(0:26, 0:26, disss, aspect_ratio=1.0, title="Sawtooth wave on BP tritaves")
Out[23]:

0 5 10 15 20 25 0 5 10 15 20 25 Sawtooth wave on BP tritaves 0.5 1.0 1.5 2.0 2.5

Marginal Net Income in Germany

We investigate the tax burden and insurance cost for different scenarios of gross income using a simplified model of the situation in Germany. The goal is to compute the marginal net income for additional gross income. From that, we can also derive the net hourly income for part-time employment.

This analysis is done using Julia in a Jupyter notebook.

In [1]:
# assume a typical monthly income
gross_monthly = 37873 / 12
Out[1]:
3156.0833333333335

Social costs

For simplicity, we ignore the possibility of private health insurance and assume that our typical employee participates in the public social system. This includes health insurance, unemployment insurance and payments to the pension system.

These payments are split evenly between the employer and employee. This means that the labor cost of the employer is higher than the gross salary paid.

Source: nettolohn.de

In [2]:
# only includes the employee's share
pension = 0.186
unemployment = 0.03
health = 0.146
longterm_care = 0.028
base_rate = 0.5 * (pension + unemployment + health + longterm_care)
health_extra = 0.011
employer_rate = base_rate
employee_rate = base_rate + health_extra
@show employer_rate, employee_rate;

(employer_rate, employee_rate) = (0.195, 0.20600000000000002)
In [3]:
# effective labor cost for employer
real_monthly = (1.0 + employer_rate) * gross_monthly
Out[3]:
3771.5195833333337

However, these relative rates are only applied up to specific bound for the gross income.
So, let’s define some functions to compute the correct costs.

In the case of self-employment, there is also a virtual minimum income that is used as a reference, which we will ignore for simplicity.

Source: Beitragsbemessungsgrenze

In [4]:
function social_cost(gross_monthly)
    pension       = 0.186 * min(gross_monthly, 6500.0)
    unemployment  = 0.030 * min(gross_monthly, 6500.0)
    health        = 0.146 * min(gross_monthly, 4425.0)
    longterm_care = 0.146 * min(gross_monthly, 4425.0)
    pension + unemployment + health + longterm_care
end

employer_cost(gross_monthly) = 0.5 * social_cost(gross_monthly)
employee_cost(gross_monthly) = 0.5 * social_cost(gross_monthly) + 0.011 * gross_monthly
total_cost(gross_monthly) = employer_cost(gross_monthly) + employee_cost(gross_monthly)
employer_total(gross_monthly) = employer_cost(gross_monthly) + gross_monthly;

Let’s visualize the total social cost relative to the real monthly labor cost.

In [5]:
using Plots
gr()
Out[5]:
Plots.GRBackend()
In [6]:
gross_range = 400.0:20.0:10000.0
real_income = employer_total.(gross_range)
relative_cost = total_cost.(gross_range) ./ real_income
plot(real_income, relative_cost, xlim=(0,11000), ylim=(0.0, 0.5), legend=false,
     xlabel="effective monthly income", ylabel="rel social cost")
Out[6]:

0 2500 5000 7500 10000 0.0 0.1 0.2 0.3 0.4 0.5 effective monthly income rel social cost

So, higher gross (or effective) income leads to a smaller relative social cost.
We can repeat that plot for the employee’s point of view:

In [7]:
plot(gross_range, employee_cost.(gross_range)./gross_range, xlim=(0,10000), ylim=(0.0, 0.3), legend=false,
     xlabel="gross monthly income", ylabel="employee's share of social cost")
Out[7]:

0 2000 4000 6000 8000 10000 0.00 0.05 0.10 0.15 0.20 0.25 gross monthly income employee’s share of social cost

The income tax only applies to the part of the gross income of which the social cost has been subtracted already.
Further, some portion of that remainder is also free from taxes.

In [8]:
taxable_income(gross_monthly) = gross_monthly - employee_cost(gross_monthly)
Out[8]:
taxable_income (generic function with 1 method)

Income tax

The tax rate is a piece-wise defined function.
We assume an unmarried employee with no kids who is not taxable by any church.

The source below actually contains flow charts with details conditions, thresholds and rounding of intermediate results.
I can’t be bothered to understand all of that, so I will try to extract the essential information.

Source: BMF Steuerrechner

In [9]:
function income_tax(yearly)
    if yearly <= 9000
        return 0
    elseif yearly < 13997
        y = (yearly - 9000) / 10000
        rw = y * 997.8 + 1400
        return ceil(rw * y)
    elseif yearly < 54950
        y = (yearly - 13996) / 10000
        rw = y * 220.13 + 2397
        return ceil(rw * y + 948.49)
    elseif yearly < 260533
        return ceil(yearly * 0.42 - 8621.75)
    else
        return ceil(yearly * 0.45 - 16437.7)
    end
end
Out[9]:
income_tax (generic function with 1 method)
In [10]:
yearly_range = 5000:500:100000
plot(yearly_range, income_tax.(yearly_range), legend=false, xlim=(0,100000), ylim=(0,35000),
    xlabel="yearly taxable income", ylabel="income tax")
Out[10]:

0 2×10 4 4×10 4 6×10 4 8×10 4 1×10 5 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 2.5×10 4 3.0×10 4 3.5×10 4 yearly taxable income income tax

In [11]:
plot(yearly_range, income_tax.(yearly_range) ./ yearly_range, legend=false, xlim=(0,100000), ylim=(0,0.4),
    xlabel="yearly taxable income", ylabel="income tax rate")
Out[11]:

0 2×10 4 4×10 4 6×10 4 8×10 4 1×10 5 0.0 0.1 0.2 0.3 0.4 yearly taxable income income tax rate

Net income

Now, let’s combine the social costs and taxes to compute the net income.

In [12]:
function net(gross_monthly)
    taxable = taxable_income(gross_monthly)
    taxes = income_tax(ceil(12 * taxable)) / 12
    taxable - taxes
end

net_income = net.(gross_range)

plot(gross_range, net_income, legend=false, xlim=(0,10000), ylim=(0,6000),
     xlabel="monthly gross income", ylabel="monthly net income")
Out[12]:

0 2000 4000 6000 8000 10000 0 1000 2000 3000 4000 5000 6000 monthly gross income monthly net income

This looks surprisingly linear. Let’s also plot the relative net income compared to the effective cost.

In [13]:
rel_net = net_income ./ real_income
plot(real_income, rel_net, legend=false, xlim=(0,12000), ylim=(0,0.7),
     xlabel="effective monthly income", ylabel="rel net income")
Out[13]:

0 2.0×10 3 4.0×10 3 6.0×10 3 8.0×10 3 1.0×10 4 1.2×10 4 0.0 0.1 0.2 0.3 0.4 0.5 0.6 effective monthly income rel net income

Interestingly, this curve is not monotonically decreasing. In fact, there seems to be a minimum at a gross income of 4425 (real income of about 5300), which is the upper reference value for the health insurances.

Marginal income

The analysis above can be slightly misleading. If we are currently earning a certain income and have an opportunity to raise the income, this might also decrease our relative net income. However, higher rates for social cost and income tax are applied equally to the previous and additional income.

If we want to show how much net income we can retain for each additional EUR earned, we have take some more care.
Here, we approximate the slope of the net income as a function of real income using a symmetric difference quotient.
There are some jumps in that curve since our net income is not smooth.

In [14]:
finite_diffs = (net_income[3:end] - net_income[1:end-2]) ./ (real_income[3:end] - real_income[1:end-2])
plot(real_income[3:end], [finite_diffs rel_net[3:end]], xlim=(0,12000), ylim=(0,0.65),
    xlabel="effective monthly income", label=["marginal income" "relative net income"],
    legend=:bottomright)
Out[14]:

0 2.0×10 3 4.0×10 3 6.0×10 3 8.0×10 3 1.0×10 4 1.2×10 4 0.0 0.1 0.2 0.3 0.4 0.5 0.6 effective monthly income marginal income relative net income

Part-time work and hourly income

We have seen that a lower gross income often corresponds to a higher relative net income.
If we have the option to reduce our working hours, we should be able to increase our income per hour.

Let’s assume a 40 hour work-week and 4 weeks per month.

In [15]:
monthly_hours = 40.0 * 4.0
Out[15]:
160.0
In [16]:
# median income
median_net = net(gross_monthly)
Out[16]:
1929.0545833333333
In [17]:
median_net / monthly_hours
Out[17]:
12.056591145833334
In [18]:
using ColorBrewer
In [19]:
# hourly income for different monthly gross incomes
hours = 10:40
steps = 6
gross = range(1000, stop=8000, length=steps)'
colors = ColorBrewer.palette("YlGnBu", steps + 3)[4:end]
parttime_net = net.((hours / 40)*gross)
hourly_net = parttime_net ./ (hours * 4)
plot(hours, hourly_net, labels=gross, ylim=(0,40),
     color_palette=colors,
     xlabel="hours per week", ylabel="hourly net income",
     title="Hourly net income for part-time work")
Out[19]:

10 20 30 40 0 10 20 30 40 Hourly net income for part-time work hours per week hourly net income 1000.0 2400.0 3800.0 5200.0 6600.0 8000.0

As we can see above, we can increase our hourly income by working fewer hours every week.
This effect is more clear if we look at the hourly net income of part-time work relative the hourly net income for full-time work:

In [20]:
rel_hourly = hourly_net ./ hourly_net[end, :]';
In [21]:
plot(hours, rel_hourly, labels=gross, color_palette=colors,
     xlabel="hours per week", ylabel="relative hourly net income",
     title="Relative hourly net income for part-time work")
Out[21]:

10 20 30 40 1.00 1.05 1.10 1.15 1.20 Relative hourly net income for part-time work hours per week relative hourly net income 1000.0 2400.0 3800.0 5200.0 6600.0 8000.0

For low values of gross income, reducing the weekly work hours will not affect the hourly income any more.
Medium values show the most promise, with the largest inrease in hourly income.
With a high enough income, reducing the work load yields less and less increase in hourly income.