Edit 5.2.2014:
Using this as a starting point:
$$\Lambda(n)=\lim\limits_{s \rightarrow 1} \zeta(s)\sum\limits_{d|n} \frac{\mu(d)}{d^{(s-1)}}.$$
Where $\Lambda(n)$ is the von Mangoldt function.
Edit 25.3.2018 All of the pictures I posted above are either slightly or almost completely wrong. I deleted the last two. The correct plot starts from the following symbolic relationships:
Let $\mu(n)$ be the Möbius function, then:
$$a(n) = \sum\limits_{d|n} d \cdot \mu(d)$$
$$T(n,k)=a(GCD(n,k))$$
$$T = \left( \begin{array}{ccccccc} +1&+1&+1&+1&+1&+1&+1&\cdots \\ +1&-1&+1&-1&+1&-1&+1 \\ +1&+1&-2&+1&+1&-2&+1 \\ +1&-1&+1&-1&+1&-1&+1 \\ +1&+1&+1&+1&-4&+1&+1 \\ +1&-1&-2&-1&+1&+2&+1 \\ +1&+1&+1&+1&+1&+1&-6 \\ \vdots&&&&&&&\ddots \end{array} \right)$$
$$\sum\limits_{k=1}^{\infty}\sum\limits_{n=1}^{\infty} \frac{T(n,k)}{n^c \cdot k^s} = \sum\limits_{n=1}^{\infty} \frac{\lim\limits_{z \rightarrow s} \zeta(z)\sum\limits_{d|n} \frac{\mu(d)}{d^{(z-1)}}}{n^c} \;\;\;\;\;\;\;\;(1)$$
$$\sum\limits_{n=1}^{\infty} \frac{\lim\limits_{z \rightarrow s} \zeta(z)\sum\limits_{d|n} \frac{\mu(d)}{d^{(z-1)}}}{n^c} = \frac{\zeta(s) \cdot \zeta(c)}{\zeta(c + s - 1)} \;\;\;\;\;\;\;\;\;(2)$$
which is part of the limit:
$$\frac{\zeta '(s)}{\zeta (s)}=\lim_{c\to 1} \, \left( \color{Red}{\zeta (c)}-\frac{\zeta (c) \zeta (s)}{\zeta (c+s-1)}\right) \;\;\;\;\;\;\;(3)$$
$$\int \frac{\zeta '(s)}{\zeta (s)} ds = \log(\zeta(s)) \;\;\;\;\;\;\;\;\;\;(4)$$
And the zeta zero counting function is roughly:
$$f(t)=\frac{\vartheta (t)+\Im\left(\log \left(\zeta \left(i t+\frac{1}{2}\right)\right)\right)}{\pi }\;\;\;\;\;\;\;\;\;\;(5)$$
where $\vartheta (t)$ is the Riemann Siegel Theta function.
Since we are looking for the Fourier series like plot we instead write right hand side of $(3)$ as:
$$\frac{\zeta '(s)}{\zeta (s)} \approx \sum\limits_{n=1}^{\infty} \left( \color{Red}{\frac{1}{n^c}}-\frac{\lim\limits_{z \rightarrow s} \zeta(z)\sum\limits_{d|n} \frac{\mu(d)}{d^{(z-1)}}}{n^c}\right)\;\;\;\;\;\;\;\;\;\;(6)$$
Integrating $(6)$ we achieve:
$$\log(\zeta(s))\approx \int\sum\limits_{n=1}^{\infty} \left( \color{Red}{\frac{1}{n^c}}-\frac{\lim\limits_{z \rightarrow s} \zeta(z)\sum\limits_{d|n} \frac{\mu(d)}{d^{(z-1)}}}{n^c}\right)ds\;\;\;\;\;\;\;\;\;\;(7)$$
The colored $\color{Red}{\zeta(c)}$ in $(3)$ and $\color{Red}{\frac{1}{n^c}}$ in $(6)$ is equivalent to the first column in matrix $T$.
This integral can be solved symbolically with the Euler-Maclaurin formula for the Riemann zeta function, and combined with Riemann Siegel Theta in $(5)$ it is something like:
$f(t)=\small \frac{\vartheta (t)+\Re\left(\sum _{\text{nn}=1}^{\text{nnn}} \left(\frac{1}{\text{nn}^c}-\frac{\sum _{d=1}^{\text{nn}} \text{If}\left[\text{nn} \bmod d=0,\sum _{r=1}^{q-1} \left(\sum _{i=0}^{r-2} -\frac{B_r d^{1-\left(\frac{1}{2}+i t\right)} k^{-r-\left(\frac{1}{2}+i t\right)} \left|S_r^{(i)}\right| \text{Expand}\left[\sum _{m=0}^i \left(\frac{1}{2}+i t\right)^m (\log (d)+\log (k))^m N\left[\frac{i!}{m!}\right]\right]}{r! (\log (d)+\log (k))^{i+1}}\right)+\mu (d) \sum _{n=1}^k \text{If}\left[n=1\lor (n=1\land d=1),\frac{1}{2}+i t,-\frac{d^{1-\left(\frac{1}{2}+i t\right)}}{n^{\frac{1}{2}+i t} (\log (d)+\log (n))}\right]+\text{Ei}\left(-\left(\left(i t+\frac{1}{2}\right)-1\right) \log (d)-\left(\left(i t+\frac{1}{2}\right)-1\right) \log (k)\right),0\right]}{\text{nn}^c}\right)\right)}{\pi } \;(8)$
where $s$ has been replaced with $1/2+it$. $f(t)$ in $(8)$ is not computationally feasible to evaluate so instead we resort to looking at the numerical integral.
If we instead look at the plot of the numerical integral in the right hand side in $(9)$ below, then for $c=1$;
$$1+\frac{\vartheta (t)+\Im\left(\log \left(\zeta \left(i t+\frac{1}{2}\right)\right)\right)}{\pi } \\
\approx \frac{\vartheta (t)+\Re\left(\int\sum\limits_{n=1}^{\infty} \left( \frac{1}{n^c}-\frac{\lim\limits_{z \rightarrow \frac{1}{2} + it} \zeta(z)\sum\limits_{d|n} \frac{\mu(d)}{d^{(z-1)}}}{n^c}\right)ds\right)}{\pi}\;\;\;\;\;(9)$$
something much more Fourier series like appears:
where there is a jump of size one at every Riemann zeta zero.
The Mathematica program for the plot of the numerical integral in $(9)$ is:
(*start*)
Clear[n, s, c, t, nn, delta];
c = 1;
delta = 40;
nn = 60;
a = Sum[1/N[n]^c -
Zeta[1/2 + I*t]*
Total[MoebiusMu[Divisors[n]]/Divisors[n]^(1/2 + I*t - 1)]/
n^c, {n, 1, 600}];
ListLinePlot[(Table[
delta*RiemannSiegelTheta[t], {t, 0, nn, 1/delta}] +
Accumulate[Table[Re[a], {t, 0, nn, 1/delta}]])/delta/Pi,
PlotStyle -> Thickness[0.003], DataRange -> {0, nn}]
(*end*)
Edit 15.8.2018:
The correct way to integrate the Euler-Maclaurin formula is:
$$\boxed{N(t)=\frac{1}{\pi}\left(\vartheta (t)-\Re\left(\sum _{n=1}^{\text{nn}} \frac{1}{n^c} \left(\underset{d \mid n}{\sum\limits_{d=1}^{n}} \left(f(\frac{1}{2}+it, d)-f(\frac{1}{2}+i0, d) \right)\right)\right)\right)}$$
where $\vartheta (t)$ is the Riemann-Siegel theta function,
and where:
$$f(s,d)=-i \mu (d) \left(\sum _{n=0}^{q-1} \left(\sum _{i=0}^{2 n+1} -\frac{d B_{2 (n+1)} k^{-2 n-1} \left|S_{2 n+1}^{(i)}\right| \log ^{-i-1}(d k) \Gamma (i+1,s \log (k d))}{(2 (n+1))!}\right)+\sum _{n=2}^k -\frac{d^{1-s} n^{-s}}{\log (d)+\log (n)}+\text{If}\left[d=1 \text{ then } 0\text{ else }-\frac{1^{-s} d^{1-s}}{\log (d)+\log (1)}\right]+\frac{d^{1-s} k^{-s}}{2 (\log (d)+\log (k))}+\text{Ei}(-(s-1) \log (d)-(s-1) \log (k))\right)$$
and where $\mu (d)$ is the Möbius function.
As a Mathematica program the symbolic integral is:
(* Integral *)
(*start*)
Clear[k, n, s, i, d, r, q, f, a, b, integral, c, nn, n, k];
nn = 30;
delta = 1/5;
k = 30;
q = 10;
c = 1;
rr = 60;
f[s_, d_] = -I*
MoebiusMu[
d]*(If[d == 1, 0, -((d^(1 - s) 1^-s)/(Log[d] + Log[1]))] +
Sum[-((d^(1 - s) n^-s)/(Log[d] + Log[n])), {n, 2, k}] +
ExpIntegralEi[-(-1 + s) Log[d] - (-1 + s) Log[k]] + (
d^(1 - s) k^-s)/(2 (Log[d] + Log[k])) +
Sum[Sum[-(BernoulliB[2*(n + 1)]/((2*(n + 1))!))*(Abs[
StirlingS1[2*n + 1, i]]) d k^(-2*n - 1)
Gamma[1 + i, s Log[k*d]] Log[k*d]^(-1 - i), {i, 0,
2*n + 1}], {n, 0, q - 1}]);
integral[a_, b_, d_] = f[b, d] - f[a, d];
Print["Counting to ", rr];
Monitor[ListLinePlot[
Table[(RiemannSiegelTheta[t] -
Re[Sum[Sum[
If[Mod[n, d] == 0, integral[1/2 + I*0, 1/2 + I*t, d]/n^c,
0], {d, 1, n}], {n, 1, nn}]])/Pi, {t, 0, rr, delta}],
PlotStyle -> Thickness[0.003], DataRange -> {0, rr},
ImageSize -> Large], Floor[t]]
(* end *)