diff --git a/CHANGES b/CHANGES
--- a/CHANGES
+++ b/CHANGES
@@ -69,3 +69,7 @@
 		renamed vector_deriv to vector_diff_double
 		added vector_diff_float
 
+0.1.3.6:
+		cleaned up mi_phase
+		generate H(X;X) for reference value
+
diff --git a/hsignal.cabal b/hsignal.cabal
--- a/hsignal.cabal
+++ b/hsignal.cabal
@@ -1,5 +1,5 @@
 Name:               hsignal
-Version:            0.1.3.5
+Version:            0.1.3.6
 License:            BSD3
 License-file:       LICENSE
 Copyright:          (c) A.V.H. McPhail 2010, 2011
diff --git a/lib/Numeric/Signal/Multichannel.hs b/lib/Numeric/Signal/Multichannel.hs
--- a/lib/Numeric/Signal/Multichannel.hs
+++ b/lib/Numeric/Signal/Multichannel.hs
@@ -191,13 +191,13 @@
 -- | map a function executed concurrently
 mapArrayConcurrently :: Ix i => (a -> b)    -- ^ function to map
                      -> I.Array i a         -- ^ input
-                     -> I.Array i b         -- ^ output
+                     -> (I.Array i b)     -- ^ output
 mapArrayConcurrently f d = unsafePerformIO $ do
-                                             let b = I.bounds d
-                                             results <- replicateM (rangeSize b) newEmptyMVar
-                                             mapM_ (forkIO . applyFunction f) $ zip results (I.assocs d)
-                                             vectors <- mapM takeMVar results
-                                             return $ I.array (I.bounds d) vectors
+  let b = I.bounds d
+  results <- replicateM (rangeSize b) newEmptyMVar
+  mapM_ (forkIO . applyFunction f) $ zip results (I.assocs d)
+  vectors <- mapM takeMVar results
+  return $ I.array b vectors
     where applyFunction f' (m,(j,e)) = putMVar m (j,f' e)
 
 {-
@@ -260,14 +260,15 @@
            Multichannel a      -- ^ input data
          -> Matrix Double
 mi_phase m = let d = fmap double $ _data m
+                 (bl,bu) = I.bounds d
+                 br = ((bl,bl),(bu,bu))
                  histarray = mapArrayConcurrently (H.fromLimits 128 (-pi,pi)) d
-                 c = channels m
-                 pairs = I.array ((1,1),(c,c)) $ map (\(a,b) -> ((a,b),((a,b),d I.! a,d I.! b))) (range ((1,1),(c,c)))
+                 pairs = I.array br $ map (\(a,b) -> ((a,b),((a,b),d I.! a,d I.! b))) (range br)
                  hist2array = mapArrayConcurrently (\(j,x,y) -> (j,H2.addVector (H2.emptyLimits 128 128 (-pi,pi) (-pi,pi)) x y)) pairs
                  mi = mapArrayConcurrently (doMI histarray d) hist2array
              in fromArray2D mi
     where doMI histarray d ((x,y),h2) 
-              | x < y     = SI.mutual_information h2 (histarray I.! x) (histarray I.! y) (d I.! x,d I.! y)
+              | x <= y     = SI.mutual_information h2 (histarray I.! x) (histarray I.! y) (d I.! x,d I.! y)
               | otherwise = 0
 
 -----------------------------------------------------------------------------