Harmonic retrieval has been one of the fundamental problems across a diverse range of important fields, from source localization, radar and sonar signal processing, channel modeling in wireless communication systems, to medical areas. Generally speaking, harmonic retrieval algorithms can be divided into parametric and non-parametric approaches. Parametric based algorithms assume that the received signal follows a specific model. These kinds of estimators provide a superior estimation performance in terms of accuracy, computational efficiency, resolution and/or identifiability, but will fail when the assumption of the data model is incorrect. Non-parametric based algorithms, however, can generate the estimated spectrum without any assumption of a data model. In this thesis, several algorithms are devised and analyzed in terms of the harmonic retrieval problem, from 1-dimensional (1-D) to R-D where R ≥ 2. Synonyms of harmonic retrieval include frequency estimation, parameter estimation and spectrum/spectral analysis. The commonly used algorithms for 1-D harmonic retrieval of multiple sinusoids can be divided into two groups: subspace based and iterative approaches. The subspace based technique is usually computationally efficient but the accuracy is not very satisfying. Iterative algorithms, however, can achieve a high level of accuracy when the signal to noise ratio (SNR) is sufficiently high, but will be costly when the data size is large. This thesis introduces a new approach, the principal-singular-vector utilization for model analysis (PUMA), which combines the advantages of both subspace and iterative techniques, for the frequency estimation of complex-valued and real-valued sinusoids. The main idea is to use a matrix without repeated elements to represent the observed signal and apply the weighted least squares (WLS) on its signal subspace to retrieve parameters accurately. It is proved that for small error conditions, the frequency estimate is approximately unbiased. Furthermore, its variance equals to Cramér-Rao lower bound (CRLB) for the single complex cisoid problem. Although the 1-D problem is the most common, the R-D harmonic retrieval with R ≥ 2 has recently been receiving increased interests. Conventional R-D parameter estimation techniques include the maximum likelihood (ML) based and subspace approaches. The ML based methods are only feasible when the problem is 2-D due to their extremely high computational requirements, while the accuracy and/or computational efficiency of subspace methodology remains unsatisfactory. Furthermore, the subspace approaches can be divided into matrix based and tensor based ones. The latter one usually give a more accurate subspace estimation with a higher complexity. In this thesis, the tensor-PUMA (TPUMA) approach for R-D parameter estimation of a summation of F ≥ 2 sinusoids in additive white Gaussian noise is devised. Firstly, with the use of an iterative procedure, the sinusoidal parameters at one dimension are accurately estimated from its signal subspace, which is retrieved using a decomposition method. The parameters in the remaining dimensions are then solved such that the R-D harmonics are automatically paired-up. Algorithm modification for a single R-D tone is made and it is shown that the frequency estimates are asymptotically unbiased while their variances approach CRLB at sufficiently high SNR conditions. It is also shown that the TPUMA algorithms are computationally simpler than conventional parameter estimators. In the special case of R-D harmonic retrieval, the single cisoid in additive white Gaussian noise is studied, and two proposals are then developed. By exploiting the correlation of the data samples, R single-tone sequences, which contain the R-D frequencies, are constructed, and accurate frequency estimation is attained from each constructed sequence. The two devised estimators are proved to be asymptotically unbiased while their variances achieve CRLB when the SNR and/or data length tend towards infinity, and also shown to be more computationally attractive than the conventional schemes. Sometimes identical frequencies exist along each dimension. When this occurs, the TPUMA scheme cannot work. To overcome this problem, a tensor eigenvector-based (TEV) estimator is proposed. Our underlying idea is to utilize the tensorial structure of the received data and then apply a higher-order singular value decomposition (HOSVD) to obtain the tensor-based signal subspace. It is decomposed into a set of single-tone tensors from which single-tone vectors can be constructed by another HOSVD later on. In doing so, the R-D multiple sinusoids are converted to a set of single-tone sequences whose frequencies are individually estimated according to structure least squares (SLS). The bias and variance of the subspace-based method are also analyzed and it is proved that the frequency estimates are approximately unbiased. It is shown that the proposed algorithm can provide superior estimation performance. Recently, an iterative adaptive approach (IAA), which employs the WLS technique to update the whole spectrum, is proposed for high-resolution spectral estimation. However, the IAA algorithm is computationally expensive for 2-D spectrum imaging. To address this issue, a new fast implementation for IAA utilizing subspace technique is devised. The main idea is to utilize the vectors of the subspaces to retrieve the spectrum. It is observed that the whole spectrum can be written as a summation of several sub-spectra, and each sub-spectrum can be updated in an iterated 1-D format based on the dominant singular vectors instead of 2-D format based on the original data samples. Fast computation is thus achieved. Furthermore, the number of signals, F, can be known or unknown for this procedure, which makes this algorithm superior to the algorithms where information of F is required. It is shown that the proposed fast IAA can achieve a performance similar to IAA with less computational complexity.