How did Heisenberg come up with the Canonical Commutation relation ($\hat X \hat P-\hat P\hat X=i\hbar$)? All answers to questions like this dodge the question by saying it's a postulate of Matrix Mechanics, so let me rephrase it. Instead of how to derive the CCR, how does it follow from Heisenberg's matrices:
$$\hat X=\{X_{mn}e^{i\omega_{mn}t}\}\quad\hat P=\{P_{mn}e^{i\omega_{mn}t}\}$$
Where $\omega_{mn}=\frac{E_m-E_n}{\hbar}$, such that
$$\hat X(t)=\left(\begin{matrix}
X_{11}&X_{12}e^{i\omega_{12}t}&X_{13}e^{i\omega_{13}t}&\cdots\\
X_{21}e^{i\omega_{21}t}&X_{22}&X_{23}e^{i\omega_{23}t}&\cdots\\
X_{31}e^{i\omega_{31}t}&X_{32}e^{i\omega_{32}t}&X_{33}\cdots\\
\vdots&\vdots&\vdots&\ddots
\end{matrix}\right)$$
$$\hat P(t)=\left(\begin{matrix}
P_{11}&P_{12}e^{i\omega_{12}t}&P_{13}e^{i\omega_{13}t}&\cdots\\
P_{21}e^{i\omega_{21}t}&P_{22}&P_{23}e^{i\omega_{23}t}&\cdots\\
P_{31}e^{i\omega_{31}t}&P_{32}e^{i\omega_{32}t}&P_{33}\cdots\\
\vdots&\vdots&\vdots&\ddots
\end{matrix}\right)$$
And from the Old quantum condition:
$$\oint pdx=nh$$
I know Heisenberg came up with it through analogies between the matrices and the classical observables, but I'd like to know how he did it.
 A: I see you properly posted your question and got your answer in the History of science SE. To the extent you have not deleted the question from here, yet, I gather you are asking for help with the evident logic detailed in the WP article. Any good matrix mechanics text, including Herbert Green's superb and elegant book cited in WP will walk you through it. Here, I'll just help you recognize the pattern, you almost certainly already know from the Fock space of the oscillator. We are certainly not going to solve the full sudoku riddle Born's intuitive genius guessed the answer to, to crucially advance our civilization! (He "convinced [himself] that the off-diagonal matrix elements had to vanish"...)
You already saw that the old Bohr condition enforces periodic motions on a circle, and thus the discrete Fourier mode structure you wrote down, with, significantly, 
$\omega_{mn}=-\omega_{nm}$ and real coefficients $X_{mn}$ and $P_{mn}$. As a result, reflection across the diagonal will yield the opposite phase in the exponential of the mode, and matrix multiplication will preserve the structure, i.e. product matrices will have the same phases in a given matrix location as each  factor did. (I'm sure there is a respectable mathematical term for this, but I don't know it. It is a feature of the unitary similarity transformation with respect to the time-evolution operators.) 
For the linear spectrum of the oscillator then, all these phases are integer multiples of the fundamental frequency ω, and energy is "conserved" in every matrix multiplication. Integrating over the cycle in Bohr's condition kills any and all time-dependent harmonics and projects onto the zero-frequency subspace.
You now know (e.g. from Messiah's vI, XII, §5, or your favorite text) the Fock representation diagonal in the number operator, $a^\dagger a$=diag(0,1,2,3,..), where 
$$
a= 
 \left(\begin{matrix}
0&\sqrt{1}&0&\cdots\\
 0& 0& \sqrt{2}&\cdots\\
 0 &0&0  &\cdots\\
\vdots&\vdots&\vdots&\ddots
\end{matrix}\right)$$  such that, directly, together with its Hermitian conjugate, they dictate
$$\hat X(t)=\sqrt{\frac{\hbar}{m\omega}} \left(\begin{matrix}
0&\sqrt{1} e^{i\omega t}&0&\cdots\\
\sqrt{1} e^{-i\omega t}&0&\sqrt{2} e^{i\omega t}&\cdots\\
0&\sqrt{2} e^{-i\omega t}&0\cdots\\
\vdots&\vdots&\vdots&\ddots
\end{matrix}\right)$$
the time-evolute of (5.37) in Merzbacher's book, and, likewise, 
$$\hat P(t)=-im\omega \sqrt{\frac{\hbar}{2m\omega}} \left(\begin{matrix}
0&\sqrt{1} e^{i\omega t}&0&\cdots\\
-\sqrt{1} e^{-i\omega t}&0&\sqrt{2} e^{i\omega t}&\cdots\\
0&-\sqrt{2} e^{-i\omega t}&0\cdots\\
\vdots&\vdots&\vdots&\ddots
\end{matrix}\right),$$
both Hermitian. The quantum phase of the ground state has washed out in the overall similarity transformation and I have flipped the sign of the common frequency, for simplicity, immaterially for our purpose here.
Now, if you liked complication, you might multiply your original matrices in two different orders, but, to illustrate the point, let us do this with the oscillator ones, here,
$$\hat X(t)\hat P(t) =-i  \frac{\hbar}{2 }  \left(\begin{matrix}
0&\sqrt{1} e^{i\omega t}&0&\cdots\\
\sqrt{1} e^{-i\omega t}&0&\sqrt{2} e^{i\omega t}&\cdots\\
0&\sqrt{2} e^{-i\omega t}&0\cdots\\
\vdots&\vdots&\vdots&\ddots
\end{matrix}\right) \cdot  \left(\begin{matrix}
0&\sqrt{1} e^{i\omega t}&0&\cdots\\
-\sqrt{1} e^{-i\omega t}&0&\sqrt{2} e^{i\omega t}&\cdots\\
0&-\sqrt{2} e^{-i\omega t}&0\cdots\\
\vdots&\vdots&\vdots&\ddots
\end{matrix}\right)\\
 =-i  \frac{\hbar}{2 }   \left(\begin{matrix}
-1&0&\sqrt{2} e^{i2\omega t}&\cdots\\
0&-1& 0&\cdots\\
-\sqrt{2} e^{-i2\omega t}& 0&-1\cdots\\
\vdots&\vdots&\vdots&\ddots
\end{matrix}\right).$$
Here, contributions from outside the 3×3 subspace have been truncated, as they will not behave differently, below, than the ones with indices 1,2,3, kept.
Subtracting the Hermitian conjugate from this nets the celebrated 
$$\hat X(t)\hat P(t)  - \hat P(t) \hat X(t) =i \hbar \operatorname{diag} (1,1,1,...),$$
where you are invited to note the automatic cancellation of all the non-diagonal (time-dependent) terms!
This is illustrated here for the 3×3 subspace of the ground and first two excited states you are depicting explicitly, but the systematics should be evident. The structure was familiar to the best mathematical physicists of that generation who handled normal modes in vibrating drums, etc... The Hilbert space techniques for this were being developed "next door", at Hilbert's math department, where Jordan had already had his apprenticeship with Courant.
