  Pulsar timing arrays (PTAs) are searching for gravitational waves from
supermassive black hole binaries (SMBHBs). Here we show how future PTAs could
use a detection of gravitational waves from individually resolved SMBHB sources
to produce a purely gravitational wave-based measurement of the Hubble
constant. This is achieved by measuring two separate distances to the same
source from the gravitational wave signal in the timing residual: the
luminosity distance $D_L$ through frequency evolution effects, and the parallax
distance $D_\mathrm{par}$ through wavefront curvature (Fresnel) effects. We
present a generalized timing residual model including these effects in an
expanding universe. Of these two distances, $D_\mathrm{par}$ is challenging to
measure due to the pulsar distance wrapping problem, a degeneracy in the
Earth-pulsar distance and gravitational wave source parameters that requires
highly precise, sub-parsec level, pulsar distance measurements to overcome.
However, in this paper we demonstrate that combining the knowledge of two SMBHB
sources in the timing residual largely removes the wrapping cycle degeneracy.
Two sources simultaneously calibrate the PTA by identifying the distances to
the pulsars, which is useful in its own right, and allow recovery of the source
luminosity and parallax distances which results in a measurement of the Hubble
constant. We find that, with optimistic PTAs in the era of the Square Kilometer
Array, two SMBHB sources within a few hundred Mpc could be used to measure the
Hubble constant with a relative uncertainty on the order of 10 per cent.
