{-# LANGUAGE BangPatterns #-}
{-# 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 #-}
module Numeric.HHT (
HHT(..), HHTLine(..)
, hhtEmd
, hht
, ihhtEmd
, ihht
, hhtSpectrum, hhtSparseSpectrum, hhtDenseSpectrum
, meanMarginal, marginal, instantaneousEnergy, degreeOfStationarity
, expectedFreq, dominantFreq
, foldFreq
, EMDOpts(..), defaultEO, BoundaryHandler(..), SiftCondition(..), defaultSC, SplineEnd(..)
, hilbert
, hilbertIm
, hilbertPolar
, hilbertMagFreq
) where
import Control.DeepSeq
import Data.Complex
import Data.Finite
import Data.Fixed
import Data.Foldable
import Data.Proxy
import Data.Semigroup
import GHC.Generics (Generic)
import GHC.TypeNats
import Numeric.EMD
import Numeric.HHT.Internal.FFT
import qualified Data.Binary as Bi
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
data HHTLine v n a = HHTLine
{
hlMags :: !(SVG.Vector v (n + 1) a)
, hlFreqs :: !(SVG.Vector v n a)
, hlInitPhase :: !a
}
deriving (Show, Eq, Ord, Generic)
instance (VG.Vector v a, KnownNat n, Bi.Binary (v a), Bi.Binary a) => Bi.Binary (HHTLine v n a)
instance (NFData (v a), NFData a) => NFData (HHTLine v n a)
data HHT v n a = HHT
{
hhtLines :: [HHTLine v n a]
, hhtResidual :: SVG.Vector v (n + 1) a
}
deriving (Show, Eq, Ord, Generic)
instance (VG.Vector v a, KnownNat n, Bi.Binary (v a), Bi.Binary a) => Bi.Binary (HHT v n a)
instance (NFData (v a), NFData a) => NFData (HHT v n 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
hht eo = hhtEmd . emd eo
hhtEmd
:: 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) emdResidual
where
go i = HHTLine m f φ0
where
(m, (f, φ0)) = hilbertMagFreq i
ihhtEmd
:: (VG.Vector v a, Floating a)
=> HHT v n a
-> EMD v (n + 1) a
ihhtEmd HHT{..} = EMD (map go hhtLines) hhtResidual
where
go HHTLine{..} = SVG.zipWith (\m θ -> m * cos θ) hlMags θs
where
θs = SVG.scanl' (+) hlInitPhase ((* (2 * pi)) `SVG.map` hlFreqs)
ihht
:: (VG.Vector v a, Floating a)
=> HHT v n a
-> SVG.Vector v (n + 1) a
ihht = iemd . ihhtEmd
foldFreq
:: forall v u n a b c. (VG.Vector v a, VG.Vector u c, KnownNat n, Monoid b)
=> (a -> a -> b)
-> (b -> c)
-> HHT v n a
-> SVG.Vector u n c
foldFreq f g = pullBack
. foldl' (SV.zipWith (<>)) (SV.replicate mempty)
. map split
. hhtLines
where
split :: HHTLine v n a -> SV.Vector n b
split HHTLine{..} = SVG.generate $ \i ->
f (hlFreqs `SVG.index` i) (hlMags `SVG.index` weaken i)
{-# INLINE split #-}
pullBack :: SV.Vector n b -> SVG.Vector u n c
pullBack v = SVG.generate $ \i -> g (v `SV.index` i)
{-# INLINE pullBack #-}
{-# INLINE foldFreq #-}
newtype SumMap k a = SumMap { getSumMap :: M.Map k a }
instance (Ord k, Num a) => Semigroup (SumMap k a) where
SumMap x <> SumMap y = SumMap $ M.unionWith (+) x y
instance (Ord k, Num a) => Monoid (SumMap k a) where
mempty = SumMap M.empty
hhtSpectrum
:: forall v n a k. (VG.Vector v a, KnownNat n, Ord k, Num a)
=> (a -> k)
-> HHT v n a
-> SV.Vector n (M.Map k a)
hhtSpectrum f = foldFreq (\k -> SumMap . M.singleton (f k)) getSumMap
hhtSparseSpectrum
:: forall v n a k. (VG.Vector v a, KnownNat n, Ord k, Num a)
=> (a -> k)
-> HHT v n a
-> M.Map (Finite n, k) a
hhtSparseSpectrum f = M.unionsWith (+) . concatMap go . hhtLines
where
go :: HHTLine v n a -> [M.Map (Finite n, k) a]
go HHTLine{..} = flip fmap (finites @n) $ \i ->
M.singleton (i, f $ hlFreqs `SVG.index` i) $
hlMags `SVG.index` weaken i
hhtDenseSpectrum
:: forall v n m a. (VG.Vector v a, KnownNat n, KnownNat m, Num a)
=> (a -> Finite m)
-> HHT v n a
-> SV.Vector n (SV.Vector m a)
hhtDenseSpectrum f h = SV.generate $ \i -> SV.generate $ \j ->
M.findWithDefault 0 (i, j) ss
where
ss = hhtSparseSpectrum f h
marginal
:: forall v n a k. (VG.Vector v a, KnownNat n, Ord k, Num a)
=> (a -> k)
-> HHT v n a
-> M.Map k a
marginal f = M.unionsWith (+) . concatMap go . hhtLines
where
go :: HHTLine v n a -> [M.Map k a]
go HHTLine{..} = flip fmap (finites @n) $ \i ->
M.singleton (f $ hlFreqs `SVG.index` i) (hlMags `SVG.index` weaken i)
meanMarginal
:: forall v n a k. (VG.Vector v a, KnownNat n, Ord k, Fractional a)
=> (a -> k)
-> HHT v n a
-> M.Map k a
meanMarginal f = fmap (/ n) . marginal f
where
n = fromIntegral $ natVal (Proxy @n)
expectedFreq
:: forall v n a. (VG.Vector v a, KnownNat n, Fractional a)
=> HHT v n a
-> SVG.Vector v n a
expectedFreq = foldFreq (\x y -> (Sum (x * y), Sum y)) (\(Sum x, Sum y) -> x / y)
dominantFreq
:: forall v n a. (VG.Vector v a, KnownNat n, Ord a)
=> HHT v n a
-> SVG.Vector v n a
dominantFreq = foldFreq comb proj
where
comb :: a -> a -> Maybe (Max (Arg a a))
comb x y = Just $ Max $ Arg y x
proj :: Maybe (Max (Arg a a)) -> a
proj Nothing = errorWithoutStackTrace
"Numeric.HHT.dominantFreq: HHT was formed with no Intrinsic Mode Functions"
proj (Just (Max (Arg _ x))) = x
instantaneousEnergy
:: forall v n a. (VG.Vector v a, KnownNat n, Num a)
=> HHT v n a
-> SVG.Vector v n a
instantaneousEnergy = foldFreq (\_ x -> Sum (x * x)) getSum
degreeOfStationarity
:: forall v n a k. (VG.Vector v a, KnownNat n, Ord k, Fractional a, Eq a)
=> (a -> k)
-> HHT v n a
-> M.Map k a
degreeOfStationarity f h = fmap (/ n)
. M.unionsWith (+)
. concatMap go
. hhtLines
$ h
where
n = fromIntegral $ natVal (Proxy @n)
meanMarg = meanMarginal f h
go :: HHTLine v n a -> [M.Map k a]
go HHTLine{..} = flip fmap (finites @n) $ \i ->
let fr = f $ hlFreqs `SVG.index` i
mm = meanMarg M.! fr
in M.singleton fr $
if mm == 0
then 0
else (1 - (hlMags `SVG.index` weaken i / mm)) ^ (2 :: Int)
hilbertMagFreq
:: 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, a))
hilbertMagFreq v = (hilbertMag, (hilbertFreq, φ0))
where
v' = hilbert v
hilbertMag = SVG.map magnitude v'
hilbertPhase = SVG.map phase v'
φ0 = SVG.head hilbertPhase
hilbertFreq = SVG.map ((`mod'` 1) . (/ (2 * pi))) $ SVG.tail hilbertPhase - SVG.init hilbertPhase
hilbertPolar
:: 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)
where
hilbertMag :: SVG.Vector v (n + 1) a
hilbertFreq :: SVG.Vector v n a
(hilbertMag, (hilbertFreq, φ0)) = hilbertMagFreq v
hilbertPhase :: SVG.Vector v (n + 1) a
hilbertPhase = SVG.scanl' (+) φ0 ((* (2 * pi)) `SVG.map` hilbertFreq)
hilbert
:: 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 = 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)
hilbertIm
:: 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 = SVG.map imagPart . hilbert