Plug-and-play (PnP) approaches currently achieve state-of-the-art quality in image restoration. These methods rely on a Gaussian denoiser, often parametrized by an artificial neural network (ANN) that learns the image key features. This work adapts the PnP approach to bivariate time series, with an emphasis on preserving the polarization, that is, the geometrical dependence between the two components of the signal. It designs an ANN denoiser in the time-frequency domain, using exclusively a synthetic dataset and data-augmentation operations. The interest of the approach is demonstrated with a non-trivial application to gravitational wave astronomy. Up to our knowledge, this work is one of the first applications of a PnP approach to inverse problems involving (multivariate) time series.