concentration of functions of Gaussian processes

Question

Let $\mathcal{C}\in\mathbb{R}^n$ be a subset of the unit ball. Also let $\mathbf{a}_1,\mathbf{a}_2,\ldots,\mathbf{a}_m\in\mathbb{R}^n$ be i.i.d. random Gaussian vectors $\mathcal{N}(\mathbf{0},\mathbf{I})$ then by Gordon's lemma we know \begin{align*} \underset{\mathbf{y}\in\mathcal{C}}{\sup}\bigg|\frac{1}{m}\sum_{k=1}^m (\mathbf{a}_k^T\mathbf{y})^2-\|\mathbf{y}\|_{\ell_2}^2\bigg|\le \delta, \end{align*} holds with high probability as long as \begin{align*} m\ge c\frac{\omega^2(\mathcal{C})}{\delta^2}, \end{align*} for a fixed numerical constant $c$ (in fact I think $c$ is at most $9$). Here $\omega(\mathcal{C})$ is the Gaussian width defined as \begin{align*} \omega(\mathcal{C})=\mathbb{E}\big[\sup_{\mathbf{y}\in\mathcal{C}}\mathbf{g}^T\mathbf{y}\big], \end{align*} with $\mathbf{g}$ is a random Gaussian vector $\mathcal{N}(\mathbf{0},\mathbf{I})$. With this introduction here is my question. I want a similar result to hold about \begin{align*} \underset{\mathbf{y}\in\mathcal{C}}{\sup}\bigg|\frac{1}{m}\sum_{k=1}^m f(\mathbf{a}_k^T\mathbf{y})-E_{\mathbf{g}}[f(\mathbf{g}^T\mathbf{y})]\bigg|\le \delta, \end{align*} with $\mathbf{g}$ is a random Gaussian vector $\mathcal{N}(\mathbf{0},\mathbf{I})$. The only change is replacing $f(x)=x^2$ in Gordon's result with a general $f$. Obviously this is not true in general but assume that the function $f$ is bounded or its Lipschitz is the results above still true with with high probability as long as \begin{align*} m\ge c\frac{\omega^2(\mathcal{C})}{\delta^2}, \end{align*} with $c$ a constant that only depends on properties of the function $f$ e.g. bound, Lipschitz constant etc.