Structural health monitoring assesses structural integrity by processing the measured responses of structures. One particular group in the structural health monitoring research is to conduct the operational modal analysis and then to extract the dynamic characteristics of structures from vibrational responses. These characteristics include natural frequencies, damping ratios, and mode shapes. Deviations in these characteristics represent the changes in structural properties and also imply possible damage to structures. In this study, a new stochastic system identification is developed using multivariate time-frequency distributions. These time-frequency distributions are derived from the short-time Fourier transform and subsequently yield a time-frequency matrix by stacking them with respect to time. As the derivation in the data-driven stochastic subspace system identification, the future time-frequency matrix is projected onto the past time-frequency matrix. By exploiting the singular value decomposition, the system and measurement matrices of a stochastic state-space representation are derived. Consequently, the dynamic characteristics of a structure are obtained. As compared to the time-domain stochastic subspace system identification, the proposed method utilizes the past and future matrices with a lower dimension in projection. A spectral magnitude envelope can be applied to the time-frequency matrix to highlight the major frequency components as well as to eliminate the components with less influence. To validate the proposed method, a numerical example is developed. This method is also applied to experimental data in order to evaluate its effectiveness. As a result, performance of the proposed method is superior to the time-domain stochastic subspace system identification.