diff --git a/CHANGELOG.md b/CHANGELOG.md --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,6 +1,16 @@ Changelog ========= +Version 0.1.7.0 +--------------- + +*September 24, 2019* + + + +* Rewrite `hilbert` using the *fft* library, matching the matlab + implementation. This means that the library now depends on *fftw*. + Version 0.1.6.0 --------------- diff --git a/emd.cabal b/emd.cabal --- a/emd.cabal +++ b/emd.cabal @@ -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 diff --git a/src/Numeric/HHT.hs b/src/Numeric/HHT.hs --- a/src/Numeric/HHT.hs +++ b/src/Numeric/HHT.hs @@ -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 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 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 + diff --git a/src/Numeric/HHT/Internal/FFT.hs b/src/Numeric/HHT/Internal/FFT.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/HHT/Internal/FFT.hs @@ -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)