We demonstrate a novel algorithm for generating stationary stochastic signals with a specified power spectral density (or equivalently, via the Wiener-Khinchin relation, a specified autocorrelation function) while satisfying constraints on the signal's probability density function. A tightly related problem has already been essentially solved by methods involving nonlinear filtering, however we use a fundamentally different approach involving optimization and stochastic interchange which immediately generalizes to generating signals with a broader range of statistics. This combination of optimization and stochastic interchange eliminates drawbacks associated with either method in isolation, improving the best-case scaling in runtime to generate a signal of length $n$ from $\mathcal{O}(n^2)$ for stochastic interchange on its own to $\mathcal{O}(n \: \text{log} \: n)$ without parallelization or $\mathcal{O}(n)$ with full parallelization. We demonstrate this speedup experimentally, and furthermore show that the signals we generate match the desired autocorrelation more accurately than those generated by stochastic interchange on its own. We observe that the signals we produce, unlike those generated by optimization on its own, are stationary.