Pulsars, a small portion of celestial sources that emit radiation at varying intensity, provide new possible navigation algorithms which are different from steady point sources. Time-delay estimation is one of the key aspects of pulsar-based navigation technology. Previous work for pulse phase estimation uses a maximum likelihood estimator (MLE) for the phase-in time domain, which is seen as one of the most useful phase estimators. However, the analytic solution for phase cannot be found using MLE. As a result, a brute-force method involving nested, iterative grid-searches is needed to solve this MLE issue, which leads to lots of computations. In order to solve this problem, a multiple harmonic X-ray pulsar signal phase estimation (MHSPE) method is proposed. This method uses the times of arrivals (TOAs) measured pulsar signal to estimate the time-delay in the frequency domain. In this paper, firstly we use the arrival time to derive the maximum-likehood (ML) estimation of phase-delay by fundamental frequency, then an analytic expression for the fundamental frequency phase is obtained. The MHSPE method based on the fundamental frequency phase equation, calculates different harmonic phases by generalizing the analytic expression of fundamental frequency phase, and the normalized amplitude of each harmonic in the spectrum is used as the weight of each harmonic phase. Finally, the weighted average of harmonic phases, which is given by the final analytic expression, is used as the estimation of the moment. To evaluate the MHSPE method, the error and variance equations are calculated and the MHSPE method is demonstrated to be unbiased and consistent. Moreover, by comparing with the ML estimation of the first harmonic, if the amplitudes of harmonic in the spectrum are almost the same, the signal-to-noise ratio (SNR) of MHSPE improves m1/2 times when the number of harmonic waves is m. Compared with the Cramer-Rao bound of pulsar time-delay estimation, the derivative of pulsar signal in the time domain reflects the number of harmonic waves in the frequency domain. Hence, the MHSPE can greatly approximate to the Cramer-Rao bound for the estimation of pulsar signal timedelay when the harmonic number increases. Finally, we utilize the TOAs of Crab pulsar, observed by Rossi X-ray timing explorer (RXTE) spacecraft, to verify the performance of MHSPE. The results show that for low SNR of pulsar signal, MHSPE can obtain high precision phase estimations. When the amplitude of the harmonic in the spectrum is larger, the estimation variance of the harmonic phase tends to be smaller. The projection orbit determined by MHSPE method can match the projection of RXTE in the direction of Crab pulsar, with the observed time increasing, the estimation accuracy converges rapidly to Cramer-Rao bound.