This paper addresses the identification of models using multiple sinusoidal. forcing subject to practical difficulties such as unknown initial states, offsets, slow and periodic disturbances, noise, and unknown model structures. A linear regression equation is derived by integrating the system equation excited by a single sinusoid and extended to the case of concurrent multiple sinusoids. The regression equation can be used to estimate the model parameters including time delay in a least-squares fashion. A simple scheme based on the condition number of the matrix formed by the regression vector is presented to infer the best model order. Two-stage sinusoidal testing, represented as two sets of sinusoids applied sequentially, is then developed to generate response data that are informative enough to deal with the aforementioned identification difficulties. This requires only the application of the regression equation for concurrent multiple sinusoids modified in a clever and sequential manner.