Band-limited random signal with arbitrary distribution?

Question

I'd like to generate a random discrete-time signal that is band-limited to some bandwidth B (by means of a digital filter, ie in MATLAB). The catch is that I'd like this signal to have an arbitrary PDF, for example, a uniform distribution.

If I start with a uniform distribution, and then filter to band-limit the signal to B, I end up with something that's basically gaussian. Since my band-limiting FIR filter has a large number of taps (say 128), this makes sense given the Central Limit Theorem, as I am basically summing a number of IID random variables.

What I don't understand, is if it is possible to now transform this signal to a new distribution, such as uniform, while maintaining the band-limitedness. I understand that I can transform a RV from one PDF to another by basically integrating/mapping over the CDF (hand-wavy, it's been a while!), but I believe this significantly changes the frequency content of the band-limited signal. For example, if I have a gaussian distribution, and I want to transform to a uniform, I basically map values "near the mean" of the gaussian to something "closer to the max" in the uniform RV. But now the frequency content has changed, and my spectrum will be smeared.

Does what I'm describing make sense, and has anyone looked at something similar? Perhaps what I need to do is look at transforming the signal as a vector, not a scalar (ie treat it as a random process, not individual random variables). But this is where I get a little out of my depth.

Answer 1

Consider this approach:

1. General a white Gaussian random sequence.
2. Filter the white Gaussian sequence. The output, which we'll call $z$, will be Gaussian because a linear combination of Gaussians is a Gaussian.
3. 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 & Nelson, Modeling and Generating Random Vectors with Arbitrary Marginal Distributions and Correlation Matrix

Answer 2

I always found the meaning of the probability density function (PDF) of a discrete-time signal a bit puzzling. If the discrete-time signal is converted to a continuous-time signal by band-limited interpolation or by a digital to analog converter (DAC), then the resulting PDF will in general not be exactly the same as before.

Perhaps you could just resample to (or generate the signal at) a sampling frequency that is twice the band limit. Then you can waveshape the samples to a PDF of your choice, and you still retain the bandlimit of half the sampling frequency, by definition.