emd 0.1.6.0 → 0.1.7.0

raw patch · 4 files changed

+101/−33 lines, 4 filesdep +arraydep +carraydep +fftPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies added: array, carray, fft

API changes (from Hackage documentation)

Files

• CHANGELOG.md +10/−0
• emd.cabal +7/−3
• src/Numeric/HHT.hs +39/−30
• src/Numeric/HHT/Internal/FFT.hs +45/−0

@@ -1,6 +1,16 @@ Changelog ========= +Version 0.1.7.0+---------------++*September 24, 2019*++<https://github.com/mstksg/emd/releases/tag/v0.1.7.0>++*   Rewrite `hilbert` using the *fft* library, matching the matlab+    implementation.  This means that the library now depends on *fftw*.+ Version 0.1.6.0 --------------- 

@@ -4,10 +4,10 @@ -- -- see: https://github.com/sol/hpack ----- hash: ff41309ce78de40674bacd4a2d67e4f7be21c4e08e289bc0faff8f38d771e3fc+-- hash: 9e35255466a474e11d7cfdc5ca94ab126753f2fe938c00366f389a8ee69a989f  name:           emd-version:        0.1.6.0+version:        0.1.7.0 synopsis:       Empirical Mode Decomposition and Hilbert-Huang Transform description:    Empirical Mode decomposition and Hilbert-Huang Transform in pure                 Haskell.@@ -38,15 +38,19 @@   other-modules:       Numeric.EMD.Internal.Tridiagonal       Numeric.EMD.Internal.Extrema+      Numeric.HHT.Internal.FFT   hs-source-dirs:       src   ghc-options: -Wall -Wredundant-constraints -Wcompat   build-depends:-      base >=4.10 && <5+      array+    , base >=4.10 && <5     , binary+    , carray     , containers     , data-default-class     , deepseq+    , fft     , finite-typelits     , ghc-typelits-knownnat     , ghc-typelits-natnormalise

@@ -2,10 +2,13 @@ {-# LANGUAGE DataKinds                                #-} {-# LANGUAGE DeriveGeneric                            #-} {-# LANGUAGE FlexibleContexts                         #-}+{-# LANGUAGE MultiWayIf                               #-} {-# LANGUAGE RecordWildCards                          #-} {-# LANGUAGE ScopedTypeVariables                      #-} {-# LANGUAGE TypeApplications                         #-}+{-# LANGUAGE TypeFamilies                             #-} {-# LANGUAGE TypeOperators                            #-}+{-# LANGUAGE ViewPatterns                             #-} {-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-} {-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise       #-} @@ -24,12 +27,6 @@ -- 'hhtEmd'.  See "Numeric.EMD" for information on why this module uses -- "sized vectors", and how to convert unsized vectors to sized vectors. ----- Note that the Hilbert Transform implementation in this module is--- slightly naive and is essentially O(n^2) on the length of the vector.--- However, computation time for the full Hilbert-Huang Transform is--- typically dominated by Empirical Mode Docomposition, which is--- approximately O(n).--- -- @since 0.1.2.0  module Numeric.HHT (@@ -62,12 +59,14 @@ import           GHC.Generics              (Generic) import           GHC.TypeNats import           Numeric.EMD+import           Numeric.HHT.Internal.FFT import qualified Data.Binary               as Bi import qualified Data.List.NonEmpty        as NE import qualified Data.Map                  as M import qualified Data.Vector.Generic       as VG import qualified Data.Vector.Generic.Sized as SVG import qualified Data.Vector.Sized         as SV+import qualified Math.FFT.Base             as FFT  -- | A Hilbert Trasnform of a given IMF, given as a "skeleton line". data HHTLine v n a = HHTLine@@ -107,7 +106,7 @@ -- | Directly compute the Hilbert-Huang transform of a given time series. -- Essentially is a composition of 'hhtEmd' and 'emd'.  See 'hhtEmd' for -- a more flexible version.-hht :: forall v n a. (VG.Vector v a, KnownNat n, RealFloat a)+hht :: forall v n a. (VG.Vector v a, VG.Vector v (Complex a), KnownNat n, FFT.FFTWReal a)     => EMDOpts a     -> SVG.Vector v (n + 1) a     -> HHT v n a@@ -116,7 +115,7 @@ -- | Compute the Hilbert-Huang transform from a given Empirical Mode -- Decomposition. hhtEmd-    :: forall v n a. (VG.Vector v a, KnownNat n, RealFloat a)+    :: forall v n a. (VG.Vector v a, VG.Vector v (Complex a), KnownNat n, FFT.FFTWReal a)     => EMD v (n + 1) a     -> HHT v n a hhtEmd EMD{..} = HHT $ map go emdIMFs@@ -277,14 +276,14 @@ -- anything faster given the discretization, and we exclude negative values -- as physically unmeaningful for an IMF. hilbertMagFreq-    :: forall v n a. (VG.Vector v a, KnownNat n, RealFloat a)+    :: forall v n a. (VG.Vector v a, VG.Vector v (Complex a), KnownNat n, FFT.FFTWReal a)     => SVG.Vector v (n + 1) a     -> (SVG.Vector v (n + 1) a, SVG.Vector v n a) hilbertMagFreq v = (hilbertMag, hilbertFreq)   where-    v'           = hilbertIm v-    hilbertMag   = SVG.zipWith (\x x' -> magnitude (x :+ x')) v v'-    hilbertPhase = SVG.zipWith (\x x' -> phase (x :+ x')) v v'+    v'           = hilbert v+    hilbertMag   = SVG.map magnitude v'+    hilbertPhase = SVG.map phase v'     hilbertFreq  = SVG.map ((`mod'` 1) . (/ (2 * pi))) $ SVG.tail hilbertPhase - SVG.init hilbertPhase  -- | The polar form of 'hilbert': returns the magnitude and phase of the@@ -304,7 +303,7 @@ -- -- @since 0.1.6.0 hilbertPolar-    :: forall v n a. (VG.Vector v a, KnownNat n, RealFloat a)+    :: forall v n a. (VG.Vector v a, VG.Vector v (Complex a), KnownNat n, FFT.FFTWReal a)     => SVG.Vector v (n + 1) a     -> (SVG.Vector v (n + 1) a, SVG.Vector v (n + 1) a) hilbertPolar v = (hilbertMag, hilbertPhase)@@ -320,31 +319,41 @@ -- series.  Creates a "helical" form of the original series that rotates -- along the complex plane. ----- Numerically assumes that the signal is zero everywhere outside of the--- vector, instead of the periodic assumption taken by matlab's version.+-- Note that since /0.1.7.0/, this uses the same algorithm as the matlab+-- implementation <https://www.mathworks.com/help/signal/ref/hilbert.html> hilbert-    :: forall v n a. (VG.Vector v a, VG.Vector v (Complex a), KnownNat n, Floating a)+    :: forall v n a.+      ( VG.Vector v a+      , VG.Vector v (Complex a)+      , KnownNat n+      , FFT.FFTWReal a+      )     => SVG.Vector v n a     -> SVG.Vector v n (Complex a)-hilbert v = SVG.zipWith (:+) v (hilbertIm v)+hilbert v = ifft u'+  where+    v' = SVG.map (:+ 0) v+    u  = fft v'+    u' = flip SVG.imap u $ \(fromIntegral->i) x ->+      if | i == 0 || i == (n `div` 2) -> x+         | i < (n `div` 2)            -> 2 * x+         | otherwise                  -> 0+    n  = natVal (Proxy @n)  -- | Hilbert transformed series.  Essentially the same series, but -- phase-shifted 90 degrees.  Is so-named because it is the "imaginary -- part" of the proper hilbert transform, 'hilbert'. ----- Numerically assumes that the signal is zero everywhere outside of the--- vector, instead of the periodic assumption taken by matlab's version.+-- Note that since /0.1.7.0/, this uses the same algorithm as the matlab+-- implementation <https://www.mathworks.com/help/signal/ref/hilbert.html> hilbertIm-    :: forall v n a. (VG.Vector v a, KnownNat n, Floating a)+    :: forall v n a.+      ( VG.Vector v a+      , VG.Vector v (Complex a)+      , KnownNat n+      , FFT.FFTWReal a+      )     => SVG.Vector v n a     -> SVG.Vector v n a-hilbertIm v = SVG.generate $ \i -> getSum . foldMap (Sum . go i) $ finites @n-  where-    -- NOTE: Can be made faster using an FFT and iFFT combo-    go :: Finite n -> Finite n -> a-    go i j-        | even k    = 0-        | otherwise = 2 * (v `SVG.index` j) / pi / fromIntegral k-      where-        k :: Int-        k = fromIntegral i - fromIntegral j+hilbertIm = SVG.map imagPart . hilbert+

@@ -0,0 +1,45 @@+{-# LANGUAGE FlexibleContexts #-}++module Numeric.HHT.Internal.FFT (+    fft+  , ifft+  ) where++import           Data.Complex+import qualified Data.Array.CArray         as CA+import qualified Data.Array.IArray         as IA+import qualified Data.Ix                   as Ix+import qualified Data.Vector.Generic       as VG+import qualified Data.Vector.Generic.Sized as SVG+import qualified Foreign.Storable          as FS+import qualified Math.FFT                  as FFT+import qualified Math.FFT.Base             as FFT++fft :: (FFT.FFTWReal a, VG.Vector v (Complex a))+    => SVG.Vector v n (Complex a)+    -> SVG.Vector v n (Complex a)+fft = SVG.withVectorUnsafe $+        fromCA+      . FFT.dft+      . toCA++ifft+    :: (FFT.FFTWReal a, VG.Vector v (Complex a))+    => SVG.Vector v n (Complex a)+    -> SVG.Vector v n (Complex a)+ifft = SVG.withVectorUnsafe $+        fromCA+      . FFT.idft+      . toCA++fromCA+    :: (FS.Storable a, VG.Vector v (Complex a))+    => CA.CArray Int (Complex a)+    -> v (Complex a)+fromCA v = VG.generate (Ix.rangeSize (IA.bounds v)) (v IA.!)++toCA+    :: (FS.Storable a, VG.Vector v (Complex a))+    => v (Complex a)+    -> CA.CArray Int (Complex a)+toCA v = IA.listArray (0, VG.length v - 1) (VG.toList v)