Classical second-order spectral analysis, which is based on the Fourier transform of the autocovariance functions, focuses on summarizing the oscillatory behaviors of a time series. However, this type of analysis is subject to two major limitations: first, being covariance-based, it cannot captures oscillatory information beyond the second moment, such as time-irreversibility and kurtosis, and cannot accommodate heavy-tail dependence and infinite variance; second, focusing on a single time series, it is unable to quantify the association between multiple time series and other covariates of interests. In this article, we propose a novel nonparametric approach to the spectral analysis of multiple time series and the associated covariates. The procedure is based on the copula spectral density kernel, which inherits the robustness properties of quantile regression and does not require any distributional assumptions such as the existence of finite moments. Copula spectral density kernels of different pairs are modeled jointly as a matrix to allow flexible smoothing. Through a tensor-product spline model of Cholesky components of the conditional copula spectral density matrix, the approach provides flexible nonparametric estimates of the copula spectral density matrix as nonparametric functions of frequency and covariate while preserving geometric constraints. Empirical performance is evaluated in simulation studies and illustrated through an analysis of stride interval time series.