The script generates spatial data with a scale-invariant power spectrum (1/f noise) and a normal error distribution.
The spectral density of the data is proportional to f^BETA, where f is the frequency and BETA is the spectral exponent (BETA=0 is white noise, BETA=-2 is Brown noise).
To generate Brown noise on a 10x10 grid you type
which gives a 10x10 matrix containihng the data.
1/f noise is not stationary, and so doesn't have a well defined variance. A variogram of the data show variance linearly increasing with increasing lag.
Time series data can also be generated. For example, a pink noise time series of lenth 1000 is produced by
Jon Yearsley (2021). Generate spatial data (https://www.mathworks.com/matlabcentral/fileexchange/5091-generate-spatial-data), MATLAB Central File Exchange. Retrieved .
Marcin is correct, the current version of the function (with BETA/2) gives a f^BETA power spectral density. To generate noise with a 1/f FFT *amplitude*, use BETA=-2.
However, I also agree with Marcin that the power spectral density looks very unusual, and is very noisy. For an alternative, see the generatepinknoise.m function from https://github.com/kendrickkay/knkutils
It has been pointed out to me that I have done my analysis wrong (and, obviously, so did @Christopher). I believe we have both plotted the signal amplitude (fft(x)) vs. frequency as opposed to the signal power spectral density (fft(x).^2) vs. frequency on a log-log plot (see e.g. https://uk.mathworks.com/help/signal/ug/power-spectral-density-estimates-using-fft.html). Anyway, the line 'S_f = (u.^2 + v.^2).^(BETA/2)' is in fact correct. I have also corrected my FileExchange submission to reflect that (https://uk.mathworks.com/matlabcentral/fileexchange/59305-randnd).
I have however found some other problems with this function. The way noise is generated (flat power spectral density rotated by a random phase into imaginary) causes the power spectral density to look very unnatural. Furthermore, from an unbeknownst to me reason, the noise generated by this function does not have a unit standard deviation. You can see that by comparing with 'randnd':
beta = 0; N = 2^12;
A1 = spatialPattern([1 N],beta);
A2 = randnd(beta,[1 N]);
A1fft = abs(fft(A1)).^2; A1fft = 2*A1fft(2:N/2);
A2fft = abs(fft(A2)).^2; A2fft = 2*A2fft(2:N/2);
k = 1:length(A1fft);
loglog(k,A1fft,k,A2fft); grid on; legend('spacialPattern','randnd');
xlim([min(k) max(k)]); xlabel('Wavenumber'); ylabel('Power/Wavenumber');
Before I've rated this function with 4 stars. Now:
+1 since the beta scaling is in fact correct
-1 for the unnatural power spectral density
-1 for the unnatural standard deviation
+1 for me feeling bad for badly assessing the 1/beta scaling beforehand
therefore still 4 stars.
Thanks, it's very useful.
I agree with @Christopher. It should be just BETA. Also; I wrote an updated version which works in an arbitrary number of dimensions (see Acknowledgments).
if I want generate a 2-D 1/f noise ,is it right ?spatialPattern([256,256],-1)?
I believe the line S_f = (u.^2 + v.^2).^(BETA/2); is not correct - since S_f is raised to the 0.5 power during the ifft, I believe this line should instead read S_f = (u.^2 + v.^2).^BETA. Double-checking the spectra by taking the fft of the resulting image seems to confirm that this gives the correct exponent.
IT is great. Thanks so much！
Great, thank you very much.
I use it for audio noise generation. It works great. It would be beneficial to eliminate non-audible components (below 20Hz) in order to improve the power density of the signal.
fine program - easy to use for pink noise for auditory purposes, haven't tried the spatial noises.
fantastic, works just like you want it to. It's a quick edit to make the default parameters a bit more user friendly, like so that Dim can be one number, and that automatically BETA sets to -2, for 1/f noise
Its of great use.
Inspired: Fractional Octave Band and A, B, C Weighting Filters DF2T SOS IIR Matlab and limited Labview, filt2 2D geospatial data filter, randnd, Noise-Power Spectrum, Nth_Oct_Hand_Arm_&_AC_Filter_Tool_Box
Find the treasures in MATLAB Central and discover how the community can help you!Start Hunting!