diff --git a/README.rst b/README.rst
index 2baae82..b38389c 100644
--- a/README.rst
+++ b/README.rst
@@ -51,3 +51,24 @@ Example
 	#plt.grid(True)
 	#plt.show()
 	
+	
+.. code:: python
+    # generate several time series of independent indentically distributed variables 
+    # repeat the simulation of each variable multiple times
+	import colorednoise as cn
+    n_repeats   = 10   # repeat simulatons
+    n_variables = 5    # independent variables in each simulation
+    timesteps   = 1000 # number of timesteps for each variable
+    y = cn.powerlaw_psd_gaussian(1, (n_repeats, n_variables, timesteps))
+    
+    # the expected variance of for each variable is 1, but each realisation is different
+    print(y.std(axis=-1))
+    
+.. code:: python
+    # generate a broken power law spectrum: white below a frequency of 
+	import colorednoise as cn
+    y = cn.powerlaw_psd_gaussian(1, 10**5, fmin=.05)
+    s, f = mlab.psd(y, NFFT=2**9)
+    #plt.loglog(f,s)
+    #plt.grid(True)
+    #plt.show()
diff --git a/colorednoise.py b/colorednoise.py
index 051818e..fae2392 100644
--- a/colorednoise.py
+++ b/colorednoise.py
@@ -37,8 +37,15 @@ def powerlaw_psd_gaussian(exponent, size, fmin=0):
 
     fmin : float, optional
         Low-frequency cutoff.
-        Default: 0 corresponds to original paper. It is not actually
-        zero, but 1/samples.
+        Default: 0 corresponds to original paper. 
+        
+        The power-spectrum below fmin is flat. fmin is defined relative
+        to a unit sampling rate (see numpy's rfftfreq). For convenience,
+        the passed value is mapped to max(fmin, 1/samples) internally
+        since 1/samples is the lowest possible finite frequency in the
+        sample. The largest possible value is fmin = 0.5, the Nyquist
+        frequency. The output for this value is white noise.
+
 
     Returns
     -------
@@ -67,9 +74,14 @@ def powerlaw_psd_gaussian(exponent, size, fmin=0):
     # Use fft functions for real output (-> hermitian spectrum)
     f = rfftfreq(samples)
     
+    # Validate / normalise fmin
+    if 0 <= fmin <= 0.5:
+        fmin = max(fmin, 1./samples) # Low frequency cutoff
+    else:
+        raise ValueError("fmin must be chosen between 0 and 0.5.")
+    
     # Build scaling factors for all frequencies
-    s_scale = f
-    fmin = max(fmin, 1./samples) # Low frequency cutoff
+    s_scale = f    
     ix   = npsum(s_scale < fmin)   # Index of the cutoff
     if ix and ix < len(s_scale):
         s_scale[:ix] = s_scale[ix]
diff --git a/setup.py b/setup.py
index 1fd00ed..3aaffde 100644
--- a/setup.py
+++ b/setup.py
@@ -2,7 +2,7 @@
 
 setup(
     name='colorednoise',
-    version='1.1.1',
+    version='1.2.0',
     description='Generate Gaussian (1/f)**beta noise (e.g. pink noise)',
     long_description="""
         Generate Gaussian distributed noise with a power law spectrum.