Colouring Noise – Generating Noise with a Specific Power Spectral Density(blog.ioces.com)
blog.ioces.com
Colouring Noise – Generating Noise with a Specific Power Spectral Density
https://blog.ioces.com/matt/posts/colouring-noise/
9 comments
We use a similar method to simulate 2D atmospheric turbulence for adaptive optics simulations in astronomy. The structure of turbulence is described by an analytical power spectral density function (classically Kolmogorov spectrum). We generate random coefficients in the Fourier space and then weight them by the PSD. After an inverse Fourier transform we obtain a beautiful draw of turbulence with the desired spatial structure.
Time domain methods to make pink noise:
https://www.firstpr.com.au/dsp/pink-noise/ http://stenzel.waldorfmusic.de/post/a-new-shade-of-pink/
https://www.firstpr.com.au/dsp/pink-noise/ http://stenzel.waldorfmusic.de/post/a-new-shade-of-pink/
The "integration to zero" effect the author observes is a consequence of the fact that the Fourier Transform is meant only for periods signals. Practically, we all often toss any time series data into an FFT to get a look at it in frequency space, but the math is only exactly correct with the time series samples used from one cycle of a period signal that repeats forever.
By generating the frequency domain Fourier coefficients and then transforming them to the time domain, the author is guaranteed to get a signal that starts and ends at the same value, because the transformation assumes it's a periodic signal and so the first and last values of the period must be the same.
By generating the frequency domain Fourier coefficients and then transforming them to the time domain, the author is guaranteed to get a signal that starts and ends at the same value, because the transformation assumes it's a periodic signal and so the first and last values of the period must be the same.
The author identified the cause of the integration-to-zero behavior - they've set the DC (0Hz) component to 0, so the signal must sum to 0.
Any length-N signal can be thought of as one period of an N-periodic signal, so the FFT will be exact. It's not the first and last samples that need to be identical, it's the 1st sample needs to match the (N+1)th sample (1-indexed), which comes right after the last sample of the original finite signal.
For example if you have the length-8 signal 12345678, the FFT is valid - you can think of it as a slice of the signal ...123456781234567812345678...
Any length-N signal can be thought of as one period of an N-periodic signal, so the FFT will be exact. It's not the first and last samples that need to be identical, it's the 1st sample needs to match the (N+1)th sample (1-indexed), which comes right after the last sample of the original finite signal.
For example if you have the length-8 signal 12345678, the FFT is valid - you can think of it as a slice of the signal ...123456781234567812345678...
Also, it seems like the author might not be aware that the FT of a Gaussian is another Gaussian. I implemented this method (Timmer and Konig) a few months ago, and I think there are some other statistical tricks that work out nice for Gaussians.
In cosmology we use a similar procedure to generate initial conditions of sample universes. Just then the noise is in 3D. The article raises a question about the Gaussian nature of the generated noise. To generate non-gaussian noise, you need a way to correlate the phases, otherwise the central-limit theorem will always kick in.
The 0 freq component is just the general, time independent bias of the signal. If it is 0 then of course the integral of something which is neither positively nor negatively biased overall will average out to null.
The apparent convergence to zero happens because the amount of sample points is exactly the amount of data that can be represented with the number of frequency components you are calculating in your FT. Obviously net error in electronic readings don't really go to zero when you have exactly 1 million of them. Use fewer frequencies in your fft, see the convergence sooner. Use more, see it later.
The libraries aren't that impressive. The FFT is simple enough to write in any coding language with basic trig and arithmetic. It could certainly fit into a similar amount of code to what the author has written.
The apparent convergence to zero happens because the amount of sample points is exactly the amount of data that can be represented with the number of frequency components you are calculating in your FT. Obviously net error in electronic readings don't really go to zero when you have exactly 1 million of them. Use fewer frequencies in your fft, see the convergence sooner. Use more, see it later.
The libraries aren't that impressive. The FFT is simple enough to write in any coding language with basic trig and arithmetic. It could certainly fit into a similar amount of code to what the author has written.
If you like this, you may find the allantools library useful.
https://github.com/aewallin/allantools
https://github.com/aewallin/allantools
That's exactly the library I'm using for my Allan deviation plots. The documentation and interface are written more in relation to testing frequency standards, but works equally as well for gyro data.