How do I create band-limited (100-640 Hz) white Gaussian noise?

Question

I would like to create 500 ms of band-limited (100-640 Hz) white Gaussian noise with a (relatively) flat frequency spectrum. The noise should be normally distributed with mean = ~0 and 99.7% of values between ± 2 (i.e. standard deviation = 2/3). My sample rate is 1280 Hz; thus, a new amplitude is generated for each frame.

duration  = 500e-3;
rate      = 1280;
amplitude = 2;

npoints   = duration * rate;
noise     = (amplitude/3)* randn( 1, npoints ); 
% Gaus distributed white noise; mean = ~0; 99.7% of amplitudes between ± 2.
time      = (0:npoints-1) / rate

Could somebody please show me how to filter the signal for the desired result (i.e. 100-640 Hz)? In addition, I was hoping somebody could also show me how to generate a graph to illustrate that the frequency spectrum is indeed flat.

I intend on importing the waveform to Signal (CED) to output as a form of transcranial electrical stimulation.

Answer 1

The following is Matlab implementation of the method alluded to by "Some Guy" in a comment to your question.

% In frequency domain, white noise has constant amplitude but uniformly
% distributed random phase. We generate this here. Only half of the
% samples are generated here, the rest are computed later using the complex
% conjugate symmetry property of the FFT (of real signals).
X         = [1; exp(i*2*pi*rand(npoints/2-1,1)); 1]; % X(1) and X(NFFT/2) must be real 

% Identify the locations of frequency bins. These will be used to zero out
% the elements of X that are not in the desired band
freqbins  = (0:npoints/2)'/npoints*rate;

% Zero out the frequency components outside the desired band 
X(find((freqbins < 100) | (freqbins > 640))) = 0;

% Use the complex conjugate symmetry property of the FFT (for real signals) to
% generate the other half of the frequency-domain signal
X         = [X; conj(flipud(X(2:end-1)))];

% IFFT to convert to time-domain
noise     = real(ifft(X));

% Normalize such that 99.7% of the times signal lies between ±2
noise     = 2*noise/prctile(noise, 99.7);

Statistical analysis of around a million samples generated using this method results in the following spectrum and distribution:

Firstly, the spectrum (using Welch method) is, as expected, flat in the band of interest:

Also, the distribution, estimated using histogram of the signal, matches the Gaussian PDF quite well.