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-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-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 Kilometre Array, two fortuitous 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.