It turns out that the trapezoidal method is inaccurate when used on the linear data. The piece-wise linear approximation is only accurate on the log-log graph. Since I havn't found a nice detailed source on the subject, I did the calulations.
Function \$\mathscr{L}(f)\$ passing through \$(f_1,\mathscr{L}_1)\$ and \$(f_2,\mathscr{L}_2)\$, as a line on a log-log figure:
\begin{equation}
log_{10}(\mathscr{L}(f))=\frac{1}{10}\frac{\mathscr{L}_{2|dBc}-\mathscr{L}_{1|dBc}}{log_{10}(f_2)-log_{10}(f_1)}log_{10}(f)+\frac{1}{10}\frac{log_{10}(f_2)\mathscr{L}_{1|dBc}-log_{10}(f_1)\mathscr{L}_{2|dBc}}{log_{10}(f_2)-log_{10}(f_1)}
\end{equation}
This can be rewritten as:
\begin{equation}
\begin{aligned}
\mathscr{L}(f) &= 10^{a \cdot log_{10}(f)+b}=10^{log_{10}(f^a)}10^b=f^a10^b \\
a &= \frac{\mathscr{L}_{2|dBc}-\mathscr{L}_{1|dBc}}{10(log_{10}(f_2)-log_{10}(f_1))} \qquad b = \frac{log_{10}(f_2)\mathscr{L}_{1|dBc}-log_{10}(f_1)\mathscr{L}_{2|dBc}}{10(log_{10}(f_2)-log_{10}(f_1))}
\end{aligned}
\end{equation}
Integration of \$\mathscr{L}(f)\$ from \$f_1\$ to \$f_2\$:
\begin{equation}
\int_{f_1}^{f_2} 2\mathscr{L}(f)df=2\cdot10^b\int_{f_1}^{f_2} f^a df=2\frac{f_2^{a+1}-f_1^{a+1}}{a+1}10^b
\end{equation}
With multiple segments:
\begin{equation}
\begin{aligned}
&\int 2\mathscr{L}(f)df = 2 \sum \frac{f_{i+1}^{a_i+1}-f_i^{a_i+1}}{a_i+1}10^{b_i} \\
a_i &= \frac{\mathscr{L}_{i+1|dBc}-\mathscr{L}_{i|dBc}}{10(log_{10}(f_{i+1})-log_{10}(f_i))} \qquad b_i = \frac{log_{10}(f_{i+1})\mathscr{L}_{i|dBc}-log_{10}(f_i)\mathscr{L}_{i+1|dBc}}{10(log_{10}(f_{i+1})-log_{10}(f_i))}
\end{aligned}
\end{equation}
The python code is changed to:
import numpy as np
f = np.array([1e0,1e1, 1e2, 1e3, 1e4, 1e5])
Lf = np.array([-50,-70,-113,-128,-135,-140])
f0 = 10e6
# f_{i}, f_{i+1}, L_{i}, L_{i+1}
fi = f[:-1]
fip1 = f[1:]
Li = Lf[:-1]
Lip1 = Lf[1:]
ai = (Lip1-Li) / (np.log10(fip1) - np.log10(fi)) /10.0
bi = (Li*np.log10(fip1)-Lip1*np.log10(fi)) / (np.log10(fip1) - np.log10(fi)) /10.0
Sphi = 2*np.sum( 10.0**(bi) * (fip1**(ai+1)-fi**(ai+1))/(ai+1) )
jitter = 1.0/2/np.pi/f0 * np.sqrt(Sphi)
print jitter * 1e12,
print 'ps'
raw_input('done')
The output is 68.66 ps, it matches the online tools.
I'm still interested if somebody have a more complete source on this subject or a simpler method.
