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.

I think you can take Fourier tranform of the noise using fft, reject the parts (substiture zero for the power) you don't want and take inverse fft using ifft of the resulting power spectrum. That should theoretically give you the desired signal. And if you take the fft of the resulting signal you should see that power is zero for the rejected frequencies. — Some Guy

aksadv aksadv · Accepted Answer · 2017-01-26T01:19:42

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.

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

1 Answers