Band-limited random signal with arbitrary distribution?
Consider this approach:


General a white Gaussian random sequence.
Filter the white Gaussian sequence.  The output, which we'll call $z$, will be Gaussian because a linear combination of Gaussians is a Gaussian.
Pass the output $z$ from stage 2 through the transformation $g(z)=F_Y^{-1}(F_z(z))$, where $F_Z(z)$ is the Cumulative Distribution Function (CDF) of $z$ (the Gaussian CDF) and $F_Y^{-1}(y)$ is the inverse CDF of your target non-Gaussian distribution.  The output $y$ will now have the desired non-Gaussian distribution.


The only question is how to choose the filter taps in stage 2 so that the power spectral density at the output of the nonlinearity in stage 3 is what you want.  The sloppy but simple and often effective answer is that stage 3 typically doesn't change the power spectral density much, and you can often get away with just filtering $z$ so it has the desired power spectrum.  The more rigorous approach is to note that the autocorrelation of $y$ at lag $\tau$ $R_y(\tau)$ depends only on the autocorrelation of $z$ at the same lag $R_z(\tau)$.  Thus the problem of finding the needed $R_z(\tau)$ to achieve the target $R_y(\tau)$ reduces to separate optimization problems for each lag $\tau$.  It also turns out that $R_y(\tau)$ is a non-decreasing function of $R_z(\tau)$, so you can solve for $R_z(\tau)$ numerically fairly easily using very crude techniques.  Then, once you have the needed $R_z(\tau)$, you can design a filter to operate on the Gaussian sequence to achieve it.

Justification for the above statements can be found by searching for the following in Google Scholar: Cario &amp; Nelson, Modeling and Generating Random Vectors with Arbitrary Marginal Distributions and Correlation Matrix