My impression is that for most choices of H this is a hard question. In general, the largest possible cardinality of a set of independent Poisson commuting functions on a 2n-dimensional symplectic manifold (where by independent I mean that their differentials are linearly independent at each point in their domain--in particular H and f(H) are not independent) is n. Sketch proof: Given k independent Poisson commuting functions, their Hamiltonian vector fields span a k-dimensional isotropic subspace (i.e., the symplectic form dxdp vanishes on the subspace), and it's an easy linear algebra exercise to show that the largest possible dimension of an isotropic subspace is n.
If there is an independent set of n functions, including H, which mutually Poisson commute, then the Hamiltonian system associated to H is called integrable, and can be solved exactly using action-angle coordinates--see Arnold's book Mathematical Methods of Classical Mechanics. The example you gave is integrable, but in general this is a fairly-rarely satisfied condition.
There are some obstructions to H having independent functions that Poisson commute with it. For example, Poincare noticed that if the Hamiltonian vector field of H has a periodic orbit, so if the orbit has period T the time-T map F of the flow has some fixed point y, then looking at the derivative of F at y gives the following restriction: if there are k Poisson-commuting functions (including H) which are independent along the orbit through y, then the derivative of F at y has to have the eigenvalue 1 with multiplicity at least 2k. The generic situation is that 1 occurs as an eigenvalue with multiplicity just 2, so in this sense the typical Hamiltonian doesn't have any other functions that Poisson commute with it along its periodic orbits.
Another point that occurs to me is that if almost every level set of H contains a dense orbit for the Hamiltonian flow of H, then since any K such that {K,H}=0 is constant along the Hamiltonian flow of H it would follow that any such K would be constant on almost every level set of H. But then K would be constant on every level set of H just by continuity, and this is only possible if K has the form f(H) that you mentioned. However I'm not sure if the assumption that almost every level set of H has a dense orbit is realistic--in particular I can't think of any examples where this holds on TnxRn.
One positive statement that can be made is the following. If you look at the proof of Darboux's theorem that's given in Arnold's above mentioned book, you'll see that it proves the following slightly stronger statement: if y is any point at which dH is nonvanishing, then there are local coordinates (Hi,Ki) around in terms of which the symplectic form is given by \sum dHidKi and H1=H. So (H=H1,H2,...,Hn) is an n-tuple of Poisson commuting functions, just defined near the given point y. So locally there are certainly functions that Poisson commute with H--they just can't generally be extended globally.