how to generate an abitrary length colored noise when the PSD is know !

I want to generate colored noise samples of arbitrary length, I have the PSD of the noise in a vector of 512 elements !

[Patrik Ek]

It was a time since I came in contact with this, but if I recall correctly is Not the power spectral density the square magnitude of the transfer function? This means that to generate a signal from a psd, could you not take the square root (elementwise) of the psd and then generate a random phase. Then you know that H(w) = abs(H(w))*arg(H(w)), so you multiply the phase with the square root of the psd to get H(w). Then just take the ifft of H(w).

This will of course introduce some extra noise from ifft.

I am not sure if this is the best way to do things, but it should at least work I guess.

[abraham]

Thanks for the response, but the point is that doing as you explained , I can only generate a noise with the same length as the H(w) (the same length as the psd vector). I would like to have a large noise vector, greater than the length of the psd vector.

[Bjorn Gustavsson]

Presumably your signal is then bandlimited to frequencies corresponing to your 512-samples wide PSD, meaning that the PSD at higer frequencies is zero - so you should be able to simply pad the PSD with zeros up to the desired length.

Answer 1

[Image Analyst]

You need to use inverse transform sampling. By doing that you can get any distribution you want from a uniform distribution. See http://en.wikipedia.org/wiki/Inverse_transform_sampling

I've also attached a demo where I pull random values from a Rayleigh distribution. Feel free to modify it.

You might also be interested in this: http://www.mathworks.com/matlabcentral/fileexchange/7309-randraw

[Image Analyst]

abraham's "Answer" moved here:

Could you explain the concept of inverse transform sampling ?

[Image Analyst]

Think of a look up table or graph of "out vs. in" If it was a straight 45 degree line, then the output equals the input. If the slope was 2 then the output equals twice the input. What if it's some weird-shaped curve? Then you send in an input and read off the output from the curve. What if the curve was the cumulative distribution function of some probability distribution? Like, say, a Gaussian? The curve is flat to start out, then steep, then flat again. OK, but let's say we invert/flip that so that the curve is steep for low "in" values, then flat for middle "in" values, then steep again for high "in" values. OK, much better. Now, what happens if you we send in a uniform distribution and use that as a look up table? Let's say you send in a thousand low (bottom third) values, a thousand middle (middle third) values, and a thousand high (upper third) values, just like you'd get if you pulled values from a uniform distribution. Well the thousand low values get sent all over the lower, say 40%, and the upper thousand input values get spread over the upper, say, 40%, but the middle 1000 values get all crammed in the 40-60% range because the inverse CDF curve is flatter there in the middle. So that's what we want for a Gaussian - more samples per unit distance in the middle, and less at the tails. Does that help explain how it works? It helps to draw the curve and axes and see where input values get mapped to.

[Image Analyst]

Abraham: Not sure why you never responded. I just made up another demo for someone else that does an inverse power law, and thought you'd also be interested. See attached demo.

With this demo, you can adapt it to any PDF you want since it calculates the CDF numerically and you don't need to be able to integrate the PDF formula analytically with calculus. So you should be able to modify it for your PDF. Please mark the answer as accepted if it helps you.

[Bjorn Gustavsson]

That would lead to a random-number sequence with a given probability distribution function, but the OP asked for a given power spectral density, as far as I understood the question.

Answer 2

[Youssef Khmou]

You proceed using Inverse Fourier transform , here is how you can start , let us consider zero mean random variable of length 200 :

 x=randn(200,1);
 f=fft(x); %  Discret Fourier Transform with same length 
 p=f.*conj(f); % Power
 figure; plot(p)
 C=real(ifft(f)); %  inverse transform of the power 
 figure; plot(C,x) %

the original signal or the ones resulting from inverse transformation are similar, so suppose you have vector of Power spectrum, you ifft with the desired length N :

 X=real(ifft(x,N));