Terahertz (THz) radiation is the electromagnetic spectrum lying between the microwave and infrared portions, and covers wavelength between 3mm and 30um. Terahertz spectroscopy has been used in qualitative and quantitative analysis since molecular vibrations and rotations of a great large number of materials lie in this region. Many techniques, such as principal component regression (PCR) and partial least square regression (PLSR), have been applied in quantitative analysis of terahertz time-domain spectroscopy. However, those “soft-modeling” methods, though effective in many cases, all require a large number of samples and have poor extrapolation capability, and the utilization is also restricted in some cases (samples lack or requirement of extrapolation). A parametric method of spectrum, Indirect Hard Modeling (IHM), was proposed to resolve terahertz spectral matrix and quantify the mixtures. Savitzky-Golay method was used to eliminate the influence of noise. Moreover, the features of terahertz spectra were considered and the artificial baseline was also resolved as a signal spectrum to improve spectrum extraction accuracy. Subsequently, the analysis of IHM was carried out, and the ability of prediction and accuracy of the model in the case of low training samples were discussed. In order to validate the capability of the proposed method, quaternary systems of Rifampicin, Lactose monohydrate, microcrystalline cellulose and magnesium stearate were tested. The regression correlation coefficient (R) and root mean square error (RMSE) were employed to estimate the performance of the model. Finally, the comparison of IHM and PLSR were carried out, and theoretical analysis and experimental results indicated that the number of samples required by IHM modeling could be reduce to two during calibration, at the same time, the accuracy of quantitative and extrapolation capabilities were also improved.