I don't have Datta now available, but if I remember correctly, he showed how to derive the "Bloch" Hamiltonian for general systems. If you understand what he does, you should be able to reproduce the same calculation for graphene. You know how to derive the bandstructure for a unit cell of two carbon atoms, right? Exactly the same principle holds if you choose a larger unit cell. I cannot give you any reference, since I doubt that anyone has published or showed such calculations just for fun... I understand that for example for systems with defects, it may be necessary to introduce a larger unit cell, but for pristine graphene, it does not really make any sense.
But to make sure that we understand what zone-folding means, let us briefly study the simplest example: a 1D TB chain with N atoms and periodic boundary conditions. The discrete Schrödinger equation for each site is
<br />
-t\psi_{j-1}-t\psi_{j+1} = E \psi_j,<br />
which can be solved using the Bloch ansatz \psi_j=ue^{ikj}, leading to
<br />
E_k=-t(e^{-ik}+e^{ik})=-2t\cos(k).<br />
The distinct k-values are k_n=2\pi n/N, \ n=0,1,\dots,N-1. The Bloch Hamiltonian is a 1*1 matrix
<br />
H_{eff}(k)=-t(e^{-ik}+e^{ik}),<br />
since from any site one can hop to the right (exponential factor eik from the Bloch Ansatz) or to the left (exponential factor e-ik from the Bloch Ansatz).
Now suppose that we use a unit cell of two atoms. We know that this should lead to a identical bandstructure, but with "folded zones". The number of unit cells is then N/2, and our Ansatz would be
<br />
\left(<br />
\begin{array}{c}<br />
\psi_j \\<br />
\psi_{j+1} \\<br />
\end{array}<br />
\right)<br />
=<br />
\left(<br />
\begin{array}{c}<br />
u_1 \\<br />
u_2\\<br />
\end{array}<br />
\right)<br />
e^{ik(j+1)/2}, \ j=1,3,5,\dots,N-1,<br />
where I chose the exponential factor simply so that the first unit cell gets factor eik, the second e2ik etc. The distinct values of k are now (assuming N is even for this to make sense) k_n=4\pi n/N,\ n=0,1,\dots,N/2-1. As expected, the number of allowed k-values in the Brillouin zone was halved due to the bigger unit cell. This leads to the Schrödinger equation
<br />
\left(<br />
\begin{array}{cc}<br />
0 & -t(1+e^{ik}) \\<br />
-t(1+e^{-ik}) & 0 \\<br />
\end{array}<br />
\right)<br />
\left(<br />
\begin{array}{c}<br />
u_1 \\<br />
u_2\\<br />
\end{array}<br />
\right)<br />
=<br />
E<br />
\left(<br />
\begin{array}{c}<br />
u_1 \\<br />
u_2\\<br />
\end{array}<br />
\right).<br />
Diagonalization gives the eigenvalues
<br />
E_k=\pm t \sqrt{(1+\cos(k))^2+\sin(k)^2}=\pm t\sqrt{2+2\cos(k)}=\pm 2t\cos(k/2).<br />
(Here one should be a little careful with the trigonometry, I guess.) As expected, we got two branches that meet each other at k=pi. Rescaling the new k-variable by a factor of two, it would be clear how the original bandstructure Ek=-2tcos(k) changes to Ek=pm 2tcos(k): plotting shows that the old bandstructure in the Brillouin zone (-pi,pi) got folded with respect to the points -pi/2 and pi/2. This scaling also halves the Brillouin zone to (-pi/2,pi/2).
Now, one can directly see how this can be generalized to more complicated systems: the only problem is that in higher dimensions and in more complicated lattices, the intuition of folding is easily lost. But in principle, the analysis should be straightforward. Just try to do it yourself and ask if you run into any problems.