In this paper, we introduce the \textbf{G}eneralized \textbf{L}inear \textbf{S}pectral \textbf{S}tatistics (GLSS) of a high-dimensional sample covariance matrix $\bm{S}_n$, denoted as $\operatorname{tr}f(\bm{S}_n)\bm{B}_n$, which effectively captures distinct spectral properties of $\bm{S}_n$ by incorporating an ancillary matrix $\bm{B}_n$ and a test function $f$. The joint asymptotic normality of GLSS associated with different test functions is established under mild assumptions on $\bm{B}_n$ and the underlying distribution, when the dimension $n$ and sample size $N$ are comparable. The convergence rate of GLSS is determined by $\sqrt{{N}/{\operatorname{rank}(\bm{B}_n)}}$. Subsequently, we propose a novel functional projection approach based on GLSS for hypothesis testing on eigenspaces of ``population-spiked'' covariance matrices, showcasing a universality phenomenon in the magnitude of the spikes. The theoretical accuracy of our results established for GLSS and the advantages of the newly suggested testing procedure are demonstrated through various numerical studies.