H-FISTA (Hierarchical Fast Iterative Shrinkage Thresholding Algorithm) retrieves the phases of the wavefield from intensity measurements for pulsar spectroscopy. The code accepts input data in ASCII format as produced by PSRchive's (ascl:1105.014) psrflux function, a FITS file, or a pickle. If using a notebook, any custom reader can be used as long as the data ends up in a NumPy array. H-FISTA obtains sparse models of the wavefield in a hierarchical approach with progressively increasing depth. Once the tail of the noise distribution is reached, the hierarchy terminates with a final unregularized optimization, resulting in a fully dense model of the complex wavefield that permits the discovery of faint signals by appropriate averaging.