diff --git a/CHANGES b/CHANGES
--- a/CHANGES
+++ b/CHANGES
@@ -21,3 +21,12 @@
 		added mapArrayConcurrently
 		improved performance of BDF readData
 		changed BDF to store data as Floats, not Double
+
+0.1.2.1:
+		fixed mapConcurrently to wait for all threads
+		fixed Signal.detrend window sizes 
+		fixed Signal.detrend mathematics
+		fixed broadband_filter to call the correct filter (Numeric.Signal)
+		fixed signal-aux.c filter scaling error
+		fixed hilbert
+		changed Multichannel.fromList to createMultichannel to avoid name conflicts
diff --git a/hsignal.cabal b/hsignal.cabal
--- a/hsignal.cabal
+++ b/hsignal.cabal
@@ -1,5 +1,5 @@
 Name:               hsignal
-Version:            0.1.1.1
+Version:            0.1.2.1
 License:            GPL
 License-file:       LICENSE
 Author:             Vivian McPhail
@@ -10,7 +10,7 @@
 Description:        Purely functional interface to signal processing based on hmatrix
                        provides data types for manipulating EEG data,
                        including reading from BDF
-Category:           Math
+Category:           Math, 
 tested-with:        GHC ==6.12.1
 
 cabal-version:      >=1.8
@@ -27,7 +27,8 @@
                         array,
                         bytestring, storable-complex, ghc-binary,
                         hmatrix >= 0.10.0,
-                        hmatrix-gsl-stats >= 0.1.1.4
+                        hmatrix-gsl-stats >= 0.1.1.5,
+                        hstatistics >= 0.2.0.5
 
     Extensions:         ForeignFunctionInterface
 
diff --git a/lib/Numeric/Signal.hs b/lib/Numeric/Signal.hs
--- a/lib/Numeric/Signal.hs
+++ b/lib/Numeric/Signal.hs
@@ -29,7 +29,7 @@
 
 import Numeric.GSL.Fitting.Linear
  
-import Complex
+import Data.Complex
 
 import qualified Data.List as L
 
@@ -53,9 +53,12 @@
        -> Int           -- ^ sampling rate
        -> Vector Double -- ^ input signal
        -> Vector Double -- ^ output signal
-filter b a s v = let sd = (fromIntegral s) * 2.0
-                     sc = recip $ (sd / (fromIntegral $ max (dim b) (dim a))) * sd
-                     in scale sc $ S.filter b a v
+filter b a s v = let len = dim v
+                     w = min s len
+                     start = (negate . fromList . reverse . toList . subVector 0 w) v
+                     finish = (negate . fromList . reverse . toList . subVector (len-w) w) v
+                     v' = join [start,v,finish]
+                 in subVector s len $ S.filter b a v'
 
 -----------------------------------------------------------------------------
                      
@@ -100,7 +103,7 @@
                  -> Vector Double -- ^ input signal
                  -> Vector Double -- ^ output signal
 broadband_filter s f v = let b = broadband_fir s f
-                         in S.filter b (constant 1.0 1) v
+                         in filter b (scalar 1.0) s v
                                 
 -----------------------------------------------------------------------------
 
@@ -208,15 +211,15 @@
         -> Vector Double   -- ^ data to be detrended
         -> Vector Double   -- ^ detrended data
 detrend w v = let windows = dim v `div` w
-                  ws = takesV ((replicate windows w) ++ [dim v - (windows * w) + 1]) v
+                  ws = takesV ((replicate windows w) ++ [dim v - (windows * w)]) v
                   ds = map detrend' ws
                   windows' = (dim v - (w `div` 2)) `div` w
-                  ws' = takesV (((w `div` 2):(replicate windows w)) ++ [dim v - (w `div` 2) - (windows' * w) + 1]) v
+                  ws' = takesV (((w `div` 2):(replicate windows' w)) ++ [dim v - (w `div` 2) - (windows' * w)]) v
                   ds' = map detrend' ws'
-              in (join ds + join ds') / 2
+              in (join ds + join ds') / 2 
     where detrend' x = let ln = dim x
-                           y = linspace ln (1.0,fromIntegral ln)
-                           (c0,c1,_,_,_,_) = linear x y
-                       in x - (fromList [c0]) - (scale c1 x)
+                           t = linspace ln (1.0,fromIntegral ln)
+                           (c0,c1,_,_,_,_) = linear t x
+                       in x - (scale c1 t + scalar c0)
 
 -----------------------------------------------------------------------------
diff --git a/lib/Numeric/Signal/EEG/BDF.hs b/lib/Numeric/Signal/EEG/BDF.hs
--- a/lib/Numeric/Signal/EEG/BDF.hs
+++ b/lib/Numeric/Signal/EEG/BDF.hs
@@ -279,7 +279,7 @@
              (bdf,bs') <- runStateT readBDF bs
              m <- case bdf of
                            (Just b) -> do
-                                       return $ Just $ M.fromList ((head $ samples b) * (num_records b)) 24 (data_ b)
+                                       return $ Just $ M.createMultichannel (head $ samples b) 24 (data_ b)
                            _        -> do
                                        putStrLn "File not read"
                                        return Nothing
diff --git a/lib/Numeric/Signal/Internal.hs b/lib/Numeric/Signal/Internal.hs
--- a/lib/Numeric/Signal/Internal.hs
+++ b/lib/Numeric/Signal/Internal.hs
@@ -92,13 +92,12 @@
 -- | Hilbert transform with original vector as real value, transformed as imaginary
 hilbert :: Vector Double -> Vector (Complex Double)
 hilbert v = unsafePerformIO $ do
-            r <- createVector (dim v)
-            let v' = complex v
+            let r = complex v
             -- could use (comp v) to make a complex vector in haskell rather than C
-            app2 signal_hilbert vec v' vec r "hilbert"
+            app1 signal_hilbert vec r "hilbert"
             return r
 
-foreign import ccall "signal-aux.h hilbert" signal_hilbert :: CInt -> PC -> CInt -> PC -> IO CInt
+foreign import ccall "signal-aux.h hilbert" signal_hilbert :: CInt -> PC -> IO CInt
 
 -----------------------------------------------------------------------------
 
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
@@ -18,13 +18,14 @@
 
 module Numeric.Signal.Multichannel (
                        Multichannel,readMultichannel,writeMultichannel,
-                       fromList,
+                       createMultichannel,
                        sampling_rate,precision,channels,samples,
                        detrended,filtered,
                        getChannel,getChannels,
                        mapConcurrently,
                        detrend,filter,
-                       slice
+                       slice,
+                       mi_phase
                 ) where
 
 -----------------------------------------------------------------------------
@@ -34,6 +35,7 @@
 --import Complex
 
 import qualified Data.Array.IArray as I
+import Data.Ix
 
 import Control.Concurrent
 --import Control.Concurrent.MVar
@@ -42,7 +44,8 @@
 
 --import qualified Data.List as L
 
-import Data.Packed.Vector hiding(fromList)
+import Data.Packed.Vector
+import Data.Packed.Matrix
 --import Data.Packed(Container(..))
 
 import Data.Binary
@@ -54,9 +57,14 @@
 
 --import qualified Numeric.GSL.Fourier as F
 
+import qualified Numeric.GSL.Histogram as H
+import qualified Numeric.GSL.Histogram2D as H2
+
+import qualified Numeric.Statistics.Information as SI
+
 import Prelude hiding(filter)
 
---import Control.Monad
+import Control.Monad(replicateM)
 
 
 {-
@@ -145,11 +153,11 @@
 -----------------------------------------------------------------------------
 
 -- | create a multichannel data type
-fromList :: Storable a => Int -- ^ sampling rate
-         -> Int               -- ^ bits of precision
-         -> [Vector a]        -- ^ data
-         -> Multichannel a    -- ^ datatype
-fromList s p d = let c = length d
+createMultichannel :: Storable a => Int -- ^ sampling rate
+                   -> Int               -- ^ bits of precision
+                   -> [Vector a]      -- ^ data
+                   -> Multichannel a    -- ^ datatype
+createMultichannel s p d = let c = length d
                  in MC s p c (dim $ head d) False Nothing (I.listArray (1,c) d)
 
 -- | the sampling rate
@@ -187,19 +195,30 @@
 -----------------------------------------------------------------------------
 
 -- | map a function executed concurrently
-mapArrayConcurrently :: (a -> b)            -- ^ function to map
-                     -> I.Array Int a       -- ^ input
-                     -> I.Array Int b       -- ^ output
+mapArrayConcurrently :: Ix i => (a -> b)    -- ^ function to map
+                     -> I.Array i a         -- ^ input
+                     -> I.Array i b         -- ^ output
 mapArrayConcurrently f d = unsafePerformIO $ do
-                                             results <- newMVar []
-                                             mapM_ (forkIO . applyFunction results f) $ I.assocs d
-                                             vectors <- takeMVar results
+                                             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
-    where applyFunction results f' (j,e) = do
-                                          let o = f' e
-                                          modifyMVar_ results (\x -> return ((j,o):x))
+    where applyFunction f' (m,(j,e)) = putMVar m (j,f' e)
 
+{-
 -- | map a function executed concurrently
+mapListConcurrently :: (a -> b)             -- ^ function to map
+                    -> [a]                  -- ^ input
+                    -> [b]                  -- ^ output
+mapListConcurrently f d = unsafePerformIO $ do
+                                            results <- replicateM (length d) newEmptyMVar
+                                            mapM_ (forkIO . applyFunction f) zip results d
+                                            mapM takeMVar results
+    where applyFunction f' (m,e) = putMVar m (f' e)
+-}
+
+-- | map a function executed concurrently
 mapConcurrently :: Storable b => (Vector a -> Vector b)  -- ^ the function to be mapped 
                 -> Multichannel a                        -- ^ input data
                 -> Multichannel b                        -- ^ output data
@@ -229,11 +248,34 @@
 -----------------------------------------------------------------------------
 
 -- | extract a slice of the data
-slice :: Int                 -- ^ starting sample number
+slice :: Storable a => Int   -- ^ starting sample number
       -> Int                 -- ^ length
-      -> Multichannel Double 
-      -> Multichannel Double
+      -> Multichannel a 
+      -> Multichannel a
 slice j w m = let m' = mapConcurrently (subVector j w) m
               in m' { _length = w }
+
+-----------------------------------------------------------------------------
+
+-- | calculate the mutual information of the phase between pairs of channels
+-- |      the lower half of the matrix displays mutual information between pairs of channels
+-- |      the upper half of the matrix displays mutual information between the data (row) and a sample
+-- |          from the same distribution as the column
+mi_phase :: Multichannel Double      -- ^ input data
+         -> Matrix Double
+mi_phase m = let d = _data m
+                 s = samples m
+                 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)))
+                 hist2array = mapArrayConcurrently (\(j,x,y) -> (j,H2.addVector (H2.emptyLimits 128 128 (-pi,pi) (-pi,pi)) x y)) pairs
+                 mi = mapArrayConcurrently (doMI histarray d s) hist2array
+             in fromArray2D mi
+    where doMI histarray d s ((x,y),h2) 
+              | x < y     = SI.mutual_information h2 (histarray I.! x) (histarray I.! y) (d I.! x,d I.! y)
+              | otherwise = let hypdf = H.fromHistogram (histarray I.! y)
+                                ys = replicate s (H.sample hypdf) 
+                                dy' = fromList ys
+                            in SI.mutual_information h2 (histarray I.! x) (histarray I.! y) (d I.! x,dy')
 
 -----------------------------------------------------------------------------
diff --git a/lib/Numeric/Signal/signal-aux.c b/lib/Numeric/Signal/signal-aux.c
--- a/lib/Numeric/Signal/signal-aux.c
+++ b/lib/Numeric/Signal/signal-aux.c
@@ -2,6 +2,7 @@
 
 #include <gsl/gsl_math.h>
 #include <gsl/gsl_fft_real.h>
+#include <gsl/gsl_fft_halfcomplex.h>
 #include <gsl/gsl_fft_complex.h>
 #include <gsl/gsl_vector.h>
 #include <gsl/gsl_blas.h>
@@ -59,21 +60,19 @@
 
   for (i = 0; i < vs; i++) {
     r[i] = 0;
-    for (j = 1; j < N; j++) {
-      if (i - j > 0) r[i] -= (l[j+1]/L)*v[i-j];
+    for (j = 0; j < N; j++) {
+      if (i - j > 0) r[i] -= (l[j+1])*v[i-j];
     }
     for (j = 0; j < M; j++) {
-      if (i - j > 0) r[i] += (k[j+1]/K)*r[i-j];
+      if (i - j > 0) r[i] += (k[j+1])*r[i-j];
     }
   }
   return 0;
 }
 
-int hilbert(int vs, const gsl_complex* v, int rs, gsl_complex* r)
+int hilbert(int rs, gsl_complex* r)
 {
-  if (vs != rs) return 2000; // BAD_SIZE
-
-  int s = vs;
+  int s = rs;
 
   gsl_fft_complex_wavetable * wavetable = gsl_fft_complex_wavetable_alloc (s);
   gsl_fft_complex_workspace * workspace = gsl_fft_complex_workspace_alloc (s);
@@ -81,17 +80,22 @@
   // forward fourier transform
   gsl_fft_complex_forward ((double*)r, 1, s, wavetable, workspace);
   // zero negative coefficients and double positive
+
   int i;
-  for (i = 0; i < s; i++) {
-    if (i <= s/2) {
-      r[i].dat[0] *= sqrt(2.0);
-      r[i].dat[1] *= sqrt(2.0);
+  int m = s/2;
+  for (i = 1; i < s; i++) {
+    if (i <= m) {
+      r[i].dat[0] *= 2;
+      r[i].dat[1] *= 2;
     }
+    else if (s % 2 == 0 && i == m+1) {
+    }
     else {
-      r[i].dat[0] = 0.0;
-      r[i].dat[1] = 0.0;
+      r[i].dat[0] = 0;
+      r[i].dat[1] = 0;
     }
   }
+
   // inverse fourier transform
   gsl_fft_complex_inverse ((double*)r, 1, s, wavetable, workspace);