What function describes the frequency for each unique ratio for all possible expansions n over d where n<d?

Question

I am hoping to solve the following problem for a scientific investigation, which relies on the probabilites of all possible expansions.

What function $f(r)$ describes the frequency for each ratio for all possible combinations of positive integers*

$$\frac{n}{d}$$

where $n$ and $d$ are finite integers $$n\wedge d \in\mathbb Z$$ ($n$ is the numerator and $d$ is the denominator) such that:

$$n<d$$

Example

When $d_{max} = 16$ the ratio $\frac{n}{d}=0.5$ can be expressed as $\frac{8}{16},\frac{4}{8},\frac{2}{4},\frac{1}{2}$ as well as $\frac{3}{6},\frac{5}{10},\frac{6}{12},\frac{7}{14}$

The number of possible expansions sum to $8$. Thus, the frequency of possible expansions when $r=0.5$ is $$f(0.5)=8$$

---

Given $n<d$ and $d_{max}=16$

What is $f(r)$?

_{*For this particular problem a solution with $d_{max}=16$ is sufficient}

---

Background Reading for the Mathematically Inclined

I found this paper by Everett (1946), where f-expansions for the special case $x/p$ are mentioned, but it is too complex for someone with my limited background to digest. It does describe "correspondence from real numbers to sequences of integers mod p" (p.1).

---

Programmatical Solution in MATLAB

_EDIT

Inspired by the reply to Number of common divisors between two given numbers I modified the code to work with MATLAB. The resulting output is presented in the figure below the code.

% MATLAB CODE
% Based on C++ math.stackexchange.com question 8611 (VelvetThunder)

dmax = 16;

distribution = zeros(dmax,1);

for j = 1:(dmax-1)

    [N,C,D] = gcd(j,dmax);
    answer = 0;
    sqt = sqrt(N);
    for i = 1:sqt
        if mod(N,i) == 0
              answer = answer + 2;
        end
        if i == N/i
            answer = answer - 1;
        end
        distribution(j) = answer;
    end

    fprintf('%d\\%d = %d\n', j, dmax, answer);

    bar(1:length(distribution),distribution);

end

Commentary

The programmatical (inductive) solution appears to yield, given the mentioned constraints, a modal peak when $n = \frac{d_{max}}{2}$ and a stepwise up-and-down progression. However, for $d_{max}=128$ the pattern becomes more intricate:

For inductive purposes, let us not forget to change the radix from 2, by setting $d_{max} = 96$. The $odd=1$ pattern is now gone, but the maximum number of expansions still appears to be located at the half-width point:

The question how to express this - formally - remains.