Spectral analysis of sinusoidal signals is a classical but still open problem in statistical signal processing, finding its applications in a wide range of areas. This problem consists of two parts - sinusoidal model order detection and parameter estimation. During the recent decades, the problem of analyzing the sinusoidal signals from multiple channels, which are contaminated by different undesired harmonics, has attracted considerable attention. Given the corresponding observations, the goal is to determine the unknown orders and parameters of the sinusoidal signals in the multiple channels, after which the signal parametrization is complete. This problem is of great research value not just because it is interesting and practical, but also there are two main significant advantages compared with single-channel modeling: • The multi-channel setup means more observed data and the parameter estimation refinement of the common mode sinusoidal components is expected, which makes it feasible to extract the common information in a more accurate way. • In a multi-source scenario, if the sources are overlapping in one dimension, a single-channel setup will not be able to resolve the sources. On the other hand, this issue can be alleviated with the multi-channel setup by considering joint higher dimensional modeling. In addition, decimation technique is utilized in the parameter estimation of oversampled multiple complex sinusoids for the sake of lower computational complexity and higher estimation resolution, where the decimative signals belong to a special form of multi-channel sinusoidal signals with the same amplitudes and frequencies. And accurate parameter estimation for dual-channel sinusoidal signal is extensively useful in the electronic measurement. Such applications are also the motivations of the research on spectral analysis of multi-channel sinusoidal signals. In this thesis, the multi-channel sinusoidal modeling consists of four parts, that is, oversampling parameter estimation for multiple sinusoids; accurate dual-channel sinewave parameter estimation; parametric modeling for damped sinusoids from multiple channels; and spatial-temporal modeling for harmonic signal from microphone array. In the oversampling parameter estimation for multiple complex sinusoids, the parameters of continuous-time frequencies are of interest. In signal processing, oversampling technique is the process of sampling a signal with a sampling rate significantly higher than the Nyquist frequency of the signal being sampled. Oversampling is utilized to obtain more data in a fixed duration, and is expected to improve the estimation accuracy. Nevertheless, two problems occur in spectral estimation, that is the problems of smaller frequency separation and higher computational complexity. To alleviate these problems, the oversampling weighted least squares frequency estimator with decimation is proposed. For the problem of sinusoidal parameter estimation at two channels, the parameters of common frequency, amplitudes, initial phases and possibly DC offsets, are of interest. Under the assumption of white Gaussian noise, an iterative linear leastsquares algorithm for accurate frequency estimation is devised. The remaining parameters are then determined according to linear least-squares fit with the use of the frequency estimate. The parameter of phase-difference is another key quantity. To estimate it, two algorithms have been proposed. The first one utilizes the maximum likelihood criterion to find the initial phases of dual-channel outputs, respectively, and the phase-difference estimate is then given by their difference. Algorithm extension to unknown frequency and/or noise powers is also studied. The development of the second method is based on the weighted linear prediction approach with a properly chosen sampling frequency. On parametric modeling of damped sinusoidal signals from multiple channels, it is aimed at addressing the issues of their model order detection and parameter estimation from a new and complete viewpoint via performing the parametric modeling with joint model selection and parameter estimation. It consists of three parts. Firstly, we extend the subspace-based automatic model order selection method to the multi-channel scenario, and detect the number of the distinct sinusoidal poles in the multiple channels with the multi-channel model order estimator. Secondly, we extend the iterative quadratic maximum likelihood approach to the current problem, that is, parameter estimation for the sinusoidal poles from multiple channels, which is referred to as the multi-channel iterative quadratic maximum likelihood estimator. Thirdly, sinusoidal model selection, or matching the estimated poles to their corresponding channels, is realized based on a sequence of hypothesis tests. At each test, we compute the significance of the maximum correlation between the estimation residual and a sinusoidal function, whose statistical property is derived from the extreme value theory about the distribution of the maximum of stochastic fields. We refer this scheme to as extreme value theory selector. The problem of spatial-temporal modeling for harmonic signal from microphone array is solved from two aspects. Firstly, we propose to estimate the fundamental frequency and direction of arrival (DOA) in two stages. At first, the multi-channel optimally weighted harmonic multiple signal classification estimator is devised, and the estimation of fundamental frequency is conducted. Then we make use of the spatial-temporal multiple signal classification estimator to estimate the DOAs with the estimated fundamental frequencies. Although the two-stage method is more computationally efficient, it cannot resolve the sources with overlapping frequencies or DOAs. To overcome this problem, in the second part, we perform DOA and fundamental frequency estimation in a joint way. In practice, there also occur the problems of order detection and detection of missing harmonics of each source. To take this issue into account in harmonic modeling, we propose to perform joint estimation based on optimal filtering method and with the maximum harmonic model, and then the model selection is accomplished according to the powers of the respective harmonic components of each source, and the maximum a posteriori criterion.