This work stems directly from the recent contribution by Carini, Furci, and Sommariva, who proposed a computationally efficient method for estimating the cross-power spectrum of a hidden multivariate stochastic process—such as brain activity—based on M/EEG observations. Their approach leverages sparse variational models and solves the resulting optimization problem via the Fast Iterative Shrinkage-Thresholding Algorithm (FISTA), achieving high specificity in functional connectivity estimation. Building upon their work, as a first step, we considered implementing both Armijo and adaptive backtracking strategies in the FISTA-based solver of the original method for estimating the Lipschitz constant dynamically during iterations. This significantly accelerates convergence compared to using conservative fixed step sizes, which are required in the original implementation for theoretical guarantees but are suboptimal in practice. Then, we explored the integration of Plug-and-Play (PnP) methods into this framework. Keeping the improved FISTA structure with backtracking, we consider a PnP denoiser in place of the $\ell^1$ proximal operator. PnP approaches introduce pre-trained neural network denoisers as implicit priors within proximal algorithms, effectively replacing explicit regularization terms such as the $\ell^1$ norm. This allows for model-aware regularization that adapts better to real data, while still retaining the interpretability and modular structure of classical optimization schemes. This hybrid formulation maintains the linear model fidelity while benefiting from learned priors embedded in the denoising network. Early comparisons suggest that this approach can enhance reconstruction quality and robustness without sacrificing computational tractability.