Background
If $v_1, v_2, \dots, v_n$ is an orthonormal basis in the inner product space $V$, then any vector $u\in V$ can be expressed as a linear combination
$$
u = \alpha_1 v_1 + \alpha_2 v_2 + \dots + \alpha_n v_n.\tag1
$$
Moreover, the coefficients can be computed using $\alpha_k=\langle v_k, u\rangle$, as can be seen by applying $\langle v_k, .\rangle$ to both sides of $(1)$.
Fidelity in terms of the characteristic function
The set $L(\mathcal{H})$ of linear operators on a $d$-dimensional Hilbert space $\mathcal{H}$ forms an inner product space with the inner product defined as
$$
\langle A, B\rangle = \mathrm{tr}(A^\dagger B).
$$
It is easy to check that the normalized Pauli operators $B_k = W_k/\sqrt{d}$ form an orthonormal basis in $L(\mathcal{H})$ with respect to this inner product. Therefore, any operator $\rho \in L(\mathcal{H})$ can be written as
$$
\rho = \alpha_1 B_1 + \alpha_2 B_2 + \dots + \alpha_{d^2} B_{d^2}\tag{1'}
$$
and the coefficients can be computed as
$$
\alpha_k = \langle B_k, \rho\rangle = \mathrm{tr}(B_k^\dagger \rho) = \mathrm{tr}(\rho B_k) = \chi_\rho(k).\tag2
$$
Finally, using $(1')$ and $(2)$, we find
$$
\begin{align}
\rho = & \sum_{i=1}^{d^2}\chi_\rho(i)B_i \\
\rho\sigma = & \sum_{i=1}^{d^2}\chi_\rho(i)B_i\sigma \\
\mathrm{tr}(\rho\sigma) = & \mathrm{tr}\left(\sum_{i=1}^{d^2}\chi_\rho(i)B_i\sigma\right) \\
\mathrm{tr}(\rho\sigma) = & \sum_{i=1}^{d^2}\chi_\rho(i)\,\mathrm{tr}\left(B_i\sigma\right) \\
\mathrm{tr}(\rho\sigma) = & \sum_{i=1}^{d^2}\chi_\rho(i)\chi_\sigma(i)
\end{align}
$$
which is the desired equality.
I am not aware of any connection between the characteristic function defined in the paper and the characteristic function of a random variable.
