Here is a proof of the conjecture, a proof that also shows how to compute the integrals explicitly.
The proof is somewhat similar to David Speyer's approach, but instead of using multivariable residues, I will just shift a one-variable contour. Without loss of generality, $m>0$. Eliminating the trigonometric functions yields $I(m,n)=(-1/4\pi^2) J(m,n)$, where
$$ J(m,n) := \int_{|w|=1} \int_{|z|=1} \log(4+z+z^{-1}+w+w^{-1}) \, z^{m-1} w^{n-1} \, dz \, dw.$$
For fixed $w \ne -1$ on the unit circle, the values of $z$ such that $4+z+z^{-1}+w+w^{-1}$ are a negative real number and its inverse; let $g(w)$ be the value in $(-1,0)$. The function $\log(4+z+z^{-1}+w+w^{-1})$ of $z$ extends to the complement of the closed interval $[g(w),0]$ in a disk $|z|<1+\epsilon$. Shrink the contour $|z|=1$ like a rubber band so that it hugs the closed interval. The upper and lower parts nearly cancel, leaving
$$ J(m,n) = \int_{|w|=1} \int_{z=0}^{g(w)} -2\pi i \, z^{m-1} w^{n-1} \,dz \,dw,$$
with the $-2\pi i$ coming from the discrepancy in the values of $\log$ on either side of the closed interval. Thus
$$ J(m,n) = -\frac{2 \pi i}{m} \int_{|w|=1} g(w)^m w^{n-1} \, dw.$$
Next use the rational parametrization $(g(w),w) = \left( (1-u)/(u+u^2), (u^2-u)/(1+u) \right)$ that David Speyer wrote out. There is a path $\gamma$ from $-i$ to $i$ in the right half of the $u$-plane that maps to $|w|=1$ (counterclockwise from $-1$ to itself) and gives the correct $g(w)$. Integrating the resulting rational function of $u$ gives
$$ J(m,n) = -\frac{2 \pi i}{m} \left.(f(u) + r \log u + s \log(1+u))\right|_\gamma = -\frac{2 \pi i}{m} (f(i)-f(-i) + r \pi i + s \pi i/2),$$
for some $f \in \mathbf{Q}(u)$ and $r,s \in \mathbf{Q}$. This shows that $I(m,n)=a + b/\pi$ for some $a,b \in \mathbf{Q}$.
Examples: Mathematica got the answers wrong, presumably because it doesn't understand homotopy very well! Here are some actual values, computed using the approach above:
$$I(1,1) = \frac{1}{2} - \frac{2}{\pi}$$
$$I(2,1) = -1 + \frac{10}{3\pi}$$
$$I(20,10) = -\frac{14826977006}{5} + \frac{56979906453888224582}{6116306625 \pi}.$$
Best Answer
It suffices that $f$ be (the restriction to $[0,2\pi]$ of) a completely monotone real-valued function defined on $[0,\infty)$. Indeed, then for some finite measure $\mu$ on $[0,\infty)$ and all real $x\ge0$ we have $$f(x)=\int_0^\infty\mu(da) e^{-a x},$$ whence for natural $n$ $$\hat f(n)=\int_0^\infty\mu(da) \int_0^{2\pi}dx\,\cos(nx)e^{-a x} =\int_0^\infty\mu(da) \frac{a \left(1-e^{-2 \pi a}\right)}{a^2+n^2},$$ which is obviously decreasing in $n$ (to $0$, by dominated convergence or by the Riemann--Lebesgue lemma).
Note that, if $f(x)\equiv1$ or $f(x)\equiv x$, then $\hat f(n)=0$ for all natural $n$. So, if $f$ has the desired property, then the function $[0,2\pi]\ni x\mapsto a+bx+f(x)$ also has it for any real $a$ and $b$. Also, clearly, if $f$ has the desired property, then do does the function $$[0,2\pi]\ni x\mapsto f^-(x):=f(2\pi-x)$$ -- because $\widehat{f^-}(n)=\hat f(n)$ for all natural $n$. It follows that, if $f$ and $g$ have the desired property, then the function $$[0,2\pi]\ni x\mapsto a+bx+f(x)+g(2\pi-x)$$ also has it for any real $a$ and $b$.
Added:
As noted in a comment by Fedor Petrov, if $f(x)=h(\pi-x)$ for some odd function $h$ and all $x\in[0,2\pi]$, then $\hat f(n)=0$ for all natural $n$.
It follows from this answer by fedja that, if $$f(x)=\int_1^\infty[\mu(dp) x^p+\nu(dp)(2\pi-x)^p]<\infty$$ for some measures $\mu$ and $\nu$ on $[1,\infty)$ and all $x\in[0,2\pi]$, then $f$ has the desired property.