diff --git a/Help/Makefile b/Help/Makefile
--- a/Help/Makefile
+++ b/Help/Makefile
@@ -1,7 +1,9 @@
 all:
 	ghc --make -Wall -O2 browse
-	ghc --make -Wall -O2 parse
+	ghc --make -Wall -O2 parser
 	ghc --make -Wall -O2 play
+	ghc --make -Wall -O2 plot
+	ghc --make -Wall -O2 stat
 
 clean:
-	rm *.o *.hi browse parse play
+	rm -f *.o *.hi browse parser play plot stat
diff --git a/Help/browse.hs b/Help/browse.hs
--- a/Help/browse.hs
+++ b/Help/browse.hs
@@ -1,5 +1,4 @@
 import Control.Monad
-import qualified Data.IntMap as Map
 import qualified Sound.Analysis.Meapsoft as M
 import System.Environment
 import System.IO
@@ -12,22 +11,25 @@
 readDef :: (Read a) => a -> String -> a
 readDef x s = maybe x id (readMaybe s)
 
-browse :: Map.IntMap M.Frame -> IO ()
-browse f = do putStr "feature> "
-              n <- fmap (readDef 0) getLine
-              putStr "chunk> "
-              m <- fmap (readDef 0) getLine
-              putStrLn (show ( "feature", n, M.features !! n
-                             , "chunk", m ))
-              let (_, _, d) = M.raw_data (f Map.! m)
-              putStrLn (show (d !! n))
+browse :: M.MEAP -> IO ()
+browse m = do
+  putStr "feature> "
+  i <- fmap (readDef 0) getLine
+  putStr "frame> "
+  j <- fmap (readDef 0) getLine
+  let f = M.features m !! i
+      l = M.feature_column f
+      r = l + M.feature_degree f - 1
+  putStrLn (show ( "feature", f))
+  putStrLn (show ("chunk", map (\c -> M.position m (j, c)) [l..r]))
 
 main :: IO ()
-main = do hSetBuffering stdout NoBuffering
-          a <- getArgs
-          unless (length a == 1) (error "feature-file")
-          let [f] = a
-          fs <- M.read_file f
-          putStrLn (show ("nframes", length fs, "nfeatures", length M.features))
-          let m = Map.fromList (zip [0..] fs)
-          forever (browse m)
+main = do
+  hSetBuffering stdout NoBuffering
+  a <- getArgs
+  unless (length a == 1) (error "feature-file")
+  let [f] = a
+  (Right m) <- M.read_meap f
+  putStrLn (show ("nframes", M.n_frames m,
+                  "nfeatures", length (M.features m)))
+  forever (browse m)
diff --git a/Help/parse.hs b/Help/parse.hs
deleted file mode 100644
--- a/Help/parse.hs
+++ /dev/null
@@ -1,23 +0,0 @@
-import Control.Monad
-import Data.List
-import qualified Sound.Analysis.Meapsoft as M
-import System.Environment
-import System.FilePath
-import System.IO
-
-replace :: Eq a => a -> a -> [a] -> [a]
-replace _ _ [] = []
-replace a b (x:xs) | x == a = b : replace a b xs
-                   | otherwise = x : replace a b xs
-
-main :: IO ()
-main = do
-  a <- getArgs
-  unless (length a == 1) (error "analysis-file")
-  let [fn] = a
-      vn = replace '-' '_' (dropExtensions (takeFileName fn))
-  fs <- M.read_file fn
-  putStrLn ("-- " ++ fn)
-  putStrLn "import Sound.Analysis.Meapsoft"
-  putStrLn (vn ++ " :: [Frame]")
-  putStrLn (vn ++ " = " ++ show fs)
diff --git a/Help/parser.hs b/Help/parser.hs
new file mode 100644
--- /dev/null
+++ b/Help/parser.hs
@@ -0,0 +1,12 @@
+import Control.Monad
+import qualified Sound.Analysis.Meapsoft as M
+import System.Environment
+
+main :: IO ()
+main = do
+  a <- getArgs
+  unless (length a == 1)
+         (error "analysis-file")
+  let [fn] = a
+  (Right h) <- M.read_header fn
+  putStrLn (show h)
diff --git a/Help/play.hs b/Help/play.hs
--- a/Help/play.hs
+++ b/Help/play.hs
@@ -1,4 +1,5 @@
 import Control.Monad
+import Data.Function
 import Data.List
 import qualified Sound.Analysis.Meapsoft as M
 import Sound.OpenSoundControl
@@ -14,11 +15,13 @@
           q = bufSampleRate KR b
           z = envTrapezoid 0.95 0.5 d g
           e = envGen KR 1 1 0 1 RemoveSynth z
-          p = playBuf 1 b r 1 (s * q) NoLoop
+          p = playBuf 1 b r 1 (s * q) NoLoop DoNothing
 
 mk_grain :: (Transport t) => t -> Double -> Double -> IO ()
-mk_grain fd s d = send fd (s_new "grain" (-1) AddToTail 1 [ ("start", s)
-                                                          , ("dur", d) ])
+mk_grain fd s d =
+    let xs = [ ("start", s)
+             , ("dur", d) ]
+    in send fd (s_new "grain" (-1) AddToTail 1 xs)
 
 tone :: UGen
 tone = p * e
@@ -30,9 +33,12 @@
           p = saw AR (MCE [f, f*2, f*3]) * (MCE [0.5, 0.2, 0.1])
 
 mk_tone :: (Transport t) => t -> Double -> Double -> IO ()
-mk_tone fd f d = send fd (s_new "tone" (-1) AddToTail 1 [ ("freq", f)
-                                                        , ("dur", d) ])
+mk_tone fd f d =
+    let xs = [ ("freq", f)
+             , ("dur", d) ]
+    in send fd (s_new "tone" (-1) AddToTail 1 xs)
 
+-- (Temporal_Scalar, [Segments], [Freq])
 type Spec = (Double, [(Double, Double)], [Double])
 
 plyr :: Transport t => t -> Spec -> IO ()
@@ -43,12 +49,7 @@
                                  plyr fd (t, xs, fs)
 plyr _ _ = undefined
 
-interleave :: [a] -> [a] -> [a]
-interleave [] b = b
-interleave a [] = a
-interleave (a:as) (b:bs) = a : b : interleave as bs
-
-run_grn :: FilePath -> FilePath -> ([M.Frame] -> Spec) -> IO ()
+run_grn :: FilePath -> FilePath -> (M.MEAP -> Spec) -> IO ()
 run_grn sf ff rule = withSC3 g
     where g fd = do reset fd
                     let g0 = out 0 (grain 10)
@@ -56,94 +57,145 @@
                     async fd (b_allocRead 10 sf 0 0)
                     async fd (d_recv (synthdef "grain" g0))
                     async fd (d_recv (synthdef "tone" g1))
-                    c <- M.read_file ff
-                    plyr fd (rule c)
+                    (Right m) <- M.read_meap ff
+                    plyr fd (rule m)
 
-freq_1 :: [M.Frame] -> [Double]
-freq_1 = map M.avg_freq_simple
+col :: String -> M.MEAP -> [Double]
+col n m =
+    let f = M.required_feature n (M.features m)
+        j = M.feature_column f
+    in M.column_l m j
 
-freq_2 :: [M.Frame] -> [Double]
-freq_2 = map M.avg_spec_centroid
+freq_1 :: M.MEAP -> [Double]
+freq_1 = col "AvgFreqSimple"
 
-rule_1 :: [M.Frame] -> Spec
-rule_1 c = (0.75, M.segments c, freq_1 c)
+-- forwards
+rule_1 :: M.MEAP -> Spec
+rule_1 c = (0.75, M.segments_l c, freq_1 c)
 
-rule_2 :: [M.Frame] -> Spec
-rule_2 c = (0.75, reverse (M.segments c), freq_2 c)
+freq_2 :: M.MEAP -> [Double]
+freq_2 = col "AvgSpecCentroid"
 
-rule_3 :: [M.Frame] -> Spec
-rule_3 c = (0.75, sortBy f (M.segments c), freq_1 c)
-    where f (_, d1) (_, d2) = compare d1 d2
+-- backwards
+rule_2 :: M.MEAP -> Spec
+rule_2 c = (0.5, reverse (M.segments_l c), reverse (freq_2 c))
 
-rule_4 :: [M.Frame] -> Spec
-rule_4 c = (0.75, interleave s (reverse s), freq_1 c)
-    where f (_, d1) (_, d2) = compare d1 d2
-          s = sortBy f (M.segments c)
+gen_order_1 :: M.MEAP -> String -> [Int]
+gen_order_1 m f =
+    let cs = zip [0..] (col f m)
+        cs' = sortBy (compare `on` snd) cs
+    in map fst cs'
 
-with_times :: (M.Frame -> a) -> M.Frame -> (Double, Double, a)
-with_times f c = (M.onset_time c, M.duration c, f c)
+apply_order :: [Int] -> [a] -> [a]
+apply_order o xs = map (\i -> xs !! i) o
 
-rule_5 :: [M.Frame] -> Spec
-rule_5 cs = (0.75, map g x1, map h x1)
-    where f (_, _, c1) (_, _, c2) = compare c1 c2
-          g (t, d, _) = (t, d)
-          h (_, _, fr) = fr
-          x0 = map (with_times M.avg_spec_centroid) cs
-          x1 = sortBy f x0
+srt_rule :: String -> (M.MEAP -> [Double]) -> M.MEAP -> Spec
+srt_rule n r m =
+    let o = gen_order_1 m n
+        s = apply_order o (M.segments_l m)
+        f = apply_order o (r m)
+    in (0.5, s, f)
 
-reverse_spec :: Spec -> Spec
-reverse_spec (t, x, f) = (t, reverse x, reverse f)
+-- sorted by segment length (ascending)
+rule_3 :: M.MEAP -> Spec
+rule_3 = srt_rule "chunk_length" freq_1
 
-rule_6 :: [M.Frame] -> Spec
-rule_6 = reverse_spec . rule_5
+-- sorted by frequency (ascending)
+rule_4 :: M.MEAP -> Spec
+rule_4 = srt_rule "AvgFreqSimple" freq_1
 
-rule_7 :: [M.Frame] -> Spec
-rule_7 cs = (0.75, map g x1, map h x1)
-    where f (_, _, c1) (_, _, c2) = compare c1 c2
-          g (t, d, _) = (t, d)
-          h (_, _, fr) = fr
-          x0 = map (with_times M.avg_freq_simple) cs
-          x1 = sortBy f x0
+-- sorted by spectral centroid (ascending)
+rule_5 :: M.MEAP -> Spec
+rule_5 = srt_rule "AvgSpecCentroid" freq_2
 
-rule_8 :: [M.Frame] -> Spec
-rule_8 = reverse_spec . rule_7
+-- sorted by spectral stability (ascending)
+rule_6 :: M.MEAP -> Spec
+rule_6 = srt_rule "SpectralStability" freq_2
 
+gen_sel_1 :: M.MEAP -> String -> (Double -> Bool) -> [Int]
+gen_sel_1 m f p =
+    let cs = zip [0..] (col f m)
+        cs' = filter (p . snd) cs
+    in map fst cs'
+
+sel_rule :: String -> (Double -> Bool) -> (M.MEAP -> [Double]) -> M.MEAP -> Spec
+sel_rule n p r m =
+    let o = gen_sel_1 m n p
+        s = apply_order o (M.segments_l m)
+        f = apply_order o (r m)
+    in (0.5, s, f)
+
+-- filtered by frequency (> 600)
+rule_7 :: M.MEAP -> Spec
+rule_7 = sel_rule "AvgFreqSimple" (> 600) freq_1
+
+-- filtered by spectral centroid (> 1200)
+rule_8 :: M.MEAP -> Spec
+rule_8 = sel_rule "AvgSpecCentroid" (> 1200) freq_2
+
+interleave :: [a] -> [a] -> [a]
+interleave [] b = b
+interleave a [] = a
+interleave (a:as) (b:bs) = a : b : interleave as bs
+
+-- interleaving of rules 4 & 5
+rule_9 :: M.MEAP -> Spec
+rule_9 m =
+    let (_, s4, f4) = rule_4 m
+        (_, s5, f5) = rule_5 m
+    in (0.75, interleave s4 s5, interleave f4 f5)
+
+reverse_spec :: Spec -> Spec
+reverse_spec (t, x, f) = (t, reverse x, reverse f)
+
+-- reverse of rule 9
+rule_10 :: M.MEAP -> Spec
+rule_10 = reverse_spec . rule_9
+
 scale_start_only :: Double -> Spec -> Spec
 scale_start_only n (t, x, f) = (t, map g x, f)
     where g (s, d) = (s * n, d)
 
-rule_9 :: [M.Frame] -> Spec
-rule_9 = (scale_start_only 0.75) . rule_7
+-- rule 10 with start times compressed but durations kept equal
+rule_11 :: M.MEAP -> Spec
+rule_11 = (scale_start_only 0.5) . rule_10
 
 scale_dur_only :: Double -> Spec -> Spec
 scale_dur_only n (t, x, f) = (t, map g x, f)
     where g (s, d) = (s, d * n)
 
-rule_10 :: [M.Frame] -> Spec
-rule_10 = (scale_dur_only 0.65) . rule_7
+-- rule 7 with durations compressed but start times retained
+rule_12 :: M.MEAP -> Spec
+rule_12 = (scale_dur_only 0.5) . rule_7
 
 fix_dur :: Double -> Spec -> Spec
 fix_dur n (t, x, f) = (t, map g x, f)
     where g (s, _) = (s, n)
 
-rule_11 :: [M.Frame] -> Spec
-rule_11 = (fix_dur 0.075) . rule_7
+-- rule 7 with durations fixed at 0.25 seconds
+rule_13 :: M.MEAP -> Spec
+rule_13 = (fix_dur 0.25) . rule_7
 
+rules :: [M.MEAP -> Spec]
+rules = [ rule_1
+        , rule_2
+        , rule_3
+        , rule_4
+        , rule_5
+        , rule_6
+        , rule_7
+        , rule_8
+        , rule_9
+        , rule_10
+        , rule_11
+        , rule_12
+        , rule_13 ]
+
 main :: IO ()
 main = do
   a <- getArgs
   unless (length a == 3) (error "audio-file feature-file rule-number")
   let [w_fn, f_fn, r] = a
       rn = read r - 1
-  unless (rn >= 0 && rn < 11) (error "unknown rule")
-  run_grn w_fn f_fn ( [ rule_1
-                      , rule_2
-                      , rule_3
-                      , rule_4
-                      , rule_5
-                      , rule_6
-                      , rule_7
-                      , rule_8
-                      , rule_9
-                      , rule_10
-                      , rule_11 ] !! rn)
+  unless (rn >= 0 && rn < length rules) (error "unknown rule")
+  run_grn w_fn f_fn ( rules !! rn )
diff --git a/Help/plot.hs b/Help/plot.hs
new file mode 100644
--- /dev/null
+++ b/Help/plot.hs
@@ -0,0 +1,41 @@
+-- Use hmatrix and gnuplot bindings to plot unit dimensioned features.
+
+import Control.Monad
+import Data.List
+import qualified Data.Packed.Convert as G
+import qualified Data.Packed.Matrix as G
+import qualified Data.Packed.Vector as G
+import qualified Graphics.Gnuplot.Simple as P
+import qualified Numeric.GSL.Vector as G
+import qualified Sound.Analysis.Meapsoft as M
+import System.Environment
+import System.IO
+
+normalize :: G.Vector Double -> G.Vector Double
+normalize v =
+    let l = G.vectorMin v
+        r = G.vectorMax v
+        d = r - l
+        v' = G.vectorMapValR G.AddConstant (negate l) v
+    in G.vectorMapValR G.Scale (recip d) v'
+
+main :: IO ()
+main =
+    do a <- getArgs
+       unless (length a > 1)
+              (error ("feature-file feature ..\n" ++ 
+                      show M.feature_names))
+       let (fn:ns) = a
+       (Right mp) <- M.read_meap fn
+       let fs = M.features mp
+           uf = map (\n -> M.required_feature n fs) ("onset_time":ns)
+           ar = M.uarray_data mp
+           m = G.matrixFromArray ar
+           ufc = map M.feature_column uf
+           vs = G.toRows (G.extractRows ufc (G.trans m))
+           vs' = map normalize vs
+           (ot:ls) = map G.toList vs'
+           ls' = map (zip ot) ls
+       P.plotPaths 
+        [P.XLabel "onset_time", P.YLabel (intercalate "," ns)]
+        ls'
diff --git a/Help/stat.hs b/Help/stat.hs
new file mode 100644
--- /dev/null
+++ b/Help/stat.hs
@@ -0,0 +1,35 @@
+-- Use hmatrix GSL bindings to do basic statistical analysis.
+
+import Control.Monad
+import Data.Function
+import Data.List
+import qualified Data.Packed.Convert as G
+import qualified Data.Packed.Matrix as G
+import qualified Data.Packed.Vector as G
+import qualified Numeric.GSL.Vector as G
+import qualified Sound.Analysis.Meapsoft as M
+import System.Environment
+import System.IO
+
+main :: IO ()
+main =
+    do a <- getArgs
+       unless (length a == 1)
+              (error "feature-file")
+       let [fn] = a
+       (Right mp) <- M.read_meap fn
+       let nf = fromIntegral (M.n_frames mp)
+           uf = filter ((== 1) . M.feature_degree) (M.features mp)
+           ar = M.uarray_data mp
+           m = G.matrixFromArray ar
+           ufc = map M.feature_column uf
+           vs = G.toRows (G.extractRows ufc (G.trans m))
+           r = zip5 (map M.feature_name uf)
+                    (zip (map G.vectorMin vs)
+                         (map G.vectorMinIndex vs))
+                    (zip (map G.vectorMax vs)
+                         (map G.vectorMaxIndex vs))
+                    (zip (map (G.toScalarR G.AbsSum) vs)
+                         (map (G.toScalarR G.Norm2) vs))
+                    (map ((/ nf) . G.toScalarR G.AbsSum) vs)
+       mapM_ (putStrLn . show) r
diff --git a/README b/README
--- a/README
+++ b/README
@@ -1,10 +1,11 @@
 hmeap - Haskell Parser for Meapsoft Analysis Files
 
-At present only supports analysis files with all
-features extracted.
+Supports analysis files with any number of features
+extracted.  The help files illustrate use with
+hmatrix, gnuplot and hsc3.
 
   meapsoft: http://labrosa.ee.columbia.edu/meapsoft/
   hmeap: http://slavepianos.org/rd/f/896678/
 
-(c) rohan drape, 2007-2008
-    GPL.2, http://gnu.org/copyleft/
+(c) rohan drape, 2007-2009
+    gpl.2, http://gnu.org/copyleft/
diff --git a/Sound/Analysis/Meapsoft.hs b/Sound/Analysis/Meapsoft.hs
--- a/Sound/Analysis/Meapsoft.hs
+++ b/Sound/Analysis/Meapsoft.hs
@@ -1,108 +1,96 @@
--- | A data representation of, and parser for, MEAPsoft analysis frames.
-module Sound.Analysis.Meapsoft ( Frame(..), Segment
-                               , features
-                               , read_file
-                               , segments) where
+-- | Data representation and parser for MEAPsoft analysis frames.
+module Sound.Analysis.Meapsoft
+    ( MEAP
+    , read_meap
+    , features
+    , n_frames
+    , uarray_data
+    , n_columns
+    , frame_l
+    , column_l
+    , position
+    , segments_l
+    , feature_names
+    , module Sound.Analysis.Meapsoft.Data
+    , module Sound.Analysis.Meapsoft.Header ) where
 
 import Control.Monad
-import System.IO
-
--- | The data type representing an analysis frame.
-data Frame = Frame { 
-    -- | Onset time of frame (seconds).
-    onset_time :: Double
-    -- | Duration of frame (seconds).
-    , duration :: Double
-    -- | Mean spectrum converted to the perceptually weighted Mel
-    --   frequency scale.
-    , avg_mel_spec :: [Double]
-    -- | Vector of energy distribution across each semitone of the
-    --   octave.
-    , avg_chroma :: [Double]
-    -- | Average power.
-    , avg_chunk_power :: Double
-    -- | Dominant semitone within the octave.
-    , avg_chroma_scalar :: Double
-    -- | A frequency estimation.
-    , avg_freq_simple :: Double
-    -- | Mean MFCCs, a commonly used feature in speech recognition.
-    , avg_mfcc :: [Double]
-    -- | Pitch estimation.
-    , avg_pitch_simple :: Double
-    -- | Average spectral center of mass of the frame.
-    , avg_spec_centroid :: Double
-    -- | Measure of flatness of the average spectrum.
-    , avg_spec_flatness :: Double
-    -- | Mean spectrum.
-    , avg_spec :: [Double]
-    -- | Length of analysis frame (sample frames).
-    , chunk_length :: Double
-    -- | Start time of analysis frame (sample frames).
-    , chunk_start_time :: Double
-    -- | The stability of the spectral energy.  More spectrally stable
-    --   frames are more likely to be pitched material.
-    , spectral_stability :: Double 
-    , raw_data :: (Double, Double, [[Double]]) }
-             deriving (Show, Read)
-
--- | A segment is a pair (onset_time, duration).
-type Segment = (Double, Double)
-
--- | Extract frame timing information.
-segments :: [Frame] -> [Segment]
-segments = map (\c -> (onset_time c, duration c))
-
--- | As association list indicating the number of data values
---   associated with each feature type.
-features :: [(String, Int)]
-features = [ ("AvgMelSpec", 40)
-           , ("AvgChroma", 12)
-           , ("AvgFramePower", 1)
-           , ("AvgChromaScalar", 1)
-           , ("AvgFreqSimple", 1)
-           , ("AvgMFCC", 13)
-           , ("AvgPitchSimple", 1)
-           , ("AvgSpecCentroid", 1)
-           , ("AvgSpecFlatness", 1)
-           , ("AvgSpec", 513)
-           , ("FrameLength", 1)
-           , ("FrameStartTime", 1)
-           , ("SpectralStability", 1) ]
-
--- | Read a MEAPsoft analysis file.
-read_file :: FilePath -> IO [Frame]
-read_file f =
-    do h <- openFile f ReadMode
-       t <- hGetContents h
-       let l = drop 2 (lines t)
-       return (map parse_frame l)
+import Data.Array.Unboxed
+import Data.List
+import Sound.Analysis.Meapsoft.Data
+import Sound.Analysis.Meapsoft.Header
 
-form_groups :: [Int] -> [a] -> [[a]]
-form_groups [] _ = []
-form_groups (n:ns) l = a : form_groups ns b
-    where (a, b) = splitAt n l
+-- | Data type representing a MEAPsoft anaylsis file.
+data MEAP = MEAP { -- | The list of 'Feature's contained in the
+                   --   analysis file.
+                   features :: [Feature]
+                   -- | The number of frames (rows) contained in the
+                   --   analysis file.
+                 , n_frames :: Int
+                   -- | The analysis data stored in a 'UArray'.
+                   --   Indices are of the form (row, column).
+                 , uarray_data :: UArray (Int, Int) Double }
+          
+-- | Load a MEAPsoft analysis file.
+read_meap :: FilePath -> IO (Either String MEAP)
+read_meap fn = do
+  r <- read_header fn
+  case r of
+    (Left e) -> return (Left e)
+    (Right h) -> do let nc = sum (map feature_degree h)
+                    (nf, d) <- read_data fn nc
+                    return (Right (MEAP h nf d))
 
-type Chunk = (Double, Double, [[Double]])
+-- | The number of columns at each analysis frame (row).
+n_columns :: MEAP -> Int
+n_columns = sum . map feature_degree . features
 
-make_chunk :: [Double] -> Chunk
-make_chunk (o:t:f) = (o, t, form_groups (map snd features) f)
-make_chunk _ = error "misformed data"
+-- | Extract the indicated column as a list.
+column_l :: MEAP -> Int -> [Double]
+column_l m j =
+    let d = uarray_data m
+        is = [0 .. n_frames m - 1]
+    in map (\i -> d ! (i, j)) is
 
-make_frame :: Chunk -> Frame
-make_frame c@(o, t, [ f1, f2, [f3], [f4], [f5]
-                    , f6, [f7], [f8], [f9], f10
-                    , [f11], [f12], [f13] ]) =
-    Frame o t f1 f2 f3 f4 f5 f6 f7 f8 f9 f10 f11 f12 f13 c
-make_frame _ = error "make_frame"
+-- | Extract the indicated frame (row) as a list.
+frame_l :: MEAP -> Int -> [Double]
+frame_l m i =
+    let d = uarray_data m
+        js = [0 .. n_columns m - 1]
+    in map (\j -> d ! (i, j)) js
 
-split_on :: (Eq a) => a -> [a] -> [[a]]
-split_on x xs = l' ++ r'
-    where (l, r) = span (/= x) xs
-          l' = if null l then [] else [l]
-          r' = if null r then [] else split_on x (tail r)
+-- | Extract data from the indicated frame and column.
+position :: MEAP -> (Int, Int) -> Double
+position m (i, j) = uarray_data m ! (i, j)
 
-to_f :: String -> Double
-to_f = read
+-- | Extract segmentation data as a list.  The segmentation data is
+--   given by the two columns onset_time and chunk_length.
+segments_l :: MEAP -> [(Double, Double)]
+segments_l m =
+    let fs = features m
+        ec n = feature_column (required_feature n fs)
+        o = ec "onset_time"
+        d = ec "chunk_length"
+    in zip (column_l m o) (column_l m d)
 
-parse_frame :: String -> Frame
-parse_frame = make_frame . make_chunk . map to_f . drop 1 . split_on ' '
+-- | The list of feature names generated by MEAPsoft 2.0.
+feature_names :: [String]
+feature_names =
+    ["onset_time"
+    ,"chunk_length"
+    ,"AvgChroma"
+    ,"AvgChromaScalar"
+    ,"AvgChunkPower"
+    ,"AvgFreqSimple"
+    ,"AvgMelSpec"
+    ,"AvgMFCC"
+    ,"AvgPitch"
+    ,"AvgSpec"
+    ,"AvgSpecCentroid"
+    ,"AvgSpecFlatness"
+    ,"AvgTonalCentroid"
+    ,"ChunkLength"
+    ,"ChunkStartTime"
+    ,"Entropy"
+    ,"RMSAmplitude"
+    ,"SpectralStability"]
diff --git a/Sound/Analysis/Meapsoft/Data.hs b/Sound/Analysis/Meapsoft/Data.hs
new file mode 100644
--- /dev/null
+++ b/Sound/Analysis/Meapsoft/Data.hs
@@ -0,0 +1,22 @@
+module Sound.Analysis.Meapsoft.Data ( read_data ) where
+
+import Data.Array.Unboxed
+import qualified Data.ByteString as B
+import qualified Data.ByteString.Lazy as L
+import qualified Data.ByteString.Lex.Double as B
+import qualified Text.Delimited as D
+
+to_f :: B.ByteString -> Double
+to_f s = fst (maybe undefined id (B.readDouble s))
+
+-- | Given the number of columns, reads an entire MEAPsoft data set
+--   into an 'UArray' and returns the data paired with the number of
+--   rows.
+read_data :: FilePath -> Int -> IO (Int, UArray (Int, Int) Double)
+read_data fn nc = do
+  s <- L.readFile fn
+  let (Right xs) = D.decode " \t" s
+      xs' = map (take nc . drop 1) (drop 1 xs)
+      nr = length xs'
+      ds = map to_f (concat xs')
+  return (nr, listArray ((0,0), (nr - 1, nc - 1)) ds)
diff --git a/Sound/Analysis/Meapsoft/Header.hs b/Sound/Analysis/Meapsoft/Header.hs
new file mode 100644
--- /dev/null
+++ b/Sound/Analysis/Meapsoft/Header.hs
@@ -0,0 +1,97 @@
+module Sound.Analysis.Meapsoft.Header ( Feature(..)
+                                      , read_header
+                                      , find_feature
+                                      , required_feature
+                                      , has_feature
+                                      ) where
+
+import Control.Monad
+import Data.List
+import System.IO
+import Text.ParserCombinators.Parsec
+
+-- | Data type representing a MEAPsoft analysis feature.  The
+--   'feature_column' is the integer column index for the feature in
+--   the analysis data.  The 'feature_degree' is the number of columns
+--   the feature requires.
+data Feature = Feature { feature_name :: String
+                       , feature_column :: Int
+                       , feature_degree :: Int }
+               deriving (Show)
+
+-- | Read the header of a MEAPsoft analysis file and extract the list
+--   of stored features.
+read_header :: FilePath -> IO (Either String [Feature])
+read_header fn = do
+  s <- read_header_string fn
+  let r = parse_header fn s
+  return (case r of
+            (Right h) -> Right (mk_features (normalize_header h))
+            (Left e) -> Left (show e))
+
+-- | Search for a named feature.
+find_feature :: String -> [Feature] -> Maybe Feature
+find_feature n = find (\x -> feature_name x == n)
+
+-- | A variant of 'find_feature' that raises an error if the feature
+--   is not located.  All analysis files have the features onset_time
+--   and chunk_length.
+required_feature :: String -> [Feature] -> Feature
+required_feature n fs = maybe (error n) id (find_feature n fs)
+
+-- | True iff the analysis data contains the named feature.
+has_feature :: String -> [Feature] -> Bool
+has_feature n fs = maybe False (const True) (find_feature n fs)
+
+-- Parsec parser for header string.
+
+type P a = GenParser Char () a
+
+word :: P String
+word = many1 (letter <|> oneOf "_") <?> "word"
+
+whitespace :: P String
+whitespace = many1 (oneOf " \t")
+
+in_paren :: P a -> P a
+in_paren p =
+    do { char '('
+       ; r <- p
+       ; char ')'
+       ; return r }
+
+int :: P Int
+int = liftM read (optional (char '-') >> many1 digit)
+
+feature :: P (String, Int)
+feature =
+  do { f <- word
+     ; n <- optionMaybe (try (in_paren int))
+     ; return (f, maybe 1 id n) }
+
+type Header = [(String, Int)]
+
+features :: P Header
+features = sepEndBy1 feature whitespace
+
+hash :: P Char
+hash = char '#'
+
+header :: P Header
+header = hash >> features
+
+read_header_string :: String -> IO String
+read_header_string fn = withFile fn ReadMode hGetLine
+
+parse_header :: String -> String -> Either ParseError Header
+parse_header fn s = parse header fn s
+
+-- Delete 'filename', which is a string.
+normalize_header :: Header -> Header
+normalize_header (("filename",1):xs) = xs
+normalize_header xs = xs
+
+mk_features :: Header -> [Feature]
+mk_features h =
+    let acc i (f, n) = (i+n, (Feature f i n))
+    in snd (mapAccumL acc 0 h)
diff --git a/hmeap.cabal b/hmeap.cabal
--- a/hmeap.cabal
+++ b/hmeap.cabal
@@ -1,27 +1,37 @@
 Name:              hmeap
-Version:           0.6
+Version:           0.7
 Synopsis:          Haskell Meapsoft Parser
-Description:       Parser for the anaylsis files produced by the Meapsoft 
-                   feature extractor.
+Description:       Parser for the anaylsis files produced by the
+                   Meapsoft feature extractor.
 License:           GPL
 Category:          Sound
-Copyright:         Rohan Drape, 2007-2008
+Copyright:         Rohan Drape, 2007-2009
 Author:            Rohan Drape
 Maintainer:        rd@slavepianos.org
 Stability:         Experimental
 Homepage:          http://slavepianos.org/rd/f/896678/
-Tested-With:       GHC==6.10.1
+Tested-With:       GHC == 6.8.2
 Build-Type:        Simple
-Cabal-Version:     >= 1.2
+Cabal-Version:     >= 1.6
 
 Data-files:        README
                    Help/browse.hs
                    Help/Makefile
-                   Help/parse.hs
+                   Help/parser.hs
                    Help/play.hs
+                   Help/plot.hs
+                   Help/stat.hs
+                   sh/extract.sh
 
 Library
-  Build-Depends:   base
+  Build-Depends:   array,
+                   base == 3.*,
+                   bytestring,
+                   bytestring-lexing,
+                   delimited-text >= 0.1.0,
+                   parsec
   GHC-Options:     -Wall -fwarn-tabs
   Exposed-modules: Sound.Analysis.Meapsoft
                    Sound.Analysis.Meapsoft.Measure
+                   Sound.Analysis.Meapsoft.Header
+                   Sound.Analysis.Meapsoft.Data
diff --git a/sh/extract.sh b/sh/extract.sh
new file mode 100644
--- /dev/null
+++ b/sh/extract.sh
@@ -0,0 +1,58 @@
+#!/bin/sh
+
+if test $# != 2
+then
+    echo "extract.sh meapsoft-jar audio-file"
+    exit 0
+fi
+
+if test -d $1
+then
+  echo "meapsoft dir: " $1
+else
+  echo "non-existant meapsoft dir: " $1
+  exit 1
+fi
+
+echo "file to analyse: " $2 $2.seg $2.feat
+
+if test -f $2.seg
+then
+    echo "$2.seg exists"
+else
+    java \
+        -cp $1 \
+        com.meapsoft.Segmenter \
+        -o $2.seg \
+        $2 \
+        -t 1.0 \
+        -s 0.1 \
+        -d
+fi
+
+# Xms = initial heap size
+# Xmx = maximum heap size
+
+java \
+    -Xms128m \
+    -Xmx512m \
+    -cp $1 \
+    com.meapsoft.FeatExtractor \
+    -f AvgChroma \
+    -f AvgChromaScalar \
+    -f AvgChunkPower \
+    -f AvgFreqSimple \
+    -f AvgMelSpec \
+    -f AvgMFCC \
+    -f AvgPitch \
+    -f AvgSpec \
+    -f AvgSpecCentroid \
+    -f AvgSpecFlatness \
+    -f AvgTonalCentroid \
+    -f ChunkLength \
+    -f ChunkStartTime \
+    -f Entropy \
+    -f RMSAmplitude \
+    -f SpectralStability \
+    $2.seg \
+    -o $2.feat