Fourier decomposition is a well-established technique used in the study of stellar pulsation. However, the quality of reconstructed light curves using this method is reduced when the observed data have uneven phase coverage. We use simulated annealing techniques together with Fourier decomposition to improve the quality of the Fourier decomposition for many Optical Gravitational Lensing Experiment LMC fundamental-mode Cepheids. This method restricts the range that Fourier amplitudes can take. The ranges are specified by well-sampled Cepheids in the Galaxy and Magellanic Clouds. We also apply this method to reconstructing Cepheid light curves observed by the Hubble Space Telescope (HST). These typically consist of 12 V-band and four I-band points. We employ a direct Fourier fit to the 12 V-band points using the simulated annealing method mentioned above and explicitly derive and use Fourier interrelations to reconstruct the I-band light curve. We discuss advantages and drawbacks of this method when applied to HST Cepheid data over existing template methods. Application of this method to reconstruct the light curves of Cepheids observed in NGC 4258 shows that the derived Cepheid distance (mu(0) = 29.38 +/- 0.06 mag, random error) is consistent with its geometrical distance (mu(0) = 29.28 +/- 0.09 mag) derived from observations of its water maser