We study the solution $h_\varepsilon$ of the Kardar-Parisi-Zhang (KPZ) equation for $d \geq 3$: $$ \frac{\partial}{\partial t} h_{\varepsilon} = \frac12 Δh_{\varepsilon} + \bigg[\frac12 |\nabla h_\varepsilon |^2 - C_\varepsilon\bigg]+ β\varepsilon^{\frac{d-2}2} ξ_{\varepsilon} $$ with $h_\varepsilon(0,x)=0$. Here $ξ_\varepsilon=ξ\star φ_\varepsilon$ is a spatially smoothened (at scale $\varepsilon$) Gaussian space-time white noise and $C_\varepsilon$ is a divergent constant as $\varepsilon\to 0$. When the disorder $β$ is sufficiently small and $\varepsilon\to 0$, $h_\varepsilon(t,x)- \mathfrak h^{\mathrm{st}}_{\varepsilon}(t,x)\to 0$ in probability where $\mathfrak h^{\mathrm{st}}_{\varepsilon}(t,x)$ is the {\emph stationary solution} of the KPZ equation - more precisely, $\mathfrak h^{\mathrm{st}}_{\varepsilon}$solves the above equation with a random initial condition (that is independent of the driving noise $ξ$) and its law is constant in $(\varepsilon,t,x)$. In the present article we quantify the rate of the above convergence in this regime and show that the fluctuation {\emph about} the stationary solution $$ (\varepsilon^{1-\frac d2} [h_\varepsilon(t,x) - \mathfrak h^{\mathrm{st}}_{\varepsilon}(t,x)])_{x,t} $$ converges pointwise (with finite dimensional distributions in space and time) to a Gaussian free field (GFF) evolved by the deterministic heat equation. We also identify the fluctuations {\it of} the stationary solution itself and show that the rescaled averages $\int_{\mathbb R^d} {\mathrm d} x \varphi(x) \varepsilon^{1-\frac d2} [\mathfrak h^{\mathrm{st}}_{\varepsilon}(t,x)- \mathbb E(\mathfrak h^{\mathrm{st}}_{\varepsilon}(t,x))]$ converge to that of the {\emph stationary solution} of the stochastic heat equation with additive noise, but with (random) {\emph GFF marginals} (instead of flat initial condition).