There are two ways to solve this problem - one by ergodic methods, and the other one using purely harmonic methods.

The harmonic method you are indicating is just to take the delta function of the point i (I'm looking at the locally-symmetric space, you can easily translate to the homogenuous space situation).

Expand delta in $L^{2}(\Gamma \backslash H)$ (one need to be a bit more careful about what that means).
Now, your measures (averging over $1/N$) and the measures which are achived as a push-forward of the delta by the Hecke operators are closely related (they differ by some $o_{f}(1)$ for any automorphic Schwartz function $f$), hence if you can proved equidistribution of those translates, you can prove equidistribution of the averging over the $1/N$ cycles.

Now if you work only harmonically, you will want to use Weyl's equidistribution criterion.
As you know, we can take the Hecke-Mass forms as a basis to the (cuspical part of) space.
In this point view, it is cleat that any non-trivial bound towards the Ramanujan conjecture, will give you equidistribution (at least in the cuspidal part of the spectrum), and this proof is effective, any bound will translate into an explicit rate.

You would be right to remind me that there is also a continuous spectrum, but the truth is that the Eisenstein series computation is easer than the Maass forms part.

The explicit computation appears in Ullmo's article here - http://www.math.u-psud.fr/~ullmo/Publications/coursMontrealfinal.pdf (see 2.3), or the general article by Clozel-Ullmo-Oh, or in Ullmo's article here http://www.math.u-psud.fr/~ullmo/Publications/clozel-ullmo.pdf (French).

There's another way to prove this, more ergodic theoretical, which basically uses Margulis' mixing trick, but in print it appear as an old result due to Sarnak (he disguised it by some Eisenstein series calculation).
There you average over the circles, and then uses mixing of the geodesic flow to get the required result. Such a result would force a bound towards Ramanujan, and this is basically what Venkatesh do in his subconvexity paper - see 3.2 in Venkatesh's "Sparse equidistribution problems" or the closely realted explanation 1.2 in Michel-Venkatesh.
In this way, bounds towards Ramanujan are simply way to effectivaze mixing rates.

P.S. a non-effective "easy way" to prove this equidistribution statment (i.e. without spectral theory) would be to use Ratner's theorems.