diff --git a/fsl/data/featanalysis.py b/fsl/data/featanalysis.py
index 448164cd467268d3e11eb18f51f8c51c68456a89..75fe43a61a9d2a58347455f2e50e46e3ea9e9b60 100644
--- a/fsl/data/featanalysis.py
+++ b/fsl/data/featanalysis.py
@@ -22,7 +22,6 @@ following functions are provided:
    loadDesign
    loadContrasts
    loadSettings
-   getEVNames
    getThresholds
    loadClusterResults
 
@@ -47,10 +46,12 @@ import                        logging
 import os.path             as op
 import numpy               as np
 
-import fsl.data.image      as fslimage
 import fsl.utils.path      as fslpath
 import fsl.utils.transform as transform
 
+from . import image as fslimage
+from . import          featdesign 
+
 
 log = logging.getLogger(__name__)
 
@@ -138,36 +139,6 @@ def getTopLevelAnalysisDir(path):
     return fslpath.shallowest(path, ['.ica', '.gica', '.feat', '.gfeat'])
 
 
-def loadDesign(featdir):
-    """Loads the design matrix from a FEAT directory.
-
-    Returns a ``numpy`` array containing the design matrix data, where the
-    first dimension corresponds to the data points, and the second to the EVs.
-
-    :arg featdir: A FEAT directory.
-    """
-
-    matrix    = None 
-    designmat = op.join(featdir, 'design.mat')
-
-    log.debug('Loading FEAT design matrix from {}'.format(designmat))
-
-    with open(designmat, 'rt') as f:
-
-        while True:
-            line = f.readline()
-            if line.strip() == '/Matrix':
-                break
-
-        matrix = np.loadtxt(f, ndmin=2)
-
-    if matrix is None or matrix.size == 0:
-        raise RuntimeError('{} does not appear to be a '
-                           'valid design.mat file'.format(designmat))
-
-    return matrix
-
-
 def loadContrasts(featdir):
     """Loads the contrasts from a FEAT directory. Returns a tuple containing:
     
@@ -265,6 +236,20 @@ def loadSettings(featdir):
     return settings
 
 
+def loadDesign(featdir, settings):
+    """Loads the design matrix from a FEAT directory.
+
+    :arg featdir:  A FEAT directory.
+
+    :arg settings: Dictionary containing FEAT settings (see
+                   :func:`loadSettings`).
+    
+    :returns:      a :class:`.FEATFSFDesign` instance which represents the
+                   design matrix.
+    """
+    return featdesign.FEATFSFDesign(featdir, settings)
+
+
 def getThresholds(settings):
     """Given a FEAT settings dictionary, returns a dictionary of
     ``{stat : threshold}`` mappings, containing the thresholds used
@@ -543,53 +528,3 @@ def getClusterMaskFile(featdir, contrast):
     """
     mfile = op.join(featdir, 'cluster_mask_zstat{}'.format(contrast + 1))
     return fslimage.addExt(mfile, mustExist=True)
-
-
-def getEVNames(settings):
-    """Returns the names of every EV in the FEAT analysis which has the given
-    ``settings`` (see the :func:`loadSettings` function).
-
-    An error of some sort will be raised if the EV names cannot be determined
-    from the FEAT settings.
-
-    :arg settings: A FEAT settings dictionary (see :func:`loadSettings`). 
-    """
-
-    numEVs    = int(settings['evs_real'])
-    titleKeys = [s for s in settings.keys() if s.startswith('evtitle')]
-    derivKeys = [s for s in settings.keys() if s.startswith('deriv_yn')]
-
-    def key(k):
-        return int(''.join([c for c in k if c.isdigit()]))
-        
-    titleKeys = sorted(titleKeys, key=key)
-    derivKeys = sorted(derivKeys, key=key)
-    evnames  = []
-
-    for titleKey, derivKey in zip(titleKeys, derivKeys):
-
-        # Figure out the ev number from
-        # the design.fsf key - skip over
-        # 'evtitle' (an offset of 7)
-        evnum = int(titleKey[7:])
-
-        # Sanity check - the evnum
-        # for the deriv_yn key matches
-        # that for the evtitle key
-        if evnum != int(derivKey[8:]):
-            raise RuntimeError('design.fsf seem to be corrupt')
-
-        title = settings[titleKey]
-        deriv = settings[derivKey]
-
-        if deriv == '0':
-            evnames.append(title)
-        else:
-            evnames.append(title)
-            evnames.append('{} - {}'.format(title, 'temporal derivative'))
-
-    if len(evnames) != numEVs:
-        raise RuntimeError('The number of EVs in design.fsf does not '
-                           'match the number of EVs in design.mat')
-
-    return evnames
diff --git a/fsl/data/featdesign.py b/fsl/data/featdesign.py
index d6ebf2a32dd116589e32855816cdb943ce5ddb37..f46316fbd269551f439c37c65e3acdacc64b268e 100644
--- a/fsl/data/featdesign.py
+++ b/fsl/data/featdesign.py
@@ -64,6 +64,18 @@ following types of *confound* EVs:
  - *Confound* EVs. These are any other EVs added by the user.
 
 
+Module contents
+---------------
+
+
+In addition to the :class:`FEATFSFDesign` class, this module contains a few
+other functions and classes that may be useful to advanced users.
+
+
+The :func:`loadDesignMat` function loads the ``design.mat`` file from a
+FEAT directory, and returns it as a numpy array.
+
+
 The following functions, defined in this module, will analyse a FEAT analysis
 to determine the contents of its design matrix (these functions are called by
 the :meth:`FEATFSFDesign.__init__` method, but may be called directly):
@@ -115,7 +127,7 @@ class FEATFSFDesign(object):
     with FSL 5.0.9 and older.
     """
     
-    def __init__(self, featDir, settings, designMatrix):
+    def __init__(self, featDir, settings):
         """Create a ``FEATFSFDesign``.
 
         :arg featDir:      Path to the FEAT directory.
@@ -123,14 +135,13 @@ class FEATFSFDesign(object):
         :arg settings:     A dictionary containing the FEAT analysis 
                            settings from its ``design.fsf`` file (see
                            :func:`.featanalysis.loadSettings`).
-                             
-        :arg designMatrix: The FEAT design matrix (a numpy array - see
-                           :func:`.featanalysis.loadDesign`). 
         """
 
-        # Get some information about the analysis
-        version = float(settings['version'])
-        level   = int(  settings['level'])
+        # Get the design matrix, and some
+        # information about the analysis
+        designMatrix = loadDesignMat(featDir)
+        version      = float(settings['version'])
+        level        = int(  settings['level'])
 
         # Print a warning if we're 
         # using an old version of FEAT
@@ -165,13 +176,31 @@ class FEATFSFDesign(object):
             # see the VoxelwisEV class.
             if ev.filename is not None: ev.image = fslimage.Image(ev.filename)
             else:                       ev.image = None
+
+
+    def getEVs(self):
+        """Returns a list containing the :class:`EV` instances that represent
+        each column of this ``FEATFSFDesign``.
+        """
+        return list(self.__evs)
  
     
-    def getDesign(self, x, y, z):
-        """Returns the design matrix for the specified voxel. """
+    def getDesign(self, voxel=None):
+        """Returns the design matrix for the specified voxel.
+
+        :arg voxel: A tuple containing the ``(x, y, z)`` voxel coordinates
+                    of interest. If ``None`` (the default), the design
+                    matrix is returned, with any voxelwise EV columns
+                    containing the mean voxelwise EV data.
+        """
 
         design = np.array(self.__design)
 
+        if voxel is None:
+            return design
+
+        x, y, z = voxel
+
         for ev in self.__evs:
 
             if not isinstance(ev, (VoxelwiseEV, VoxelwiseConfoundEV)):
@@ -581,3 +610,33 @@ def getHigherLevelEVs(featDir, settings, designMat):
         evs.append(VoxelwiseEV(len(evs), origIdx, title, filename))
         
     return evs
+
+
+def loadDesignMat(featdir):
+    """Loads the design matrix from a FEAT directory.
+
+    Returns a ``numpy`` array containing the design matrix data, where the
+    first dimension corresponds to the data points, and the second to the EVs.
+
+    :arg featdir: A FEAT directory.
+    """
+
+    matrix    = None 
+    designmat = op.join(featdir, 'design.mat')
+
+    log.debug('Loading FEAT design matrix from {}'.format(designmat))
+
+    with open(designmat, 'rt') as f:
+
+        while True:
+            line = f.readline()
+            if line.strip() == '/Matrix':
+                break
+
+        matrix = np.loadtxt(f, ndmin=2)
+
+    if matrix is None or matrix.size == 0:
+        raise FSFError('{} does not appear to be a '
+                       'valid design.mat file'.format(designmat))
+
+    return matrix
diff --git a/fsl/data/featimage.py b/fsl/data/featimage.py
index 1f4218abca4669aafb5704692847aea623c786c4..f2b60a3cf17ea782fa603abb64ea5678e4001e3e 100644
--- a/fsl/data/featimage.py
+++ b/fsl/data/featimage.py
@@ -29,8 +29,8 @@ class FEATImage(fslimage.Image):
         import fsl.data.featimage as featimage
 
         # You can pass in the name of the
-        # .feat/.gfeat directory, or any
-        # file contained within that directory.
+        # .feat directory, or any file
+        # contained within that directory.
         img = featimage.FEATImage('myanalysis.feat/filtered_func_data.nii.gz')
 
         # Query information about the FEAT analysis
@@ -74,10 +74,10 @@ class FEATImage(fslimage.Image):
         settings    = featanalysis.loadSettings( featDir)
 
         if featanalysis.hasStats(featDir):
-            design      = featanalysis.loadDesign(   featDir)
+            design      = featanalysis.loadDesign(   featDir, settings)
             names, cons = featanalysis.loadContrasts(featDir)
         else:
-            design      = np.zeros((0, 0))
+            design      = None
             names, cons = [], []
         
         fslimage.Image.__init__(self, path, **kwargs)
@@ -88,7 +88,6 @@ class FEATImage(fslimage.Image):
         self.__contrastNames = names
         self.__contrasts     = cons
         self.__settings      = settings
-        self.__evNames       = featanalysis.getEVNames(settings)
 
         self.__residuals     =  None
         self.__pes           = [None] * self.numEVs()
@@ -124,32 +123,47 @@ class FEATImage(fslimage.Image):
         """Returns ``True`` if the analysis for this ``FEATImage`` contains
         a statistical analysis.
         """
-        return self.__design.size > 0
+        return self.__design is not None
         
 
-    def getDesign(self):
+    def getDesign(self, voxel=None):
         """Returns the analysis design matrix as a :mod:`numpy` array
         with shape :math:`numPoints\\times numEVs`.
+        See :meth:`.FEATFSFDesign.getDesign`.
         """
-        return np.array(self.__design)
+        
+        if self.__design is None:
+            return None
+        
+        return self.__design.getDesign(voxel)
         
     
     def numPoints(self):
         """Returns the number of points (e.g. time points, number of
         subjects, etc) in the analysis.
         """
-        return self.__design.shape[0] 
+        if self.__design is None:
+            return None
+        
+        return self.__design.getDesign().shape[0]
 
     
     def numEVs(self):
         """Returns the number of explanatory variables (EVs) in the analysis.
         """
-        return self.__design.shape[1]
+        if self.__design is None:
+            return None 
+
+        return len(self.__design.getEVs())
 
 
     def evNames(self):
         """Returns a list containing the name of each EV in the analysis."""
-        return list(self.__evNames)
+        
+        if self.__design is None:
+            return None 
+        
+        return [ev.title for ev in self.__design.getEVs()]
 
     
     def numContrasts(self):
@@ -285,6 +299,9 @@ class FEATImage(fslimage.Image):
                         otherwise.
         """
 
+        if self.__design is None:
+            raise RuntimeError('No design')
+
         if not fullmodel:
             contrast = np.array(contrast)
             contrast = contrast / np.sqrt((contrast ** 2).sum())
@@ -295,7 +312,7 @@ class FEATImage(fslimage.Image):
         if len(contrast) != numEVs:
             raise ValueError('Contrast is wrong length')
         
-        X        = self.__design
+        X        = self.__design.getDesign(xyz)
         data     = self.data[x, y, z, :]
         modelfit = np.zeros(len(data))