It is not what you want, but may be worth mentioning. There is a huge branch of abstract harmonic analysis on (abelian) locally compact groups, which generalizes Fourier transformation on reals and circle. The main point about sin and cos (or rather complex exponent $e^{i n x}$) is that it is a character (continuous homomorphism from a group to a circle) and it is not hard to see that those are the only characters of the circle. That what makes Fourier transform so powerful. If you generalize it along the direction which drops characters, you'll probably get a much weaker theory.
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's a standard series computation to show that
$$ \sum_{n \ge 1} \frac{x^n}{n} = \log \frac{1}{1 - x} $$
Now substitute $x = e^{i t}$ and take the real part.
(As an aside, the reason I write the identity this way is that this is the version which is combinatorially significant.)