-- SPDX-License-Identifier: MPL-2.0

{- | CPU reference DSP. Sample rates/frequencies are in Hz, times in seconds,
angular frequencies in radians/sample. Filtering uses causal zero-state
initial conditions and transposed direct form II; state is explicit and may
be carried across blocks. Convolution is full linear convolution; correlation
lags run -(length y-1)..length x-1 and real-valued y is reversed.

Filter prototypes use the classical bilinear transform and prewarped edges;
elliptic prototypes use Jacobi elliptic identities and the nome degree
equation (see <https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.ellipap.html>
and <https://dlmf.nist.gov/22.6>). Bessel is delay-normalized, so its supplied
edge is a prototype scale, not the Butterworth half-power frequency.
Wavelets use periodic extension; lengths must be divisible by 2^levels.
-}
module HMath.Signal (
    convolution,
    circularConvolution,
    correlation,
    autocorrelation,
    overlapAdd,
    overlapSave,
    Filter,
    iir,
    fir,
    filterCoefficients,
    filterState,
    resetFilter,
    filterSignal,
    FilterFamily (..),
    FilterBand (..),
    designIIR,
    frequencyResponse,
    notch,
    allPass,
    movingAverage,
    medianFilter,
    gaussianFilter,
    savitzkyGolay,
    wiener,
    Window (..),
    window,
    windowedSinc,
    downsample,
    upsample,
    resamplePolyphase,
    resample,
    antiAlias,
    Peak (..),
    peaks,
    zeroCrossings,
    envelope,
    unwrapPhase,
    detrend,
    normalize,
    smooth,
    noiseEstimate,
    signalToNoise,
    rms,
    GaussianState,
    gaussianState,
    stateMean,
    stateCovariance,
    kalman,
    extendedKalman,
    unscentedKalman,
    Wavelet (..),
    WaveletTree,
    dwt,
    idwt,
    waveletApproximation,
    waveletDetails,
    MotherWavelet (..),
    cwt,
    WaveletPacket,
    waveletPackets,
    inverseWaveletPackets,
    packetBands,
    waveletDenoise,
) where

import Control.Monad (foldM)
import Data.Complex (Complex ((:+)), cis, conjugate, imagPart, magnitude, realPart)
import Data.List (sort)
import HMath.Error (MathError (..))
import qualified HMath.Fourier as F
import qualified HMath.Interpolation as I
import qualified HMath.Polynomial as P
import HMath.Signal.Internal
import qualified HMath.Special as SF

-- | Direct full real linear convolution, with explicit quadratic work budget; an empty operand yields an empty result.
convolution :: [Double] -> [Double] -> Either MathError [Double]
convolution :: [Double] -> [Double] -> Either MathError [Double]
convolution [Double]
xs [Double]
ys = do
    Bool -> String -> Either MathError ()
require
        (Int -> Integer
forall a. Integral a => a -> Integer
toInteger ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
forall a. Integral a => a -> Integer
toInteger ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
ys) Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
10000000)
        String
"direct convolution work budget exceeded"
    _ <- [Double] -> Either MathError [Double]
vector ([Double]
xs [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double]
ys)
    if null xs || null ys
        then Right []
        else
            vector
                [ sum [xs !! j * ys !! (k - j) | j <- [max 0 (k - length ys + 1) .. min k (length xs - 1)]]
                | k <- [0 .. length xs + length ys - 2]
                ]

-- | Length-N periodic convolution of two equally sized nonempty real signals (N<=4096).
circularConvolution :: [Double] -> [Double] -> Either MathError [Double]
circularConvolution :: [Double] -> [Double] -> Either MathError [Double]
circularConvolution [Double]
x [Double]
y = do
    Bool -> String -> Either MathError ()
require
        (Bool -> Bool
not ([Double] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Double]
x) Bool -> Bool -> Bool
&& [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
x Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
y Bool -> Bool -> Bool
&& [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
x Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
4096)
        String
"circular convolution needs equal nonempty lengths <=4096"
    _ <- [Double] -> Either MathError [Double]
vector ([Double]
x [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double]
y)
    let n = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
x
    vector [sum [x !! j * y !! ((k - j) `mod` n) | j <- [0 .. n - 1]] | k <- [0 .. n - 1]]

-- | Full real cross-correlation, equivalent to convolution with the reversed second signal.
correlation :: [Double] -> [Double] -> Either MathError [Double]
correlation :: [Double] -> [Double] -> Either MathError [Double]
correlation [Double]
x [Double]
y = [Double] -> [Double] -> Either MathError [Double]
convolution [Double]
x ([Double] -> [Double]
forall a. [a] -> [a]
reverse [Double]
y)

-- | Unnormalized full autocorrelation with zero lag in the middle.
autocorrelation :: [Double] -> Either MathError [Double]
autocorrelation :: [Double] -> Either MathError [Double]
autocorrelation [Double]
x = [Double] -> [Double] -> Either MathError [Double]
correlation [Double]
x [Double]
x

-- | Block FFT linear convolution; block is the number of new input samples, not the padded transform length.
overlapAdd :: Int -> [Double] -> [Double] -> Either MathError [Double]
overlapAdd :: Int -> [Double] -> [Double] -> Either MathError [Double]
overlapAdd Int
block [Double]
kernel [Double]
xs = do
    Bool -> String -> Either MathError ()
require
        ( Int
block Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0
            Bool -> Bool -> Bool
&& Bool -> Bool
not ([Double] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Double]
kernel)
            Bool -> Bool -> Bool
&& Bool -> Bool
not ([Double] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Double]
xs)
            Bool -> Bool -> Bool
&& Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
block Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Int -> Integer
forall a. Integral a => a -> Integer
toInteger ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
kernel) Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
65537
        )
        String
"invalid overlap-add block/kernel/data"
    _ <- [Double] -> Either MathError [Double]
vector ([Double]
kernel [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double]
xs)
    pieces <-
        mapM
            (\[Double]
segment -> (Complex Double -> Double) -> [Complex Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Complex Double -> Double
forall a. Complex a -> a
realPart ([Complex Double] -> [Double])
-> Either MathError [Complex Double] -> Either MathError [Double]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> [Complex Double]
-> [Complex Double] -> Either MathError [Complex Double]
F.fftConvolve ((Double -> Complex Double) -> [Double] -> [Complex Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
0) [Double]
segment) ((Double -> Complex Double) -> [Double] -> [Complex Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
0) [Double]
kernel))
            (chunks block xs)
    let n = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs Int -> Int -> Int
forall a. Num a => a -> a -> a
+ [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
kernel Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1
    vector
        [ sum
            [row !! (k - start) | (start, row) <- zip [0, block ..] pieces, k >= start, k - start < length row]
        | k <- [0 .. n - 1]
        ]

-- | Streaming FFT convolution with explicit transform length >= kernel length; retained block size is N-M+1.
overlapSave :: Int -> [Double] -> [Double] -> Either MathError [Double]
overlapSave :: Int -> [Double] -> [Double] -> Either MathError [Double]
overlapSave Int
fftLength [Double]
kernel [Double]
xs = do
    let m :: Int
m = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
kernel; l :: Int
l = Int
fftLength Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1
    Bool -> String -> Either MathError ()
require
        (Int
m Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 Bool -> Bool -> Bool
&& Bool -> Bool
not ([Double] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Double]
xs) Bool -> Bool -> Bool
&& Int
fftLength Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
m Bool -> Bool -> Bool
&& Int
fftLength Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
65536)
        String
"invalid overlap-save transform length/kernel/data"
    _ <- [Double] -> Either MathError [Double]
vector ([Double]
kernel [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double]
xs)
    let padded = Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate (Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Double
0 [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double]
xs [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate (Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Double
0
        count = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1
        starts = (Int -> Bool) -> [Int] -> [Int]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
count) [Int
0, Int
l ..]
    h <- F.fft (map (:+ 0) (kernel ++ replicate (fftLength - m) 0))
    blocks <-
        mapM
            ( \Int
start -> do
                let segment :: [Double]
segment = Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
take Int
fftLength (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
start [Double]
padded [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ Double -> [Double]
forall a. a -> [a]
repeat Double
0)
                a <- [Complex Double] -> Either MathError [Complex Double]
F.fft ((Double -> Complex Double) -> [Double] -> [Complex Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
0) [Double]
segment)
                row <- F.ifft (zipWith (*) a h)
                Right (take l (drop (m - 1) (map realPart row)))
            )
            starts
    vector (take count (concat blocks))

-- | Opaque normalized numerator/denominator and caller-visible persistent delay state.
data Filter = Filter [Double] [Double] [Double] deriving (Filter -> Filter -> Bool
(Filter -> Filter -> Bool)
-> (Filter -> Filter -> Bool) -> Eq Filter
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Filter -> Filter -> Bool
== :: Filter -> Filter -> Bool
$c/= :: Filter -> Filter -> Bool
/= :: Filter -> Filter -> Bool
Eq, Int -> Filter -> ShowS
[Filter] -> ShowS
Filter -> String
(Int -> Filter -> ShowS)
-> (Filter -> String) -> ([Filter] -> ShowS) -> Show Filter
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Filter -> ShowS
showsPrec :: Int -> Filter -> ShowS
$cshow :: Filter -> String
show :: Filter -> String
$cshowList :: [Filter] -> ShowS
showList :: [Filter] -> ShowS
Show)

-- | Normalize finite ascending z^-1 coefficients by nonzero a0 and initialize a zero transposed-direct-form-II state.
iir :: [Double] -> [Double] -> Either MathError Filter
iir :: [Double] -> [Double] -> Either MathError Filter
iir [Double]
numerator [Double]
denominator = do
    Bool -> String -> Either MathError ()
require
        ( Bool -> Bool
not ([Double] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Double]
numerator)
            Bool -> Bool -> Bool
&& Bool -> Bool
not ([Double] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Double]
denominator)
            Bool -> Bool -> Bool
&& [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
numerator Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
257
            Bool -> Bool -> Bool
&& [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
denominator Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
257
            Bool -> Bool -> Bool
&& [Double]
denominator [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
0 Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
/= Double
0
        )
        String
"invalid IIR coefficients (1..257, nonzero a0)"
    _ <- [Double] -> Either MathError [Double]
vector ([Double]
numerator [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double]
denominator)
    let a0 = [Double]
denominator [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
0; order = Int -> Int -> Int
forall a. Ord a => a -> a -> a
max ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
numerator) ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
denominator) Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1
    b <- vector (map (/ a0) numerator ++ replicate (order + 1 - length numerator) 0)
    a <- vector (map (/ a0) denominator ++ replicate (order + 1 - length denominator) 0)
    Right (Filter b a (replicate order 0))

-- | Construct a causal zero-state finite-impulse-response filter; maximum 257 taps.
fir :: [Double] -> Either MathError Filter
fir :: [Double] -> Either MathError Filter
fir [Double]
b = [Double] -> [Double] -> Either MathError Filter
iir [Double]
b [Double
1]

-- | Read normalized padded (numerator,denominator) coefficients in increasing delay order.
filterCoefficients :: Filter -> ([Double], [Double])
filterCoefficients :: Filter -> ([Double], [Double])
filterCoefficients (Filter [Double]
b [Double]
a [Double]
_) = ([Double]
b, [Double]
a)

-- | Read explicit delay state after a block; the list is immutable and independently reusable.
filterState :: Filter -> [Double]
filterState :: Filter -> [Double]
filterState (Filter [Double]
_ [Double]
_ [Double]
s) = [Double]
s

-- | Preserve coefficients while replacing every delay with zero.
resetFilter :: Filter -> Filter
resetFilter :: Filter -> Filter
resetFilter (Filter [Double]
b [Double]
a [Double]
s) = [Double] -> [Double] -> [Double] -> Filter
Filter [Double]
b [Double]
a (Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
s) Double
0)

-- | Process a finite input block, returning both its output and the filter state for the next block.
filterSignal :: Filter -> [Double] -> Either MathError ([Double], Filter)
filterSignal :: Filter -> [Double] -> Either MathError ([Double], Filter)
filterSignal (Filter [Double]
b [Double]
a [Double]
state) [Double]
xs = do
    _ <- [Double] -> Either MathError [Double]
vector [Double]
xs
    let step ([Double]
acc, [Double]
s) Double
x = do
            let y :: Double
y = [Double]
b [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ if [Double] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Double]
s then Double
0 else [Double]
s [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
0
                next :: [Double]
next =
                    [ (if Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1 Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
s then [Double]
s [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) else Double
0) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ [Double]
b [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double]
a [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
y
                    | Int
i <- [Int
0 .. [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
s Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
                    ]
            _ <- [Double] -> Either MathError [Double]
vector (Double
y Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: [Double]
next)
            Right (y : acc, next)
    (ys, next) <- foldM step ([], state) xs
    Right (reverse ys, Filter b a next)

-- | Complex unit-circle response at omega radians/sample; reject poles and nonfinite arithmetic.
frequencyResponse :: Filter -> Double -> Either MathError (Complex Double)
frequencyResponse :: Filter -> Double -> Either MathError (Complex Double)
frequencyResponse (Filter [Double]
b [Double]
a [Double]
_) Double
omega = do
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
omega) String
"invalid angular frequency"
    let z :: Complex Double
z = Double -> Complex Double
forall a. Floating a => a -> Complex a
cis (-Double
omega)
        den :: Complex Double
den = [Complex Double] -> Complex Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Complex Double
coefficient Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
* Complex Double
z Complex Double -> Int -> Complex Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
j | (Int
j, Complex Double
coefficient) <- [Int] -> [Complex Double] -> [(Int, Complex Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 :: Int ..] ((Double -> Complex Double) -> [Double] -> [Complex Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
0) [Double]
a)]
        numerator :: Complex Double
numerator = [Complex Double] -> Complex Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Complex Double
coefficient Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
* Complex Double
z Complex Double -> Int -> Complex Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
j | (Int
j, Complex Double
coefficient) <- [Int] -> [Complex Double] -> [(Int, Complex Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 :: Int ..] ((Double -> Complex Double) -> [Double] -> [Complex Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
0) [Double]
b)]
    Bool -> String -> Either MathError ()
require (Complex Double -> Double
forall a. RealFloat a => Complex a -> a
magnitude Complex Double
den Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0) String
"filter pole at requested frequency"
    let value :: Complex Double
value = Complex Double
numerator Complex Double -> Complex Double -> Complex Double
forall a. Fractional a => a -> a -> a
/ Complex Double
den
    _ <- [Double] -> Either MathError [Double]
vector [Complex Double -> Double
forall a. Complex a -> a
realPart Complex Double
value, Complex Double -> Double
forall a. Complex a -> a
imagPart Complex Double
value]
    Right value

-- | Classical analog prototypes with ripple/attenuation parameters in decibels.
data FilterFamily
    = Butterworth
    | ChebyshevI Double
    | ChebyshevII Double
    | Elliptic Double Double
    | Bessel
    deriving (FilterFamily -> FilterFamily -> Bool
(FilterFamily -> FilterFamily -> Bool)
-> (FilterFamily -> FilterFamily -> Bool) -> Eq FilterFamily
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: FilterFamily -> FilterFamily -> Bool
== :: FilterFamily -> FilterFamily -> Bool
$c/= :: FilterFamily -> FilterFamily -> Bool
/= :: FilterFamily -> FilterFamily -> Bool
Eq, Int -> FilterFamily -> ShowS
[FilterFamily] -> ShowS
FilterFamily -> String
(Int -> FilterFamily -> ShowS)
-> (FilterFamily -> String)
-> ([FilterFamily] -> ShowS)
-> Show FilterFamily
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> FilterFamily -> ShowS
showsPrec :: Int -> FilterFamily -> ShowS
$cshow :: FilterFamily -> String
show :: FilterFamily -> String
$cshowList :: [FilterFamily] -> ShowS
showList :: [FilterFamily] -> ShowS
Show)

-- | Prewarped low/high edge or ordered band edges, all in Hz.
data FilterBand = Lowpass Double | Highpass Double | Bandpass Double Double | Bandstop Double Double
    deriving (FilterBand -> FilterBand -> Bool
(FilterBand -> FilterBand -> Bool)
-> (FilterBand -> FilterBand -> Bool) -> Eq FilterBand
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: FilterBand -> FilterBand -> Bool
== :: FilterBand -> FilterBand -> Bool
$c/= :: FilterBand -> FilterBand -> Bool
/= :: FilterBand -> FilterBand -> Bool
Eq, Int -> FilterBand -> ShowS
[FilterBand] -> ShowS
FilterBand -> String
(Int -> FilterBand -> ShowS)
-> (FilterBand -> String)
-> ([FilterBand] -> ShowS)
-> Show FilterBand
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> FilterBand -> ShowS
showsPrec :: Int -> FilterBand -> ShowS
$cshow :: FilterBand -> String
show :: FilterBand -> String
$cshowList :: [FilterBand] -> ShowS
showList :: [FilterBand] -> ShowS
Show)

{- | Orders 1..8; ripple/attenuation in dB; all edges strictly inside (0, fs/2).
Expanded coefficients are a bounded reference representation; high orders
should use cascaded sections in production applications.
-}
designIIR :: FilterFamily -> FilterBand -> Int -> Double -> Either MathError Filter
designIIR :: FilterFamily
-> FilterBand -> Int -> Double -> Either MathError Filter
designIIR FilterFamily
family FilterBand
band Int
order Double
rate = do
    Bool -> String -> Either MathError ()
require (Int
order Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
1 Bool -> Bool -> Bool
&& Int
order Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
8 Bool -> Bool -> Bool
&& Double -> Bool
finite Double
rate Bool -> Bool -> Bool
&& Double
rate Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0) String
"invalid IIR order/sample rate"
    let edges :: [Double]
edges = case FilterBand
band of Lowpass Double
f -> [Double
f]; Highpass Double
f -> [Double
f]; Bandpass Double
a Double
b -> [Double
a, Double
b]; Bandstop Double
a Double
b -> [Double
a, Double
b]
    Bool -> String -> Either MathError ()
require
        ((Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (\Double
f -> Double -> Bool
finite Double
f Bool -> Bool -> Bool
&& Double
f Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& Double
f Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
rate Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2) [Double]
edges Bool -> Bool -> Bool
&& [Bool] -> Bool
forall (t :: * -> *). Foldable t => t Bool -> Bool
and ((Double -> Double -> Bool) -> [Double] -> [Double] -> [Bool]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
(<) [Double]
edges (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
1 [Double]
edges)))
        String
"filter edges must increase inside Nyquist"
    (zeros, poles, gain) <- FilterFamily
-> Int
-> Either MathError ([Complex Double], [Complex Double], Double)
prototype FilterFamily
family Int
order
    let warp Double
f = Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
rate Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
tan (Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
f Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
rate)
        degree = [Complex Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Complex Double]
poles Int -> Int -> Int
forall a. Num a => a -> a -> a
- [Complex Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Complex Double]
zeros
        ratio = Complex Double -> Double
forall a. Complex a -> a
realPart ([Complex Double] -> Complex Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product ((Complex Double -> Complex Double)
-> [Complex Double] -> [Complex Double]
forall a b. (a -> b) -> [a] -> [b]
map Complex Double -> Complex Double
forall a. Num a => a -> a
negate [Complex Double]
zeros) Complex Double -> Complex Double -> Complex Double
forall a. Fractional a => a -> a -> a
/ [Complex Double] -> Complex Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product ((Complex Double -> Complex Double)
-> [Complex Double] -> [Complex Double]
forall a b. (a -> b) -> [a] -> [b]
map Complex Double -> Complex Double
forall a. Num a => a -> a
negate [Complex Double]
poles))
        quadratic Complex a
term a
center =
            let discriminant :: Complex a
discriminant = Complex a -> Complex a
forall a. Floating a => a -> a
sqrt (Complex a
term Complex a -> Complex a -> Complex a
forall a. Num a => a -> a -> a
* Complex a
term Complex a -> Complex a -> Complex a
forall a. Num a => a -> a -> a
- (a
4 a -> a -> a
forall a. Num a => a -> a -> a
* a
center a -> a -> a
forall a. Num a => a -> a -> a
* a
center a -> a -> Complex a
forall a. a -> a -> Complex a
:+ a
0))
             in [(Complex a
term Complex a -> Complex a -> Complex a
forall a. Num a => a -> a -> a
+ Complex a
discriminant) Complex a -> Complex a -> Complex a
forall a. Fractional a => a -> a -> a
/ Complex a
2, (Complex a
term Complex a -> Complex a -> Complex a
forall a. Num a => a -> a -> a
- Complex a
discriminant) Complex a -> Complex a -> Complex a
forall a. Fractional a => a -> a -> a
/ Complex a
2]
        (az, ap, ak) = case band of
            Lowpass Double
f -> let w :: Double
w = Double -> Double
warp Double
f in ((Complex Double -> Complex Double)
-> [Complex Double] -> [Complex Double]
forall a b. (a -> b) -> [a] -> [b]
map (Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
* (Double
w Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
0)) [Complex Double]
zeros, (Complex Double -> Complex Double)
-> [Complex Double] -> [Complex Double]
forall a b. (a -> b) -> [a] -> [b]
map (Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
* (Double
w Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
0)) [Complex Double]
poles, Double
gain Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
degree)
            Highpass Double
f ->
                let w :: Double
w = Double -> Double
warp Double
f
                 in ((Complex Double -> Complex Double)
-> [Complex Double] -> [Complex Double]
forall a b. (a -> b) -> [a] -> [b]
map ((Double
w Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
0) Complex Double -> Complex Double -> Complex Double
forall a. Fractional a => a -> a -> a
/) [Complex Double]
zeros [Complex Double] -> [Complex Double] -> [Complex Double]
forall a. [a] -> [a] -> [a]
++ Int -> Complex Double -> [Complex Double]
forall a. Int -> a -> [a]
replicate Int
degree Complex Double
0, (Complex Double -> Complex Double)
-> [Complex Double] -> [Complex Double]
forall a b. (a -> b) -> [a] -> [b]
map ((Double
w Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
0) Complex Double -> Complex Double -> Complex Double
forall a. Fractional a => a -> a -> a
/) [Complex Double]
poles, Double
gain Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
ratio)
            Bandpass Double
lo Double
hi ->
                let w1 :: Double
w1 = Double -> Double
warp Double
lo; w2 :: Double
w2 = Double -> Double
warp Double
hi; width :: Double
width = Double
w2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
w1; center :: Double
center = Double -> Double
forall a. Floating a => a -> a
sqrt (Double
w1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w2)
                 in ( (Complex Double -> [Complex Double])
-> [Complex Double] -> [Complex Double]
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap (\Complex Double
z -> Complex Double -> Double -> [Complex Double]
forall {a}. RealFloat a => Complex a -> a -> [Complex a]
quadratic (Complex Double
z Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
* (Double
width Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
0)) Double
center) [Complex Double]
zeros [Complex Double] -> [Complex Double] -> [Complex Double]
forall a. [a] -> [a] -> [a]
++ Int -> Complex Double -> [Complex Double]
forall a. Int -> a -> [a]
replicate Int
degree Complex Double
0
                    , (Complex Double -> [Complex Double])
-> [Complex Double] -> [Complex Double]
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap (\Complex Double
p -> Complex Double -> Double -> [Complex Double]
forall {a}. RealFloat a => Complex a -> a -> [Complex a]
quadratic (Complex Double
p Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
* (Double
width Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
0)) Double
center) [Complex Double]
poles
                    , Double
gain Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
width Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
degree
                    )
            Bandstop Double
lo Double
hi ->
                let w1 :: Double
w1 = Double -> Double
warp Double
lo; w2 :: Double
w2 = Double -> Double
warp Double
hi; width :: Double
width = Double
w2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
w1; center :: Double
center = Double -> Double
forall a. Floating a => a -> a
sqrt (Double
w1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w2)
                 in ( (Complex Double -> [Complex Double])
-> [Complex Double] -> [Complex Double]
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap (\Complex Double
z -> Complex Double -> Double -> [Complex Double]
forall {a}. RealFloat a => Complex a -> a -> [Complex a]
quadratic ((Double
width Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
0) Complex Double -> Complex Double -> Complex Double
forall a. Fractional a => a -> a -> a
/ Complex Double
z) Double
center) [Complex Double]
zeros
                        [Complex Double] -> [Complex Double] -> [Complex Double]
forall a. [a] -> [a] -> [a]
++ [[Complex Double]] -> [Complex Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat (Int -> [Complex Double] -> [[Complex Double]]
forall a. Int -> a -> [a]
replicate Int
degree [Double
0 Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
center, Double
0 Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ (-Double
center)])
                    , (Complex Double -> [Complex Double])
-> [Complex Double] -> [Complex Double]
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap (\Complex Double
p -> Complex Double -> Double -> [Complex Double]
forall {a}. RealFloat a => Complex a -> a -> [Complex a]
quadratic ((Double
width Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
0) Complex Double -> Complex Double -> Complex Double
forall a. Fractional a => a -> a -> a
/ Complex Double
p) Double
center) [Complex Double]
poles
                    , Double
gain Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
ratio
                    )
        c = Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
rate Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
0
        dz = (Complex Double -> Complex Double)
-> [Complex Double] -> [Complex Double]
forall a b. (a -> b) -> [a] -> [b]
map (\Complex Double
z -> (Complex Double
c Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
+ Complex Double
z) Complex Double -> Complex Double -> Complex Double
forall a. Fractional a => a -> a -> a
/ (Complex Double
c Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
- Complex Double
z)) [Complex Double]
az [Complex Double] -> [Complex Double] -> [Complex Double]
forall a. [a] -> [a] -> [a]
++ Int -> Complex Double -> [Complex Double]
forall a. Int -> a -> [a]
replicate ([Complex Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Complex Double]
ap Int -> Int -> Int
forall a. Num a => a -> a -> a
- [Complex Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Complex Double]
az) (-Complex Double
1)
        dp = (Complex Double -> Complex Double)
-> [Complex Double] -> [Complex Double]
forall a b. (a -> b) -> [a] -> [b]
map (\Complex Double
p -> (Complex Double
c Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
+ Complex Double
p) Complex Double -> Complex Double -> Complex Double
forall a. Fractional a => a -> a -> a
/ (Complex Double
c Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
- Complex Double
p)) [Complex Double]
ap
        dk = Double
ak Double -> Double -> Double
forall a. Num a => a -> a -> a
* Complex Double -> Double
forall a. Complex a -> a
realPart ([Complex Double] -> Complex Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product ((Complex Double -> Complex Double)
-> [Complex Double] -> [Complex Double]
forall a b. (a -> b) -> [a] -> [b]
map (Complex Double
c Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
-) [Complex Double]
az) Complex Double -> Complex Double -> Complex Double
forall a. Fractional a => a -> a -> a
/ [Complex Double] -> Complex Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product ((Complex Double -> Complex Double)
-> [Complex Double] -> [Complex Double]
forall a b. (a -> b) -> [a] -> [b]
map (Complex Double
c Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
-) [Complex Double]
ap))
        coeffs = ([Complex Double] -> Complex Double -> [Complex Double])
-> [Complex Double] -> [Complex Double] -> [Complex Double]
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl (\[Complex Double]
cs Complex Double
z -> [Complex Double] -> [Complex Double] -> [Complex Double]
forall {a}. Num a => [a] -> [a] -> [a]
zipPad [Complex Double]
cs (Complex Double
0 Complex Double -> [Complex Double] -> [Complex Double]
forall a. a -> [a] -> [a]
: (Complex Double -> Complex Double)
-> [Complex Double] -> [Complex Double]
forall a b. (a -> b) -> [a] -> [b]
map (Complex Double -> Complex Double
forall a. Num a => a -> a
negate Complex Double
z Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
*) [Complex Double]
cs)) [Complex Double
1]
        zipPad [] [a]
ys = [a]
ys; zipPad [a]
xs [] = [a]
xs; zipPad (a
x : [a]
xs) (a
y : [a]
ys) = (a
x a -> a -> a
forall a. Num a => a -> a -> a
+ a
y) a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a] -> [a] -> [a]
zipPad [a]
xs [a]
ys
        numerator = (Complex Double -> Complex Double)
-> [Complex Double] -> [Complex Double]
forall a b. (a -> b) -> [a] -> [b]
map (Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
* (Double
dk Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
0)) ([Complex Double] -> [Complex Double]
coeffs [Complex Double]
dz)
        denominator = [Complex Double] -> [Complex Double]
coeffs [Complex Double]
dp
    require (all (\Complex Double
p -> Complex Double -> Double
forall a. RealFloat a => Complex a -> a
magnitude Complex Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
1) dp) "designed digital filter is unstable"
    require
        (all (\Complex Double
z -> Double -> Double
forall a. Num a => a -> a
abs (Complex Double -> Double
forall a. Complex a -> a
imagPart Complex Double
z) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
1e-8 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
1 (Double -> Double
forall a. Num a => a -> a
abs (Complex Double -> Double
forall a. Complex a -> a
realPart Complex Double
z))) (numerator ++ denominator))
        "filter coefficients lost conjugate symmetry"
    iir (map realPart numerator) (map realPart denominator)

prototype :: FilterFamily -> Int -> Either MathError ([Complex Double], [Complex Double], Double)
prototype :: FilterFamily
-> Int
-> Either MathError ([Complex Double], [Complex Double], Double)
prototype FilterFamily
family Int
n = case FilterFamily
family of
    FilterFamily
Butterworth ->
        let poles :: [Complex Double]
poles = [Double -> Complex Double
forall a. Floating a => a -> Complex a
cis (Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n)) | Int
k <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
         in ([Complex Double], [Complex Double], Double)
-> Either MathError ([Complex Double], [Complex Double], Double)
forall a b. b -> Either a b
Right ([], [Complex Double]
poles, Complex Double -> Double
forall a. Complex a -> a
realPart ([Complex Double] -> Complex Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product ((Complex Double -> Complex Double)
-> [Complex Double] -> [Complex Double]
forall a b. (a -> b) -> [a] -> [b]
map Complex Double -> Complex Double
forall a. Num a => a -> a
negate [Complex Double]
poles)))
    ChebyshevI Double
ripple -> do
        Bool -> String -> Either MathError ()
require
            (Double -> Bool
finite Double
ripple Bool -> Bool -> Bool
&& Double
ripple Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& Double
ripple Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
100)
            String
"Chebyshev I ripple must be positive <=100 dB"
        let eps :: Double
eps = Double -> Double
forall a. Floating a => a -> a
sqrt (Double
10 Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (Double
ripple Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
10) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1)
            mu :: Double
mu = Double -> Double
forall a. Floating a => a -> a
asinh (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
eps) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
            poles :: [Complex Double]
poles =
                [ let theta :: Double
theta = Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n)
                   in (-Double -> Double
forall a. Floating a => a -> a
sinh Double
mu Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin Double
theta) Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ (Double -> Double
forall a. Floating a => a -> a
cosh Double
mu Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos Double
theta)
                | Int
k <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
                ]
            gain :: Double
gain = Complex Double -> Double
forall a. Complex a -> a
realPart ([Complex Double] -> Complex Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product ((Complex Double -> Complex Double)
-> [Complex Double] -> [Complex Double]
forall a b. (a -> b) -> [a] -> [b]
map Complex Double -> Complex Double
forall a. Num a => a -> a
negate [Complex Double]
poles)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (if Int -> Bool
forall a. Integral a => a -> Bool
even Int
n then Double -> Double
forall a. Floating a => a -> a
sqrt (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
eps Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
eps) else Double
1)
        ([Complex Double], [Complex Double], Double)
-> Either MathError ([Complex Double], [Complex Double], Double)
forall a b. b -> Either a b
Right ([], [Complex Double]
poles, Double
gain)
    ChebyshevII Double
attenuation -> do
        Bool -> String -> Either MathError ()
require
            (Double -> Bool
finite Double
attenuation Bool -> Bool -> Bool
&& Double
attenuation Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& Double
attenuation Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
300)
            String
"Chebyshev II attenuation must be positive <=300 dB"
        let mu :: Double
mu = Double -> Double
forall a. Floating a => a -> a
asinh (Double -> Double
forall a. Floating a => a -> a
sqrt (Double
10 Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (Double
attenuation Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
10) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
            poles :: [Complex Double]
poles =
                [ Complex Double -> Complex Double
forall a. Fractional a => a -> a
recip ((-Double -> Double
forall a. Floating a => a -> a
sinh Double
mu Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin Double
theta) Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ (Double -> Double
forall a. Floating a => a -> a
cosh Double
mu Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos Double
theta))
                | Int
k <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
                , let theta :: Double
theta = Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n)
                ]
            zeros :: [Complex Double]
zeros =
                [ Double
0 Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
cos Double
theta)
                | Int
k <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
                , Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
n
                , let theta :: Double
theta = Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n)
                ]
        ([Complex Double], [Complex Double], Double)
-> Either MathError ([Complex Double], [Complex Double], Double)
forall a b. b -> Either a b
Right ([Complex Double]
zeros, [Complex Double]
poles, Complex Double -> Double
forall a. Complex a -> a
realPart ([Complex Double] -> Complex Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product ((Complex Double -> Complex Double)
-> [Complex Double] -> [Complex Double]
forall a b. (a -> b) -> [a] -> [b]
map Complex Double -> Complex Double
forall a. Num a => a -> a
negate [Complex Double]
poles) Complex Double -> Complex Double -> Complex Double
forall a. Fractional a => a -> a -> a
/ [Complex Double] -> Complex Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product ((Complex Double -> Complex Double)
-> [Complex Double] -> [Complex Double]
forall a b. (a -> b) -> [a] -> [b]
map Complex Double -> Complex Double
forall a. Num a => a -> a
negate [Complex Double]
zeros)))
    Elliptic Double
ripple Double
attenuation -> do
        Bool -> String -> Either MathError ()
require
            (Double -> Bool
finite Double
ripple Bool -> Bool -> Bool
&& Double -> Bool
finite Double
attenuation Bool -> Bool -> Bool
&& Double
ripple Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& Double
attenuation Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
ripple Bool -> Bool -> Bool
&& Double
attenuation Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
150)
            String
"elliptic requires 0<ripple<attenuation<=150 dB"
        let eps :: Double
eps = Double -> Double
forall a. Floating a => a -> a
sqrt (Double
10 Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (Double
ripple Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
10) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1); m1 :: Double
m1 = Double
eps Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
eps Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
10 Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (Double
attenuation Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
10) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1)
        Bool -> String -> Either MathError ()
require (Double
m1 Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
1e-14 Bool -> Bool -> Bool
&& Double
m1 Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
1) String
"elliptic selectivity outside reference precision"
        k1 <- Double -> Either MathError Double
completeK Double
m1
        k1p <- completeK (1 - m1)
        let q = Double -> Double
forall a. Floating a => a -> a
exp (-Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
k1p Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
k1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n)
            theta2 = Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
q Double -> Double -> Double
forall a. Floating a => a -> a -> a
** Double
0.25 Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double
q Double -> Double -> Double
forall a. Floating a => a -> a -> a
** Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
* (Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)) | Int
j <- [Int
0 .. Int
100 :: Int]]
            theta3 = Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double
q Double -> Double -> Double
forall a. Floating a => a -> a -> a
** Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
j) | Int
j <- [Int
1 .. Int
100 :: Int]]
            m = (Double
theta2 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
theta3) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
4 :: Int)
        require (m > 0 && m < 1 - 1e-12) "elliptic degree equation ill-conditioned"
        kk <- completeK m
        inverseSC <- SF.sfValue <$> SF.ellipticF (atan (1 / eps)) (1 - m1)
        let v0 = Double
kk Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
inverseSC Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
k1)
        (sv, cv, dv) <- SF.jacobiElliptic v0 (1 - m)
        triples <-
            mapM
                (\Int
j -> Double -> Double -> Either MathError (Double, Double, Double)
SF.jacobiElliptic (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
j Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
kk Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n) Double
m)
                [1 - n `mod` 2, 3 - n `mod` 2 .. n - 1]
        let positiveZeros = [Double
0 Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double -> Double
forall a. Floating a => a -> a
sqrt Double
m Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
s)) | (Double
s, Double
_, Double
_) <- [(Double, Double, Double)]
triples, Double -> Double
forall a. Num a => a -> a
abs Double
s Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
1e-14]
            zeros = [Complex Double]
positiveZeros [Complex Double] -> [Complex Double] -> [Complex Double]
forall a. [a] -> [a] -> [a]
++ (Complex Double -> Complex Double)
-> [Complex Double] -> [Complex Double]
forall a b. (a -> b) -> [a] -> [b]
map Complex Double -> Complex Double
forall a. Num a => Complex a -> Complex a
conjugate [Complex Double]
positiveZeros
            halfPoles =
                [Complex Double -> Complex Double
forall a. Num a => a -> a
negate ((Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
d Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sv Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cv) Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ (Double
s Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
dv)) Complex Double -> Complex Double -> Complex Double
forall a. Fractional a => a -> a -> a
/ ((Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- (Double
d Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sv) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int)) Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
0) | (Double
s, Double
c, Double
d) <- [(Double, Double, Double)]
triples]
            poles = [Complex Double]
halfPoles [Complex Double] -> [Complex Double] -> [Complex Double]
forall a. [a] -> [a] -> [a]
++ (Complex Double -> Complex Double)
-> [Complex Double] -> [Complex Double]
forall a b. (a -> b) -> [a] -> [b]
map Complex Double -> Complex Double
forall a. Num a => Complex a -> Complex a
conjugate [Complex Double
p | Complex Double
p <- [Complex Double]
halfPoles, Double -> Double
forall a. Num a => a -> a
abs (Complex Double -> Double
forall a. Complex a -> a
imagPart Complex Double
p) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
1e-12]
            gain =
                Complex Double -> Double
forall a. Complex a -> a
realPart ([Complex Double] -> Complex Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product ((Complex Double -> Complex Double)
-> [Complex Double] -> [Complex Double]
forall a b. (a -> b) -> [a] -> [b]
map Complex Double -> Complex Double
forall a. Num a => a -> a
negate [Complex Double]
poles) Complex Double -> Complex Double -> Complex Double
forall a. Fractional a => a -> a -> a
/ [Complex Double] -> Complex Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product ((Complex Double -> Complex Double)
-> [Complex Double] -> [Complex Double]
forall a b. (a -> b) -> [a] -> [b]
map Complex Double -> Complex Double
forall a. Num a => a -> a
negate [Complex Double]
zeros))
                    Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (if Int -> Bool
forall a. Integral a => a -> Bool
even Int
n then Double -> Double
forall a. Floating a => a -> a
sqrt (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
eps Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
eps) else Double
1)
        require (length poles == n) "elliptic pole count mismatch"
        Right (zeros, poles, gain)
    FilterFamily
Bessel -> do
        let factorial :: a -> Double
factorial a
k = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Double
1 .. a -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
k] :: Double
            cs :: [Complex Double]
cs =
                [ Int -> Double
forall {a}. Integral a => a -> Double
factorial (Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
k) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
2 Double -> Double -> Double
forall a. Floating a => a -> a -> a
** Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
k) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall {a}. Integral a => a -> Double
factorial Int
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall {a}. Integral a => a -> Double
factorial (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
k)) Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
0
                | Int
k <- [Int
0 .. Int
n]
                ]
        poles <- Int
-> Double
-> Polynomial (Complex Double)
-> Either MathError [Complex Double]
P.aberthEhrlich Int
2000 Double
1e-11 ([Complex Double] -> Polynomial (Complex Double)
forall a. (Eq a, Num a) => [a] -> Polynomial a
P.polynomial [Complex Double]
cs)
        require (all ((< 0) . realPart) poles) "Bessel root solver returned nonstable pole"
        Right ([], poles, realPart (product (map negate poles)))
  where
    completeK :: Double -> Either MathError Double
completeK Double
m = do
        Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
m Bool -> Bool -> Bool
&& Double
m Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0 Bool -> Bool -> Bool
&& Double
m Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
1) String
"elliptic K parameter"
        let loop :: t -> t -> t -> Either MathError t
loop t
0 t
_ t
_ = MathError -> Either MathError t
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"elliptic AGM" Int
100)
            loop t
k t
a t
b =
                if t -> t
forall a. Num a => a -> a
abs (t
a t -> t -> t
forall a. Num a => a -> a -> a
- t
b) t -> t -> Bool
forall a. Ord a => a -> a -> Bool
<= t
1e-15 t -> t -> t
forall a. Num a => a -> a -> a
* t
a then t -> Either MathError t
forall a b. b -> Either a b
Right (t
forall a. Floating a => a
pi t -> t -> t
forall a. Fractional a => a -> a -> a
/ (t
2 t -> t -> t
forall a. Num a => a -> a -> a
* t
a)) else t -> t -> t -> Either MathError t
loop (t
k t -> t -> t
forall a. Num a => a -> a -> a
- t
1) ((t
a t -> t -> t
forall a. Num a => a -> a -> a
+ t
b) t -> t -> t
forall a. Fractional a => a -> a -> a
/ t
2) (t -> t
forall a. Floating a => a -> a
sqrt (t
a t -> t -> t
forall a. Num a => a -> a -> a
* t
b))
        Int -> Double -> Double -> Either MathError Double
forall {t} {t}.
(Ord t, Num t, Floating t, Eq t) =>
t -> t -> t -> Either MathError t
loop (Int
100 :: Int) Double
1 (Double -> Double
forall a. Floating a => a -> a
sqrt (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
m))

-- | Second-order notch or unity-magnitude all-pass with center omega in (0,pi) and positive quality factor Q.
notch, allPass :: Double -> Double -> Either MathError Filter
notch :: Double -> Double -> Either MathError Filter
notch Double
omega Double
quality = do
    Bool -> String -> Either MathError ()
require
        (Double -> Bool
finite Double
omega Bool -> Bool -> Bool
&& Double
omega Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& Double
omega Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
forall a. Floating a => a
pi Bool -> Bool -> Bool
&& Double -> Bool
finite Double
quality Bool -> Bool -> Bool
&& Double
quality Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0)
        String
"invalid notch frequency/Q"
    let alpha :: Double
alpha = Double -> Double
forall a. Floating a => a -> a
sin Double
omega Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
quality); c :: Double
c = Double -> Double
forall a. Floating a => a -> a
cos Double
omega
    [Double] -> [Double] -> Either MathError Filter
iir [Double
1, -Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
c, Double
1] [Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
alpha, -Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
c, Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
alpha]
allPass :: Double -> Double -> Either MathError Filter
allPass Double
omega Double
quality = do
    Bool -> String -> Either MathError ()
require
        (Double -> Bool
finite Double
omega Bool -> Bool -> Bool
&& Double
omega Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& Double
omega Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
forall a. Floating a => a
pi Bool -> Bool -> Bool
&& Double -> Bool
finite Double
quality Bool -> Bool -> Bool
&& Double
quality Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0)
        String
"invalid all-pass frequency/Q"
    let alpha :: Double
alpha = Double -> Double
forall a. Floating a => a -> a
sin Double
omega Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
quality); c :: Double
c = Double -> Double
forall a. Floating a => a -> a
cos Double
omega
    [Double] -> [Double] -> Either MathError Filter
iir [Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
alpha, -Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
c, Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
alpha] [Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
alpha, -Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
c, Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
alpha]

-- | Causal zero-state box FIR with supplied width; initial outputs include implicit leading zeros.
movingAverage :: Int -> [Double] -> Either MathError [Double]
movingAverage :: Int -> [Double] -> Either MathError [Double]
movingAverage Int
count [Double]
xs = do
    Bool -> String -> Either MathError ()
require (Int
count Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 Bool -> Bool -> Bool
&& Int
count Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
257) String
"moving average width must be 1..257"
    f <- [Double] -> Either MathError Filter
fir (Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate Int
count (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
count))
    fst <$> filterSignal f xs

{- | Centered finite filters clamp sample indices at endpoints (replicate
padding), preserve output length, and require odd positive support widths.
-}
localFilter :: Int -> ([Double] -> Either MathError Double) -> [Double] -> Either MathError [Double]
localFilter :: Int
-> ([Double] -> Either MathError Double)
-> [Double]
-> Either MathError [Double]
localFilter Int
width [Double] -> Either MathError Double
f [Double]
xs = do
    Bool -> String -> Either MathError ()
require
        (Int
width Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 Bool -> Bool -> Bool
&& Int -> Bool
forall a. Integral a => a -> Bool
odd Int
width Bool -> Bool -> Bool
&& Int
width Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
1025 Bool -> Bool -> Bool
&& Bool -> Bool
not ([Double] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Double]
xs))
        String
"invalid centered filter width/data"
    _ <- [Double] -> Either MathError [Double]
vector [Double]
xs
    let radius = Int
width Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2
    mapM
        (\Int
i -> [Double] -> Either MathError Double
f [[Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
0 (Int -> Int -> Int
forall a. Ord a => a -> a -> a
min ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
j)) | Int
j <- [-Int
radius .. Int
radius]])
        [0 .. length xs - 1]

-- | Centered odd-width running median using replicate endpoint padding.
medianFilter :: Int -> [Double] -> Either MathError [Double]
medianFilter :: Int -> [Double] -> Either MathError [Double]
medianFilter Int
width = Int
-> ([Double] -> Either MathError Double)
-> [Double]
-> Either MathError [Double]
localFilter Int
width (\[Double]
ys -> Double -> Either MathError Double
forall a b. b -> Either a b
Right ([Double] -> [Double]
forall a. Ord a => [a] -> [a]
sort [Double]
ys [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
ys Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2)))

-- | Centered normalized Gaussian convolution with supplied sigma in samples and odd support width.
gaussianFilter :: Int -> Double -> [Double] -> Either MathError [Double]
gaussianFilter :: Int -> Double -> [Double] -> Either MathError [Double]
gaussianFilter Int
width Double
sigma [Double]
xs = do
    Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
sigma Bool -> Bool -> Bool
&& Double
sigma Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0) String
"Gaussian sigma must be positive"
    let r :: Int
r = Int
width Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2; weights :: [Double]
weights = [Double -> Double
forall a. Floating a => a -> a
exp (-Double
0.5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
j Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
sigma) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int)) | Int
j <- [-Int
r .. Int
r]]
    Int
-> ([Double] -> Either MathError Double)
-> [Double]
-> Either MathError [Double]
localFilter Int
width (Double -> Either MathError Double
checked (Double -> Either MathError Double)
-> ([Double] -> Double) -> [Double] -> Either MathError Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double]
weights) (Double -> Double) -> ([Double] -> Double) -> [Double] -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Double] -> [Double] -> Double
dot [Double]
weights) [Double]
xs

-- | Centered local QR polynomial fit; choose width, degree, derivative order and positive sample spacing.
savitzkyGolay :: Int -> Int -> Int -> Double -> [Double] -> Either MathError [Double]
savitzkyGolay :: Int
-> Int -> Int -> Double -> [Double] -> Either MathError [Double]
savitzkyGolay Int
width Int
degree Int
derivative Double
spacing [Double]
xs = do
    Bool -> String -> Either MathError ()
require
        ( Int
degree Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0
            Bool -> Bool -> Bool
&& Int
degree Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
width
            Bool -> Bool -> Bool
&& Int
degree Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
16
            Bool -> Bool -> Bool
&& Int
derivative Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0
            Bool -> Bool -> Bool
&& Int
derivative Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
degree
            Bool -> Bool -> Bool
&& Double -> Bool
finite Double
spacing
            Bool -> Bool -> Bool
&& Double
spacing Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0
        )
        String
"invalid Savitzky-Golay polynomial/derivative/spacing"
    let r :: Int
r = Int
width Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2
        coordinates :: [Double]
coordinates = [Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
j Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
spacing | Int
j <- [-Int
r .. Int
r]]
        factorial :: Double
factorial = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Double
1 .. Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
derivative]
    Int
-> ([Double] -> Either MathError Double)
-> [Double]
-> Either MathError [Double]
localFilter
        Int
width
        ( \[Double]
ys -> do
            fit <- Int -> [(Double, Double)] -> Either MathError Approximation
I.polynomialFit Int
degree ([Double] -> [Double] -> [(Double, Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Double]
coordinates [Double]
ys)
            checked (factorial * (I.approximationCoefficients fit !! derivative))
        )
        [Double]
xs

-- | Local mean/variance Wiener shrinkage; optional noise variance, otherwise average local variance is estimated.
wiener :: Int -> Maybe Double -> [Double] -> Either MathError [Double]
wiener :: Int -> Maybe Double -> [Double] -> Either MathError [Double]
wiener Int
width Maybe Double
noise [Double]
xs = do
    means <- Int
-> ([Double] -> Either MathError Double)
-> [Double]
-> Either MathError [Double]
localFilter Int
width (\[Double]
ys -> Double -> Either MathError Double
checked ([Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double]
ys Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
ys))) [Double]
xs
    variances <-
        localFilter
            width
            ( \[Double]
ys ->
                let m :: Double
m = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double]
ys Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
ys)
                 in Double -> Either MathError Double
checked ([Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [(Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
m) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int) | Double
y <- [Double]
ys] Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
ys))
            )
            xs
    let nv = Double -> (Double -> Double) -> Maybe Double -> Double
forall b a. b -> (a -> b) -> Maybe a -> b
maybe ([Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double]
variances Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
variances)) Double -> Double
forall a. a -> a
id Maybe Double
noise
    require (finite nv && nv >= 0) "Wiener noise variance must be nonnegative"
    vector (zipWith3 (\Double
x Double
m Double
v -> Double
m Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
m) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (if Double
v Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
nv then Double
0 else Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
nv Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
v)) xs means variances)

-- | Symmetric endpoint-inclusive windows, with shape parameters beta/sigma/alpha/attenuation as appropriate.
data Window
    = Hann
    | Hamming
    | Blackman
    | Bartlett
    | Kaiser Double
    | Gaussian Double
    | Tukey Double
    | FlatTop
    | DolphChebyshev Double
    deriving (Window -> Window -> Bool
(Window -> Window -> Bool)
-> (Window -> Window -> Bool) -> Eq Window
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Window -> Window -> Bool
== :: Window -> Window -> Bool
$c/= :: Window -> Window -> Bool
/= :: Window -> Window -> Bool
Eq, Int -> Window -> ShowS
[Window] -> ShowS
Window -> String
(Int -> Window -> ShowS)
-> (Window -> String) -> ([Window] -> ShowS) -> Show Window
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Window -> ShowS
showsPrec :: Int -> Window -> ShowS
$cshow :: Window -> String
show :: Window -> String
$cshowList :: [Window] -> ShowS
showList :: [Window] -> ShowS
Show)

-- | Generate a finite symmetric window. Gaussian sigma is relative to half length; Dolph attenuation is in dB.
window :: Window -> Int -> Either MathError [Double]
window :: Window -> Int -> Either MathError [Double]
window Window
kind Int
n = do
    Bool -> String -> Either MathError ()
require (Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 Bool -> Bool -> Bool
&& Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
65536) String
"window length must be 1..65536"
    case Window
kind of
        Kaiser Double
b -> Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
b Bool -> Bool -> Bool
&& Double
b Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0 Bool -> Bool -> Bool
&& Double
b Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
100) String
"Kaiser beta 0..100"
        Gaussian Double
s -> Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
s Bool -> Bool -> Bool
&& Double
s Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0) String
"Gaussian window sigma positive"
        Tukey Double
a -> Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
a Bool -> Bool -> Bool
&& Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0 Bool -> Bool -> Bool
&& Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
1) String
"Tukey alpha 0..1"
        DolphChebyshev Double
db -> Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
db Bool -> Bool -> Bool
&& Double
db Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& Double
db Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
300 Bool -> Bool -> Bool
&& Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
4096) String
"Dolph-Chebyshev attenuation/length"
        Window
_ -> () -> Either MathError ()
forall a b. b -> Either a b
Right ()
    if Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1
        then [Double] -> Either MathError [Double]
forall a b. b -> Either a b
Right [Double
1]
        else do
            let m :: Double
m = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)
                theta :: a -> Double
theta a
j = Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* a -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
j Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
m
                i0 :: a -> a
i0 a
z = [a] -> a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take Int
100 ((a -> Int -> a) -> a -> [Int] -> [a]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl (\a
t Int
k -> a
t a -> a -> a
forall a. Num a => a -> a -> a
* a
z a -> a -> a
forall a. Num a => a -> a -> a
* a
z a -> a -> a
forall a. Fractional a => a -> a -> a
/ (a
4 a -> a -> a
forall a. Num a => a -> a -> a
* Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
k))) a
1 [Int
1 .. Int
99 :: Int]))
                cheb :: a -> a -> a
cheb a
degree a
x
                    | a -> a
forall a. Num a => a -> a
abs a
x a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
1 = a -> a
forall a. Floating a => a -> a
cos (a -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
degree a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Floating a => a -> a
acos a
x)
                    | a
x a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> a
1 = a -> a
forall a. Floating a => a -> a
cosh (a -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
degree a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Floating a => a -> a
acosh a
x)
                    | Bool
otherwise = (if a -> Bool
forall a. Integral a => a -> Bool
even a
degree then a
1 else -a
1) a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Floating a => a -> a
cosh (a -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
degree a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Floating a => a -> a
acosh (-a
x))
                coefficient :: a -> Double
coefficient a
j =
                    let x :: Double
x = a -> Double
forall {a}. Integral a => a -> Double
theta a
j; t :: Double
t = a -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
j Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
m; center :: Double
center = (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* a -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
j Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
m) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
m
                     in case Window
kind of
                            Window
Hann -> Double
0.5 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
0.5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos Double
x
                            Window
Hamming -> Double
0.54 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
0.46 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos Double
x
                            Window
Blackman -> Double
0.42 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
0.5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
0.08 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x)
                            Window
Bartlett -> Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double
forall a. Num a => a -> a
abs Double
center
                            Kaiser Double
b -> Double -> Double
forall a. Fractional a => a -> a
i0 (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sqrt (Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
0 (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
center Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
center))) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Fractional a => a -> a
i0 Double
b
                            Gaussian Double
sigma -> Double -> Double
forall a. Floating a => a -> a
exp (-Double
0.5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
center Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
sigma) Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int))
                            Tukey Double
alpha
                                | Double
alpha Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 -> Double
1
                                | Double
t Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
alpha Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2 -> Double
0.5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double
forall a. Floating a => a -> a
cos (Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
t Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
alpha Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1)))
                                | Double
t Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
alpha Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2 -> Double
0.5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double
forall a. Floating a => a -> a
cos (Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
t Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
alpha Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
alpha Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1)))
                                | Bool
otherwise -> Double
1
                            Window
FlatTop ->
                                Double
0.21557895
                                    Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
0.41663158 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos Double
x
                                    Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
0.277263158 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x)
                                    Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
0.083578947 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos (Double
3 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x)
                                    Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
0.006947368 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos (Double
4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x)
                            DolphChebyshev Double
db ->
                                let beta :: Double
beta = Double -> Double
forall a. Floating a => a -> a
cosh (Double -> Double
forall a. Floating a => a -> a
acosh (Double
10 Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (Double
db Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
20)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
m)
                                 in [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum
                                        [ Int -> Double -> Double
forall {a} {a}. (Floating a, Integral a, Ord a) => a -> a -> a
cheb (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) (Double
beta Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos (Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n))
                                            Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* (a -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
j Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
m Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n)
                                        | Int
k <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
                                        ]
                                        Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
                raw :: [Double]
raw = (Int -> Double) -> [Int] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Int -> Double
forall {a}. Integral a => a -> Double
coefficient [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
                result :: [Double]
result = case Window
kind of DolphChebyshev Double
_ -> let scale :: Double
scale = [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Double -> Double
forall a. Num a => a -> a
abs [Double]
raw) in (Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
scale) [Double]
raw; Window
_ -> [Double]
raw
            [Double] -> Either MathError [Double]
vector [Double]
result

{- | Given window kind, odd tap count in [1,257], cutoff frequency and positive sample rate,
return normalized low-pass FIR coefficients. Require 0\<cutoff\<rate/2 and finite
parameters; invalid values or normalization failure returns Left MathError.
-}
windowedSinc :: Window -> Int -> Double -> Double -> Either MathError [Double]
windowedSinc :: Window -> Int -> Double -> Double -> Either MathError [Double]
windowedSinc Window
kind Int
taps Double
cutoff Double
rate = do
    Bool -> String -> Either MathError ()
require
        ( Int
taps Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0
            Bool -> Bool -> Bool
&& Int -> Bool
forall a. Integral a => a -> Bool
odd Int
taps
            Bool -> Bool -> Bool
&& Int
taps Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
257
            Bool -> Bool -> Bool
&& Double -> Bool
finite Double
rate
            Bool -> Bool -> Bool
&& Double
rate Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0
            Bool -> Bool -> Bool
&& Double -> Bool
finite Double
cutoff
            Bool -> Bool -> Bool
&& Double
cutoff Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0
            Bool -> Bool -> Bool
&& Double
cutoff Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
rate Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2
        )
        String
"invalid windowed-sinc taps/cutoff/rate"
    weights <- Window -> Int -> Either MathError [Double]
window Window
kind Int
taps
    let center = Int
taps Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2
        fc = Double
cutoff Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
rate
        raw =
            (Double -> Int -> Double) -> [Double] -> [Int] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith
                (\Double
w Int
j -> Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
* (if Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 then Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
fc else Double -> Double
forall a. Floating a => a -> a
sin (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
fc Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
j) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
j)))
                [Double]
weights
                [-Int
center .. Int
center]
        scale = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double]
raw
    require (scale /= 0 && finite scale) "windowed sinc has zero/overflow DC gain"
    vector (map (/ scale) raw)

-- | Select indices congruent to offset modulo factor; this primitive performs no implicit anti-alias filtering.
downsample :: Int -> Int -> [Double] -> Either MathError [Double]
downsample :: Int -> Int -> [Double] -> Either MathError [Double]
downsample Int
factor Int
offset [Double]
xs = do
    Bool -> String -> Either MathError ()
require
        (Int
factor Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 Bool -> Bool -> Bool
&& Int
offset Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 Bool -> Bool -> Bool
&& Int
offset Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
factor Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite [Double]
xs)
        String
"invalid decimation factor/offset/data"
    [Double] -> Either MathError [Double]
forall a b. b -> Either a b
Right [Double
x | (Int
i, Double
x) <- [Int] -> [Double] -> [(Int, Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 :: Int ..] [Double]
xs, Int
i Int -> Int -> Int
forall a. Integral a => a -> a -> a
`mod` Int
factor Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
offset]

-- | Insert factor-1 zeros after each sample, with a checked 1048576-element output budget.
upsample :: Int -> [Double] -> Either MathError [Double]
upsample :: Int -> [Double] -> Either MathError [Double]
upsample Int
factor [Double]
xs = do
    Bool -> String -> Either MathError ()
require
        (Int
factor Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 Bool -> Bool -> Bool
&& Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
factor Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
forall a. Integral a => a -> Integer
toInteger ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs) Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
1048576 Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite [Double]
xs)
        String
"invalid upsampling factor/output budget"
    [Double] -> Either MathError [Double]
forall a b. b -> Either a b
Right ((Double -> [Double]) -> [Double] -> [Double]
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap (\Double
x -> Double
x Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate (Int
factor Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Double
0) [Double]
xs)

{- | Zero-extended centered polyphase FIR resampling. Only nonzero upsampled
phases are evaluated. Output has ceiling(N*up/down) samples, no group delay.
Kernel has odd length and is centered at index length/2; gain includes up.
-}
resamplePolyphase :: Int -> Int -> [Double] -> [Double] -> Either MathError [Double]
resamplePolyphase :: Int -> Int -> [Double] -> [Double] -> Either MathError [Double]
resamplePolyphase Int
up Int
down [Double]
kernel [Double]
xs = do
    Bool -> String -> Either MathError ()
require
        ( Int
up Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0
            Bool -> Bool -> Bool
&& Int
down Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0
            Bool -> Bool -> Bool
&& Bool -> Bool
not ([Double] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Double]
kernel)
            Bool -> Bool -> Bool
&& Int -> Bool
forall a. Integral a => a -> Bool
odd ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
kernel)
            Bool -> Bool -> Bool
&& Int -> Integer
forall a. Integral a => a -> Integer
toInteger ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
up Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
1048576
        )
        String
"invalid polyphase resampling factors/kernel/budget"
    _ <- [Double] -> Either MathError [Double]
vector ([Double]
kernel [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double]
xs)
    let output = (Int -> Integer
forall a. Integral a => a -> Integer
toInteger ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
up Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
down Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
1) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`div` Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
down
        center = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
kernel Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2
        at Int
k =
            [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum
                [ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
up Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* ([Double]
kernel [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j)
                | (Int
i, Double
x) <- [Int] -> [Double] -> [(Int, Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] [Double]
xs
                , let j :: Int
j = Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
down Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
up Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
center
                , Int
j Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0
                , Int
j Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
kernel
                ]
    vector (map at [0 .. fromInteger output - 1])

-- | Construct a Kaiser-windowed anti-alias/interpolation FIR, then evaluate only the necessary rational polyphase taps.
resample :: Int -> Int -> Int -> [Double] -> Either MathError [Double]
resample :: Int -> Int -> Int -> [Double] -> Either MathError [Double]
resample Int
up Int
down Int
taps [Double]
xs = do
    Bool -> String -> Either MathError ()
require (Int
up Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 Bool -> Bool -> Bool
&& Int
down Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 Bool -> Bool -> Bool
&& Int
up Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
1024 Bool -> Bool -> Bool
&& Int
down Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
1024) String
"invalid resampling factors"
    let cutoff :: Double
cutoff = Double
0.45 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
up Int
down)
    kernel <- Window -> Int -> Double -> Double -> Either MathError [Double]
windowedSinc (Double -> Window
Kaiser Double
8) Int
taps Double
cutoff Double
1
    resamplePolyphase up down kernel xs

-- | Centered replicate-padded low-pass before decimation by factor; cutoff is 0.45/factor cycles/sample.
antiAlias :: Int -> Int -> [Double] -> Either MathError [Double]
antiAlias :: Int -> Int -> [Double] -> Either MathError [Double]
antiAlias Int
factor Int
taps [Double]
xs = do
    Bool -> String -> Either MathError ()
require (Int
factor Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
1) String
"anti-alias factor must exceed one"
    kernel <- Window -> Int -> Double -> Double -> Either MathError [Double]
windowedSinc (Double -> Window
Kaiser Double
8) Int
taps (Double
0.45 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
factor) Double
1
    -- Centered replicate-padded convolution, followed by explicit decimation.
    localFilter taps (checked . dot kernel) xs

-- | Interior peak index, height, topographic prominence and half-prominence width in samples.
data Peak = Peak
    { Peak -> Int
peakIndex :: Int
    -- ^ Zero-based sample index of the detected peak.
    , Peak -> Double
peakHeight :: Double
    -- ^ Sample amplitude at peakIndex.
    , Peak -> Double
peakProminence :: Double
    -- ^ Peak height above its reference surrounding contour level.
    , Peak -> Double
peakWidth :: Double
    -- ^ Estimated width at the selected prominence height, measured in samples.
    }
    deriving (Peak -> Peak -> Bool
(Peak -> Peak -> Bool) -> (Peak -> Peak -> Bool) -> Eq Peak
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Peak -> Peak -> Bool
== :: Peak -> Peak -> Bool
$c/= :: Peak -> Peak -> Bool
/= :: Peak -> Peak -> Bool
Eq, Int -> Peak -> ShowS
[Peak] -> ShowS
Peak -> String
(Int -> Peak -> ShowS)
-> (Peak -> String) -> ([Peak] -> ShowS) -> Show Peak
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Peak -> ShowS
showsPrec :: Int -> Peak -> ShowS
$cshow :: Peak -> String
show :: Peak -> String
$cshowList :: [Peak] -> ShowS
showList :: [Peak] -> ShowS
Show)

{- | Interior local maxima; flat plateaus select the lower middle index. Width
is linearly interpolated at half prominence; endpoints are not peaks.
-}
peaks :: Double -> Int -> [Double] -> Either MathError [Peak]
peaks :: Double -> Int -> [Double] -> Either MathError [Peak]
peaks Double
minimumProminence Int
minimumDistance [Double]
xs = do
    Bool -> String -> Either MathError ()
require
        (Double -> Bool
finite Double
minimumProminence Bool -> Bool -> Bool
&& Double
minimumProminence Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0 Bool -> Bool -> Bool
&& Int
minimumDistance Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
1 Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite [Double]
xs)
        String
"invalid peak criteria/data"
    let n :: Int
n = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs
        candidates :: [(Int, Int)]
candidates =
            [ (Int
i, Int
j)
            | Int
i <- [Int
1 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
2]
            , [Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> [Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)
            , let j :: Int
j = Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ((Double -> Bool) -> [Double] -> [Double]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== [Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [Double]
xs))
            , Int
j Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1
            , [Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> [Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)
            ]
        make :: (Int, Int) -> Peak
make (Int
start, Int
end) =
            let i :: Int
i = (Int
start Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
end) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2
                h :: Double
h = [Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i
                left :: [Double]
left = [Double] -> [Double]
forall a. [a] -> [a]
reverse (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
take Int
start [Double]
xs)
                right :: [Double]
right = Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop (Int
end Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [Double]
xs
                leftB :: Double
leftB = [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum (Double
h Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: (Double -> Bool) -> [Double] -> [Double]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
h) [Double]
left)
                rightB :: Double
rightB = [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum (Double
h Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: (Double -> Bool) -> [Double] -> [Double]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
h) [Double]
right)
                prominence :: Double
prominence = Double
h Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
leftB Double
rightB
                level :: Double
level = Double
h Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
prominence Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2
                l :: Int
l = [Int] -> Int
forall a. HasCallStack => [a] -> a
last (Int
0 Int -> [Int] -> [Int]
forall a. a -> [a] -> [a]
: [Int
j | Int
j <- [Int
0 .. Int
start Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1], [Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
level])
                r :: Int
r = case [Int
j | Int
j <- [Int
end Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1], [Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
level] of Int
j : [Int]
_ -> Int
j; [] -> Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1
                interpolateCross :: Int -> Int -> Double
interpolateCross Int
a Int
b =
                    if [Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
a Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== [Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
b
                        then Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
a
                        else Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
level Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
a) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ ([Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
a) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
b Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
a)
                width :: Double
width = Int -> Int -> Double
interpolateCross (Int
r Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int
r Double -> Double -> Double
forall a. Num a => a -> a -> a
- Int -> Int -> Double
interpolateCross Int
l (Int
l Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)
             in Int -> Double -> Double -> Double -> Peak
Peak Int
i Double
h Double
prominence (Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
0 Double
width)
        ranked :: [Peak]
ranked = [Peak] -> [Peak]
sortByHeight ((Peak -> Bool) -> [Peak] -> [Peak]
forall a. (a -> Bool) -> [a] -> [a]
filter ((Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
minimumProminence) (Double -> Bool) -> (Peak -> Double) -> Peak -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Peak -> Double
peakProminence) (((Int, Int) -> Peak) -> [(Int, Int)] -> [Peak]
forall a b. (a -> b) -> [a] -> [b]
map (Int, Int) -> Peak
make [(Int, Int)]
candidates))
        accepted :: [Peak]
accepted =
            ([Peak] -> Peak -> [Peak]) -> [Peak] -> [Peak] -> [Peak]
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl
                ( \[Peak]
acc Peak
p -> if (Peak -> Bool) -> [Peak] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (\Peak
q -> Int -> Int
forall a. Num a => a -> a
abs (Peak -> Int
peakIndex Peak
p Int -> Int -> Int
forall a. Num a => a -> a -> a
- Peak -> Int
peakIndex Peak
q) Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
minimumDistance) [Peak]
acc then Peak
p Peak -> [Peak] -> [Peak]
forall a. a -> [a] -> [a]
: [Peak]
acc else [Peak]
acc
                )
                []
                [Peak]
ranked
    _ <- [Double] -> Either MathError [Double]
vector ([[Double]] -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [[Peak -> Double
peakHeight Peak
p, Peak -> Double
peakProminence Peak
p, Peak -> Double
peakWidth Peak
p] | Peak
p <- [Peak]
accepted])
    Right (sortByIndex accepted)
  where
    sortByHeight :: [Peak] -> [Peak]
sortByHeight = [Peak] -> [Peak]
forall a. [a] -> [a]
reverse ([Peak] -> [Peak]) -> ([Peak] -> [Peak]) -> [Peak] -> [Peak]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Peak] -> [Peak]
sortOnHeight
    sortOnHeight :: [Peak] -> [Peak]
sortOnHeight =
        (Peak -> Peak -> Ordering) -> [Peak] -> [Peak]
forall {t}. (t -> t -> Ordering) -> [t] -> [t]
sortByLocal
            (\Peak
a Peak
b -> (Double, Int) -> (Double, Int) -> Ordering
forall a. Ord a => a -> a -> Ordering
compare (Peak -> Double
peakHeight Peak
a, Int -> Int
forall a. Num a => a -> a
negate (Peak -> Int
peakIndex Peak
a)) (Peak -> Double
peakHeight Peak
b, Int -> Int
forall a. Num a => a -> a
negate (Peak -> Int
peakIndex Peak
b)))
    sortByIndex :: [Peak] -> [Peak]
sortByIndex = (Peak -> Peak -> Ordering) -> [Peak] -> [Peak]
forall {t}. (t -> t -> Ordering) -> [t] -> [t]
sortByLocal (\Peak
a Peak
b -> Int -> Int -> Ordering
forall a. Ord a => a -> a -> Ordering
compare (Peak -> Int
peakIndex Peak
a) (Peak -> Int
peakIndex Peak
b))
    sortByLocal :: (t -> t -> Ordering) -> [t] -> [t]
sortByLocal t -> t -> Ordering
_ [] = []
    sortByLocal t -> t -> Ordering
compareValues (t
p : [t]
ps) =
        (t -> t -> Ordering) -> [t] -> [t]
sortByLocal t -> t -> Ordering
compareValues [t
x | t
x <- [t]
ps, t -> t -> Ordering
compareValues t
x t
p Ordering -> Ordering -> Bool
forall a. Eq a => a -> a -> Bool
== Ordering
LT]
            [t] -> [t] -> [t]
forall a. [a] -> [a] -> [a]
++ [t
p]
            [t] -> [t] -> [t]
forall a. [a] -> [a] -> [a]
++ (t -> t -> Ordering) -> [t] -> [t]
sortByLocal t -> t -> Ordering
compareValues [t
x | t
x <- [t]
ps, t -> t -> Ordering
compareValues t
x t
p Ordering -> Ordering -> Bool
forall a. Eq a => a -> a -> Bool
/= Ordering
LT]

-- | Interpolated opposite-sign transitions; zero plateaus count once, leading/trailing zeros are ignored.
zeroCrossings :: [Double] -> Either MathError [Double]
zeroCrossings :: [Double] -> Either MathError [Double]
zeroCrossings [Double]
xs = do
    _ <- [Double] -> Either MathError [Double]
vector [Double]
xs
    -- Exact zero plateaus count once, located at their first zero if a sign
    -- change occurs across the plateau; leading/trailing zeros are ignored.
    let nonzero = [(Int
i, Double
x) | (Int
i, Double
x) <- [Int] -> [Double] -> [(Int, Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 :: Int ..] [Double]
xs, Double
x Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
/= Double
0]
    let crossing a
a a
b = let scale :: a
scale = a -> a -> a
forall a. Ord a => a -> a -> a
max (a -> a
forall a. Num a => a -> a
abs a
a) (a -> a
forall a. Num a => a -> a
abs a
b); aa :: a
aa = a -> a
forall a. Num a => a -> a
abs a
a a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
scale; bb :: a
bb = a -> a
forall a. Num a => a -> a
abs a
b a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
scale in a
aa a -> a -> a
forall a. Fractional a => a -> a -> a
/ (a
aa a -> a -> a
forall a. Num a => a -> a -> a
+ a
bb)
    vector
        [ if j == i + 1 then fromIntegral i + crossing a b else fromIntegral (i + 1)
        | ((i, a), (j, b)) <- zip nonzero (drop 1 nonzero)
        , (a < 0 && b > 0) || (a > 0 && b < 0)
        ]

-- | Magnitude of the periodic analytic-signal representation.
envelope :: [Double] -> Either MathError [Double]
envelope :: [Double] -> Either MathError [Double]
envelope [Double]
xs = (Complex Double -> Double) -> [Complex Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Complex Double -> Double
forall a. RealFloat a => Complex a -> a
magnitude ([Complex Double] -> [Double])
-> Either MathError [Complex Double] -> Either MathError [Double]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> [Double] -> Either MathError [Complex Double]
F.analyticSignal [Double]
xs

-- | Remove nearest multiples of 2*pi from adjacent phase jumps; reject overflowing finite differences before integer conversion.
unwrapPhase :: [Double] -> Either MathError [Double]
unwrapPhase :: [Double] -> Either MathError [Double]
unwrapPhase [Double]
xs = do
    _ <- [Double] -> Either MathError [Double]
vector [Double]
xs
    case xs of
        [] -> [Double] -> Either MathError [Double]
forall a b. b -> Either a b
Right []
        Double
x : [Double]
rest -> do
            let step :: ([Double], Double, Double)
-> Double -> Either MathError ([Double], Double, Double)
step ([Double]
acc, Double
previous, Double
offset) Double
v = do
                    delta <- Double -> Either MathError Double
checked (Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
previous)
                    let wrapped = Double
delta Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Integer -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Double -> Integer
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
floor ((Double
delta Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
forall a. Floating a => a
pi) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
forall a. Floating a => a
pi)) :: Integer)
                        shift = if Double
wrapped Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== (-Double
forall a. Floating a => a
pi) Bool -> Bool -> Bool
&& Double
delta Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 then Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
delta else Double
wrapped Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
delta
                    next <- checked (offset + shift)
                    value <- checked (v + next)
                    Right (value : acc, v, next)
            (ys, _, _) <- (([Double], Double, Double)
 -> Double -> Either MathError ([Double], Double, Double))
-> ([Double], Double, Double)
-> [Double]
-> Either MathError ([Double], Double, Double)
forall (t :: * -> *) (m :: * -> *) b a.
(Foldable t, Monad m) =>
(b -> a -> m b) -> b -> t a -> m b
foldM ([Double], Double, Double)
-> Double -> Either MathError ([Double], Double, Double)
step ([Double
x], Double
x, Double
0) [Double]
rest
            vector (reverse ys)

-- | Subtract a constant mean or a fitted linear trend in sample index.
detrend :: Bool -> [Double] -> Either MathError [Double]
detrend :: Bool -> [Double] -> Either MathError [Double]
detrend Bool
linear [Double]
xs = do
    Bool -> String -> Either MathError ()
require (Bool -> Bool
not ([Double] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Double]
xs)) String
"detrend needs samples"
    _ <- [Double] -> Either MathError [Double]
vector [Double]
xs
    let n = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs); mean = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
n) [Double]
xs)
    if not linear || length xs == 1
        then vector (map (subtract mean) xs)
        else do
            fit <- I.polynomialFit 1 (zip [0 ..] xs)
            let cs = Approximation -> [Double]
I.approximationCoefficients Approximation
fit
            vector [x - cs !! 0 - cs !! 1 * fromIntegral j | (j, x) <- zip [0 :: Int ..] xs]

-- | Divide by maximum absolute amplitude; reject an empty or identically zero signal.
normalize :: [Double] -> Either MathError [Double]
normalize :: [Double] -> Either MathError [Double]
normalize [Double]
xs = do
    Bool -> String -> Either MathError ()
require (Bool -> Bool
not ([Double] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Double]
xs)) String
"normalization needs samples"
    _ <- [Double] -> Either MathError [Double]
vector [Double]
xs
    let scale = [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Double -> Double
forall a. Num a => a -> a
abs [Double]
xs)
    require (scale > 0) "cannot normalize zero signal"
    vector (map (/ scale) xs)

-- | Centered odd-width arithmetic mean with replicate endpoint padding.
smooth :: Int -> [Double] -> Either MathError [Double]
smooth :: Int -> [Double] -> Either MathError [Double]
smooth Int
width = Int
-> ([Double] -> Either MathError Double)
-> [Double]
-> Either MathError [Double]
localFilter Int
width (\[Double]
ys -> Double -> Either MathError Double
checked ([Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double]
ys Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
ys)))

-- | Robust Gaussian noise sigma from median absolute adjacent differences divided by sqrt(2)*Phi^-1(0.75).
noiseEstimate :: [Double] -> Either MathError Double
noiseEstimate :: [Double] -> Either MathError Double
noiseEstimate [Double]
xs = do
    Bool -> String -> Either MathError ()
require ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
2) String
"noise estimate needs two samples"
    _ <- [Double] -> Either MathError [Double]
vector [Double]
xs
    let differences = (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Double
a Double
b -> Double -> Double
forall a. Num a => a -> a
abs (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a)) [Double]
xs (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
1 [Double]
xs)
    _ <- vector differences
    checked (sort differences !! (length differences `div` 2) / (0.6744897501960817 * sqrt 2))

-- | 20*(log10(signal RMS)-log10(noise RMS)) in dB; both energies must be positive.
signalToNoise :: [Double] -> [Double] -> Either MathError Double
signalToNoise :: [Double] -> [Double] -> Either MathError Double
signalToNoise [Double]
signal [Double]
noise = do
    Bool -> String -> Either MathError ()
require ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
signal Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
noise Bool -> Bool -> Bool
&& Bool -> Bool
not ([Double] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Double]
signal)) String
"SNR length mismatch/empty"
    s <- [Double] -> Either MathError Double
rms [Double]
signal
    n <- rms noise
    require (s > 0 && n > 0) "SNR requires positive signal and noise energies"
    checked (20 * (logBase 10 s - logBase 10 n))

-- | Scaled root-mean-square calculation that avoids overflow when the RMS itself is representable.
rms :: [Double] -> Either MathError Double
rms :: [Double] -> Either MathError Double
rms [Double]
xs = do
    Bool -> String -> Either MathError ()
require (Bool -> Bool
not ([Double] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Double]
xs)) String
"RMS requires samples"
    _ <- [Double] -> Either MathError [Double]
vector [Double]
xs
    let scale = [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Double -> Double
forall a. Num a => a -> a
abs [Double]
xs)
    checked
        ( if scale == 0
            then 0
            else scale * sqrt (sum [(x / scale) ^ (2 :: Int) | x <- xs] / fromIntegral (length xs))
        )

{- | Mean and symmetric positive semidefinite covariance, dimension 1..64.
User-supplied transition/measurement callbacks must be finite, pure and total.
-}
data GaussianState = GaussianState [Double] [[Double]] deriving (GaussianState -> GaussianState -> Bool
(GaussianState -> GaussianState -> Bool)
-> (GaussianState -> GaussianState -> Bool) -> Eq GaussianState
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: GaussianState -> GaussianState -> Bool
== :: GaussianState -> GaussianState -> Bool
$c/= :: GaussianState -> GaussianState -> Bool
/= :: GaussianState -> GaussianState -> Bool
Eq, Int -> GaussianState -> ShowS
[GaussianState] -> ShowS
GaussianState -> String
(Int -> GaussianState -> ShowS)
-> (GaussianState -> String)
-> ([GaussianState] -> ShowS)
-> Show GaussianState
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> GaussianState -> ShowS
showsPrec :: Int -> GaussianState -> ShowS
$cshow :: GaussianState -> String
show :: GaussianState -> String
$cshowList :: [GaussianState] -> ShowS
showList :: [GaussianState] -> ShowS
Show)

{- | Validate a mean and symmetric positive semidefinite covariance; Cholesky
roundoff tolerance is 1e-12*max(1,diagonal). Accepted tiny negative pivots
are clipped, and the stored covariance is reconstructed from that factor.
Thus a tolerated negative input variance never escapes through stateCovariance.
-}
gaussianState :: [Double] -> [[Double]] -> Either MathError GaussianState
gaussianState :: [Double] -> [[Double]] -> Either MathError GaussianState
gaussianState [Double]
mean [[Double]]
covariance = do
    Bool -> String -> Either MathError ()
require (Bool -> Bool
not ([Double] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Double]
mean) Bool -> Bool -> Bool
&& [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
mean Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
64) String
"Gaussian state dimension must be 1..64"
    _ <- [Double] -> Either MathError [Double]
vector [Double]
mean
    factor <- covarianceCheck (length mean) covariance
    let accepted = [[Double]] -> [[Double]] -> [[Double]]
mm [[Double]]
factor ([[Double]] -> [[Double]]
forall a. [[a]] -> [[a]]
tr [[Double]]
factor)
    _ <- vector (concat accepted)
    Right (GaussianState mean accepted)

-- | Read the posterior mean vector.
stateMean :: GaussianState -> [Double]
stateMean :: GaussianState -> [Double]
stateMean (GaussianState [Double]
mean [[Double]]
_) = [Double]
mean

-- | Read the symmetric posterior covariance matrix.
stateCovariance :: GaussianState -> [[Double]]
stateCovariance :: GaussianState -> [[Double]]
stateCovariance (GaussianState [Double]
_ [[Double]]
covariance) = [[Double]]
covariance

symmetrize :: [[Double]] -> [[Double]]
symmetrize :: [[Double]] -> [[Double]]
symmetrize [[Double]]
a = ([Double] -> [Double] -> [Double])
-> [[Double]] -> [[Double]] -> [[Double]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Double
x Double
y -> Double
x Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
y Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2)) [[Double]]
a ([[Double]] -> [[Double]]
forall a. [[a]] -> [[a]]
tr [[Double]]
a)
matrixCheck :: Int -> Int -> [[Double]] -> Either MathError ()
matrixCheck :: Int -> Int -> [[Double]] -> Either MathError ()
matrixCheck Int
rows Int
cols [[Double]]
a = do
    Bool -> String -> Either MathError ()
require
        (Int
rows Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 Bool -> Bool -> Bool
&& Int
cols Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 Bool -> Bool -> Bool
&& Int
rows Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
64 Bool -> Bool -> Bool
&& Int
cols Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
64 Bool -> Bool -> Bool
&& [[Double]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
a Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
rows Bool -> Bool -> Bool
&& ([Double] -> Bool) -> [[Double]] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all ((Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
cols) (Int -> Bool) -> ([Double] -> Int) -> [Double] -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length) [[Double]]
a)
        String
"invalid Kalman matrix dimensions (1..64)"
    _ <- [Double] -> Either MathError [Double]
vector ([[Double]] -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [[Double]]
a)
    Right ()
covarianceCheck :: Int -> [[Double]] -> Either MathError [[Double]]
covarianceCheck :: Int -> [[Double]] -> Either MathError [[Double]]
covarianceCheck Int
n [[Double]]
a = do
    Int -> Int -> [[Double]] -> Either MathError ()
matrixCheck Int
n Int
n [[Double]]
a
    Bool -> String -> Either MathError ()
require
        ( [Bool] -> Bool
forall (t :: * -> *). Foldable t => t Bool -> Bool
and
            [ Double -> Double
forall a. Num a => a -> a
abs ([[Double]]
a [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j Double -> Double -> Double
forall a. Num a => a -> a -> a
- [[Double]]
a [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
j [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
1e-10 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
1 (Double -> Double -> Double
forall a. Ord a => a -> a -> a
max (Double -> Double
forall a. Num a => a -> a
abs ([[Double]]
a [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j)) (Double -> Double
forall a. Num a => a -> a
abs ([[Double]]
a [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
j [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i)))
            | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
            , Int
j <- [Int
0 .. Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
            ]
        )
        String
"covariance must be symmetric"
    [[Double]] -> Either MathError [[Double]]
cholesky [[Double]]
a
cholesky :: [[Double]] -> Either MathError [[Double]]
cholesky :: [[Double]] -> Either MathError [[Double]]
cholesky [[Double]]
a = ([[Double]] -> Int -> Either MathError [[Double]])
-> [[Double]] -> [Int] -> Either MathError [[Double]]
forall (t :: * -> *) (m :: * -> *) b a.
(Foldable t, Monad m) =>
(b -> a -> m b) -> b -> t a -> m b
foldM [[Double]] -> Int -> Either MathError [[Double]]
row [] [Int
0 .. [[Double]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[Double]]
a Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
  where
    row :: [[Double]] -> Int -> Either MathError [[Double]]
row [[Double]]
previous Int
i = do
        values <- ([Double] -> Int -> Either MathError [Double])
-> [Double] -> [Int] -> Either MathError [Double]
forall (t :: * -> *) (m :: * -> *) b a.
(Foldable t, Monad m) =>
(b -> a -> m b) -> b -> t a -> m b
foldM ([[Double]] -> Int -> [Double] -> Int -> Either MathError [Double]
entry [[Double]]
previous Int
i) [] [Int
0 .. Int
i]
        Right (previous ++ [values ++ replicate (length a - i - 1) 0])
    entry :: [[Double]] -> Int -> [Double] -> Int -> Either MathError [Double]
entry [[Double]]
previous Int
i [Double]
current Int
j = do
        let residual :: Double
residual =
                [[Double]]
a [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j
                    Double -> Double -> Double
forall a. Num a => a -> a -> a
- [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [[Double]
current [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* (if Int
j Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
i then [Double]
current [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
k else [[Double]]
previous [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
j [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
k) | Int
k <- [Int
0 .. Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
            scale :: Double
scale = Double
1e-12 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
1 (Double -> Double
forall a. Num a => a -> a
abs ([[Double]]
a [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i))
        value <-
            if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j
                then if Double
residual Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< (-Double
scale) then MathError -> Either MathError Double
forall a b. a -> Either a b
Left (Int -> MathError
NotPositiveDefinite Int
i) else Double -> Either MathError Double
checked (Double -> Double
forall a. Floating a => a -> a
sqrt (Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
0 Double
residual))
                else
                    let pivot :: Double
pivot = [[Double]]
previous [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
j [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j
                     in if Double
pivot Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0
                            then Bool -> String -> Either MathError ()
require (Double -> Double
forall a. Num a => a -> a
abs Double
residual Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
scale) String
"singular inconsistent covariance" Either MathError ()
-> Either MathError Double -> Either MathError Double
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Double -> Either MathError Double
forall a b. b -> Either a b
Right Double
0
                            else Double -> Either MathError Double
checked (Double
residual Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
pivot)
        Right (current ++ [value])
addM, subM :: [[Double]] -> [[Double]] -> [[Double]]
addM :: [[Double]] -> [[Double]] -> [[Double]]
addM = ([Double] -> [Double] -> [Double])
-> [[Double]] -> [[Double]] -> [[Double]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(+))
subM :: [[Double]] -> [[Double]] -> [[Double]]
subM = ([Double] -> [Double] -> [Double])
-> [[Double]] -> [[Double]] -> [[Double]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-))
outer :: [Double] -> [Double] -> [[Double]]
outer :: [Double] -> [Double] -> [[Double]]
outer [Double]
x [Double]
y = [(Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
*) [Double]
y | Double
a <- [Double]
x]
solveRight :: [[Double]] -> [[Double]] -> Either MathError [[Double]]
solveRight :: [[Double]] -> [[Double]] -> Either MathError [[Double]]
solveRight [[Double]]
numerator [[Double]]
denominator = ([Double] -> Either MathError [Double])
-> [[Double]] -> Either MathError [[Double]]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
forall (m :: * -> *) a b. Monad m => (a -> m b) -> [a] -> m [b]
mapM ([[Double]] -> [Double] -> Either MathError [Double]
solve ([[Double]] -> [[Double]]
forall a. [[a]] -> [[a]]
tr [[Double]]
denominator)) [[Double]]
numerator

{- | One linear Gaussian prediction+measurement update. Uses Joseph covariance
update to preserve positive semidefiniteness under roundoff.
-}
kalman ::
    [[Double]] ->
    [[Double]] ->
    [[Double]] ->
    [[Double]] ->
    [Double] ->
    GaussianState ->
    Either MathError GaussianState
kalman :: [[Double]]
-> [[Double]]
-> [[Double]]
-> [[Double]]
-> [Double]
-> GaussianState
-> Either MathError GaussianState
kalman [[Double]]
transition [[Double]]
observation [[Double]]
processNoise [[Double]]
measurementNoise [Double]
measured (GaussianState [Double]
mean [[Double]]
covariance) = do
    let n :: Int
n = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
mean; m :: Int
m = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
measured
    Int -> Int -> [[Double]] -> Either MathError ()
matrixCheck Int
n Int
n [[Double]]
transition
    Int -> Int -> [[Double]] -> Either MathError ()
matrixCheck Int
m Int
n [[Double]]
observation
    _ <- Int -> [[Double]] -> Either MathError [[Double]]
covarianceCheck Int
n [[Double]]
processNoise
    _ <- covarianceCheck m measurementNoise
    _ <- vector measured
    let predicted = [[Double]] -> [Double] -> [Double]
mv [[Double]]
transition [Double]
mean
        pp = [[Double]] -> [[Double]] -> [[Double]]
addM ([[Double]] -> [[Double]] -> [[Double]]
mm ([[Double]] -> [[Double]] -> [[Double]]
mm [[Double]]
transition [[Double]]
covariance) ([[Double]] -> [[Double]]
forall a. [[a]] -> [[a]]
tr [[Double]]
transition)) [[Double]]
processNoise
    linearUpdate observation measurementNoise measured predicted pp

linearUpdate ::
    [[Double]] -> [[Double]] -> [Double] -> [Double] -> [[Double]] -> Either MathError GaussianState
linearUpdate :: [[Double]]
-> [[Double]]
-> [Double]
-> [Double]
-> [[Double]]
-> Either MathError GaussianState
linearUpdate [[Double]]
h [[Double]]
r [Double]
measured [Double]
predicted [[Double]]
pp = do
    let innovation :: [Double]
innovation = (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) [Double]
measured ([[Double]] -> [Double] -> [Double]
mv [[Double]]
h [Double]
predicted)
        pht :: [[Double]]
pht = [[Double]] -> [[Double]] -> [[Double]]
mm [[Double]]
pp ([[Double]] -> [[Double]]
forall a. [[a]] -> [[a]]
tr [[Double]]
h)
        s :: [[Double]]
s = [[Double]] -> [[Double]] -> [[Double]]
addM ([[Double]] -> [[Double]] -> [[Double]]
mm [[Double]]
h [[Double]]
pht) [[Double]]
r
    gain <- [[Double]] -> [[Double]] -> Either MathError [[Double]]
solveRight [[Double]]
pht [[Double]]
s
    let next = (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(+) [Double]
predicted ([[Double]] -> [Double] -> [Double]
mv [[Double]]
gain [Double]
innovation)
        factor = [[Double]] -> [[Double]] -> [[Double]]
subM (Int -> [[Double]]
eye ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
predicted)) ([[Double]] -> [[Double]] -> [[Double]]
mm [[Double]]
gain [[Double]]
h)
        covariance = [[Double]] -> [[Double]] -> [[Double]]
addM ([[Double]] -> [[Double]] -> [[Double]]
mm ([[Double]] -> [[Double]] -> [[Double]]
mm [[Double]]
factor [[Double]]
pp) ([[Double]] -> [[Double]]
forall a. [[a]] -> [[a]]
tr [[Double]]
factor)) ([[Double]] -> [[Double]] -> [[Double]]
mm ([[Double]] -> [[Double]] -> [[Double]]
mm [[Double]]
gain [[Double]]
r) ([[Double]] -> [[Double]]
forall a. [[a]] -> [[a]]
tr [[Double]]
gain))
    gaussianState next covariance

jacobianFinite :: ([Double] -> [Double]) -> [Double] -> Int -> Either MathError [[Double]]
jacobianFinite :: ([Double] -> [Double])
-> [Double] -> Int -> Either MathError [[Double]]
jacobianFinite [Double] -> [Double]
f [Double]
point Int
output = do
    columns <-
        (Int -> Either MathError [Double])
-> [Int] -> Either MathError [[Double]]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
forall (m :: * -> *) a b. Monad m => (a -> m b) -> [a] -> m [b]
mapM
            ( \Int
j -> do
                let h :: Double
h = Double
1e-5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
1 (Double -> Double
forall a. Num a => a -> a
abs ([Double]
point [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j))
                    plus :: [Double]
plus = [Double] -> [Double]
f (Int -> Double -> [Double] -> [Double]
forall a. Int -> a -> [a] -> [a]
replace Int
j ([Double]
point [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
h) [Double]
point)
                    minus :: [Double]
minus = [Double] -> [Double]
f (Int -> Double -> [Double] -> [Double]
forall a. Int -> a -> [a] -> [a]
replace Int
j ([Double]
point [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
h) [Double]
point)
                Bool -> String -> Either MathError ()
require ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
plus Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
output Bool -> Bool -> Bool
&& [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
minus Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
output) String
"Kalman callback changed dimension"
                [Double] -> Either MathError [Double]
vector ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Double
a Double
b -> (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
b) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
h)) [Double]
plus [Double]
minus)
            )
            [Int
0 .. [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
point Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
    Right (tr columns)

-- | EKF with checked central-difference Jacobians (step 1e-5*max(1,abs x)).
extendedKalman ::
    ([Double] -> [Double]) ->
    ([Double] -> [Double]) ->
    [[Double]] ->
    [[Double]] ->
    [Double] ->
    GaussianState ->
    Either MathError GaussianState
extendedKalman :: ([Double] -> [Double])
-> ([Double] -> [Double])
-> [[Double]]
-> [[Double]]
-> [Double]
-> GaussianState
-> Either MathError GaussianState
extendedKalman [Double] -> [Double]
transition [Double] -> [Double]
measurement [[Double]]
q [[Double]]
r [Double]
measured (GaussianState [Double]
mean [[Double]]
covariance) = do
    let n :: Int
n = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
mean; m :: Int
m = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
measured; predicted :: [Double]
predicted = [Double] -> [Double]
transition [Double]
mean
    Bool -> String -> Either MathError ()
require ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
predicted Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
n Bool -> Bool -> Bool
&& Int
m Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0) String
"EKF callback/output dimension mismatch"
    _ <- [Double] -> Either MathError [Double]
vector ([Double]
predicted [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double]
measured)
    _ <- covarianceCheck n q
    _ <- covarianceCheck m r
    f <- jacobianFinite transition mean n
    h <- jacobianFinite measurement predicted m
    let expected = [Double] -> [Double]
measurement [Double]
predicted
    require (length expected == m) "EKF observation dimension mismatch"
    _ <- vector expected
    let pp = [[Double]] -> [[Double]] -> [[Double]]
addM ([[Double]] -> [[Double]] -> [[Double]]
mm ([[Double]] -> [[Double]] -> [[Double]]
mm [[Double]]
f [[Double]]
covariance) ([[Double]] -> [[Double]]
forall a. [[a]] -> [[a]]
tr [[Double]]
f)) [[Double]]
q
        -- linearUpdate subtracts H*x, so shift observed data by H*x-h(x).
        adjusted = (Double -> Double -> Double -> Double)
-> [Double] -> [Double] -> [Double] -> [Double]
forall a b c d. (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
zipWith3 (\Double
z Double
linear Double
nonlinear -> Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
linear Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
nonlinear) [Double]
measured ([[Double]] -> [Double] -> [Double]
mv [[Double]]
h [Double]
predicted) [Double]
expected
    linearUpdate h r adjusted predicted pp

{- | Scaled UKF. alpha>0, beta>=0, n+kappa>0, lambda+n>0. Process
noise is incorporated before regenerating measurement sigma points.
-}
unscentedKalman ::
    Double ->
    Double ->
    Double ->
    ([Double] -> [Double]) ->
    ([Double] -> [Double]) ->
    [[Double]] ->
    [[Double]] ->
    [Double] ->
    GaussianState ->
    Either MathError GaussianState
unscentedKalman :: Double
-> Double
-> Double
-> ([Double] -> [Double])
-> ([Double] -> [Double])
-> [[Double]]
-> [[Double]]
-> [Double]
-> GaussianState
-> Either MathError GaussianState
unscentedKalman Double
alpha Double
beta Double
kappa [Double] -> [Double]
transition [Double] -> [Double]
measurement [[Double]]
q [[Double]]
r [Double]
measured GaussianState
state = do
    let n :: Int
n = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length (GaussianState -> [Double]
stateMean GaussianState
state)
        m :: Int
m = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
measured
        lambda :: Double
lambda = Double
alpha Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
alpha Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
kappa) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
        c :: Double
c = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
lambda
    Bool -> String -> Either MathError ()
require
        ((Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Double -> Bool
finite [Double
alpha, Double
beta, Double
kappa, Double
c] Bool -> Bool -> Bool
&& Double
alpha Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& Double
beta Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0 Bool -> Bool -> Bool
&& Double
c Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
1e-12 Bool -> Bool -> Bool
&& Int
m Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0)
        String
"invalid UKF scaling/output dimension"
    _ <- Int -> [[Double]] -> Either MathError [[Double]]
covarianceCheck Int
n [[Double]]
q
    _ <- covarianceCheck m r
    _ <- vector measured
    let weightsMean = Double
lambda Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
c Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate (Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
n) (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
c))
        weightsCov = (Double
lambda Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
alpha Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
alpha Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
beta) Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate (Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
n) (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
c))
        sigma (GaussianState [Double]
mu [[Double]]
p) = do
            l <- Int -> [[Double]] -> Either MathError [[Double]]
covarianceCheck Int
n [[Double]]
p
            let columns = ([Double] -> [Double]) -> [[Double]] -> [[Double]]
forall a b. (a -> b) -> [a] -> [b]
map ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Double
forall a. Floating a => a -> a
sqrt Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
*)) ([[Double]] -> [[Double]]
forall a. [[a]] -> [[a]]
tr [[Double]]
l)
            Right (mu : [zipWith (+) mu col | col <- columns] ++ [zipWith (-) mu col | col <- columns])
        transform Int
output a -> [Double]
f [a]
points = do
            let values :: [[Double]]
values = (a -> [Double]) -> [a] -> [[Double]]
forall a b. (a -> b) -> [a] -> [b]
map a -> [Double]
f [a]
points
            Bool -> String -> Either MathError ()
require (([Double] -> Bool) -> [[Double]] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all ((Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
output) (Int -> Bool) -> ([Double] -> Int) -> [Double] -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length) [[Double]]
values) String
"UKF callback changed dimension"
            _ <- [Double] -> Either MathError [Double]
vector ([[Double]] -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [[Double]]
values)
            let mu = [[Double]] -> [Double] -> [Double]
mv ([[Double]] -> [[Double]]
forall a. [[a]] -> [[a]]
tr [[Double]]
values) [Double]
weightsMean
                covariance =
                    ([[Double]] -> [[Double]] -> [[Double]])
-> [[Double]] -> [[[Double]]] -> [[Double]]
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl
                        [[Double]] -> [[Double]] -> [[Double]]
addM
                        (Int -> [Double] -> [[Double]]
forall a. Int -> a -> [a]
replicate Int
output (Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate Int
output Double
0))
                        [([Double] -> [Double]) -> [[Double]] -> [[Double]]
forall a b. (a -> b) -> [a] -> [b]
map ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
*)) ([Double] -> [Double] -> [[Double]]
outer ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) [Double]
v [Double]
mu) ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) [Double]
v [Double]
mu)) | (Double
w, [Double]
v) <- [Double] -> [[Double]] -> [(Double, [Double])]
forall a b. [a] -> [b] -> [(a, b)]
zip [Double]
weightsCov [[Double]]
values]
            _ <- vector (mu ++ concat covariance)
            Right (mu, covariance, values)
    points <- sigma state
    (pm, p0, _) <- transform n transition points
    predicted <- gaussianState pm (addM p0 q)
    observedPoints <- sigma predicted
    (zm, s0, zvalues) <- transform m measurement observedPoints
    let innovationCov = [[Double]] -> [[Double]] -> [[Double]]
addM [[Double]]
s0 [[Double]]
r
        cross =
            ([[Double]] -> [[Double]] -> [[Double]])
-> [[Double]] -> [[[Double]]] -> [[Double]]
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl
                [[Double]] -> [[Double]] -> [[Double]]
addM
                (Int -> [Double] -> [[Double]]
forall a. Int -> a -> [a]
replicate Int
n (Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate Int
m Double
0))
                [ ([Double] -> [Double]) -> [[Double]] -> [[Double]]
forall a b. (a -> b) -> [a] -> [b]
map ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
*)) ([Double] -> [Double] -> [[Double]]
outer ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) [Double]
x [Double]
pm) ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) [Double]
z [Double]
zm))
                | (Double
w, ([Double]
x, [Double]
z)) <- [Double]
-> [([Double], [Double])] -> [(Double, ([Double], [Double]))]
forall a b. [a] -> [b] -> [(a, b)]
zip [Double]
weightsCov ([[Double]] -> [[Double]] -> [([Double], [Double])]
forall a b. [a] -> [b] -> [(a, b)]
zip [[Double]]
observedPoints [[Double]]
zvalues)
                ]
    gain <- solveRight cross innovationCov
    let mean = (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(+) [Double]
pm ([[Double]] -> [Double] -> [Double]
mv [[Double]]
gain ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) [Double]
measured [Double]
zm))
        covariance = [[Double]] -> [[Double]] -> [[Double]]
subM (GaussianState -> [[Double]]
stateCovariance GaussianState
predicted) ([[Double]] -> [[Double]] -> [[Double]]
mm ([[Double]] -> [[Double]] -> [[Double]]
mm [[Double]]
gain [[Double]]
innovationCov) ([[Double]] -> [[Double]]
forall a. [[a]] -> [[a]]
tr [[Double]]
gain))
    gaussianState mean (symmetrize covariance)

{- | Orthogonal Haar, 4-tap db2, 8-tap db4, and reversible CDF 5/3 lifting.
CDF coefficients use the unnormalized lifting convention (constant detail=0).
-}
data Wavelet = Haar | Daubechies2 | Daubechies4 | Biorthogonal53 deriving (Wavelet -> Wavelet -> Bool
(Wavelet -> Wavelet -> Bool)
-> (Wavelet -> Wavelet -> Bool) -> Eq Wavelet
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Wavelet -> Wavelet -> Bool
== :: Wavelet -> Wavelet -> Bool
$c/= :: Wavelet -> Wavelet -> Bool
/= :: Wavelet -> Wavelet -> Bool
Eq, Int -> Wavelet -> ShowS
[Wavelet] -> ShowS
Wavelet -> String
(Int -> Wavelet -> ShowS)
-> (Wavelet -> String) -> ([Wavelet] -> ShowS) -> Show Wavelet
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Wavelet -> ShowS
showsPrec :: Int -> Wavelet -> ShowS
$cshow :: Wavelet -> String
show :: Wavelet -> String
$cshowList :: [Wavelet] -> ShowS
showList :: [Wavelet] -> ShowS
Show)

-- | Opaque multilevel approximation and detail coefficients, retaining family and original length.
data WaveletTree = WaveletTree Wavelet Int [Double] [[Double]] deriving (WaveletTree -> WaveletTree -> Bool
(WaveletTree -> WaveletTree -> Bool)
-> (WaveletTree -> WaveletTree -> Bool) -> Eq WaveletTree
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: WaveletTree -> WaveletTree -> Bool
== :: WaveletTree -> WaveletTree -> Bool
$c/= :: WaveletTree -> WaveletTree -> Bool
/= :: WaveletTree -> WaveletTree -> Bool
Eq, Int -> WaveletTree -> ShowS
[WaveletTree] -> ShowS
WaveletTree -> String
(Int -> WaveletTree -> ShowS)
-> (WaveletTree -> String)
-> ([WaveletTree] -> ShowS)
-> Show WaveletTree
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> WaveletTree -> ShowS
showsPrec :: Int -> WaveletTree -> ShowS
$cshow :: WaveletTree -> String
show :: WaveletTree -> String
$cshowList :: [WaveletTree] -> ShowS
showList :: [WaveletTree] -> ShowS
Show)

-- | Coarsest approximation coefficients after all requested levels.
waveletApproximation :: WaveletTree -> [Double]
waveletApproximation :: WaveletTree -> [Double]
waveletApproximation (WaveletTree Wavelet
_ Int
_ [Double]
a [[Double]]
_) = [Double]
a

-- | Detail bands ordered from finest first to coarsest last.
waveletDetails :: WaveletTree -> [[Double]]
waveletDetails :: WaveletTree -> [[Double]]
waveletDetails (WaveletTree Wavelet
_ Int
_ [Double]
_ [[Double]]
d) = [[Double]]
d

waveletLow :: Wavelet -> [Double]
waveletLow :: Wavelet -> [Double]
waveletLow Wavelet
Haar = [Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sqrt Double
2, Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sqrt Double
2]
waveletLow Wavelet
Daubechies2 =
    [ (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double
forall a. Floating a => a -> a
sqrt Double
3) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sqrt Double
2)
    , (Double
3 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double
forall a. Floating a => a -> a
sqrt Double
3) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sqrt Double
2)
    , (Double
3 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double
forall a. Floating a => a -> a
sqrt Double
3) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sqrt Double
2)
    , (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double
forall a. Floating a => a -> a
sqrt Double
3) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sqrt Double
2)
    ]
waveletLow Wavelet
Daubechies4 =
    [ Double
0.2303778133088965
    , Double
0.7148465705529154
    , Double
0.6308807679298587
    , -Double
0.02798376941685985
    , -Double
0.1870348117190931
    , Double
0.03084138183556076
    , Double
0.0328830116668852
    , -Double
0.01059740178506903
    ]
waveletLow Wavelet
Biorthogonal53 = []
waveletStep :: Wavelet -> [Double] -> ([Double], [Double])
waveletStep :: Wavelet -> [Double] -> ([Double], [Double])
waveletStep Wavelet
Biorthogonal53 [Double]
xs =
    let n :: Int
n = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2
        e :: [Double]
e = [[Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
i) | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
        o :: [Double]
o = [[Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
        d :: [Double]
d = [[Double]
o [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
- ([Double]
e [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
+ [Double]
e [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! ((Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`mod` Int
n)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2 | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
        a :: [Double]
a = [[Double]
e [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
+ ([Double]
d [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! ((Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`mod` Int
n) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ [Double]
d [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
4 | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
     in ([Double]
a, [Double]
d)
waveletStep Wavelet
family [Double]
xs =
    let n :: Int
n = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs
        lo :: [Double]
lo = Wavelet -> [Double]
waveletLow Wavelet
family
        hi :: [Double]
hi = (Int -> Double -> Double) -> [Int] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Int
j Double
h -> if Int -> Bool
forall a. Integral a => a -> Bool
even Int
j then Double
h else -Double
h) [Int
0 :: Int ..] ([Double] -> [Double]
forall a. [a] -> [a]
reverse [Double]
lo)
        apply :: [Double] -> [Double]
apply [Double]
kernel = [[Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double
h Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double]
xs [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! ((Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
j) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`mod` Int
n) | (Int
j, Double
h) <- [Int] -> [Double] -> [(Int, Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] [Double]
kernel] | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
     in ([Double] -> [Double]
apply [Double]
lo, [Double] -> [Double]
apply [Double]
hi)
inverseStep :: Wavelet -> [Double] -> [Double] -> [Double]
inverseStep :: Wavelet -> [Double] -> [Double] -> [Double]
inverseStep Wavelet
Biorthogonal53 [Double]
a [Double]
d =
    let n :: Int
n = [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
a
        e :: [Double]
e = [[Double]
a [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
- ([Double]
d [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! ((Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`mod` Int
n) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ [Double]
d [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
4 | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
        o :: [Double]
o = [[Double]
d [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
+ ([Double]
e [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
+ [Double]
e [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! ((Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`mod` Int
n)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2 | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
     in [[Double]] -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat ((Double -> Double -> [Double])
-> [Double] -> [Double] -> [[Double]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Double
x Double
y -> [Double
x, Double
y]) [Double]
e [Double]
o)
inverseStep Wavelet
family [Double]
a [Double]
d =
    let n :: Int
n = Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
a
        lo :: [Double]
lo = Wavelet -> [Double]
waveletLow Wavelet
family
        hi :: [Double]
hi = (Int -> Double -> Double) -> [Int] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Int
j Double
h -> if Int -> Bool
forall a. Integral a => a -> Bool
even Int
j then Double
h else -Double
h) [Int
0 :: Int ..] ([Double] -> [Double]
forall a. [a] -> [a]
reverse [Double]
lo)
     in [ [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum
            [ [Double]
a [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
l Double -> Double -> Double
forall a. Num a => a -> a -> a
+ [Double]
d [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
h
            | (Int
j, (Double
l, Double
h)) <- [Int] -> [(Double, Double)] -> [(Int, (Double, Double))]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] ([Double] -> [Double] -> [(Double, Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Double]
lo [Double]
hi)
            , let residue :: Int
residue = (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
j) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`mod` Int
n
            , Int -> Bool
forall a. Integral a => a -> Bool
even Int
residue
            , let i :: Int
i = Int
residue Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2
            ]
        | Int
k <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
        ]

-- | Periodic multilevel analysis; length must be divisible by 2^levels, with levels 0..20 and length <=1048576.
dwt :: Wavelet -> Int -> [Double] -> Either MathError WaveletTree
dwt :: Wavelet -> Int -> [Double] -> Either MathError WaveletTree
dwt Wavelet
family Int
levels [Double]
xs = do
    Bool -> String -> Either MathError ()
require
        ( Int
levels Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0
            Bool -> Bool -> Bool
&& Int
levels Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
20
            Bool -> Bool -> Bool
&& Bool -> Bool
not ([Double] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Double]
xs)
            Bool -> Bool -> Bool
&& [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
1048576
            Bool -> Bool -> Bool
&& Int -> Integer
forall a. Integral a => a -> Integer
toInteger ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`mod` (Integer
2 Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
levels :: Integer) Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
0
        )
        String
"DWT length must be positive and divisible by 2^levels"
    _ <- [Double] -> Either MathError [Double]
vector [Double]
xs
    let go t
0 [Double]
a [[Double]]
details = WaveletTree -> Either MathError WaveletTree
forall a b. b -> Either a b
Right (Wavelet -> Int -> [Double] -> [[Double]] -> WaveletTree
WaveletTree Wavelet
family ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs) [Double]
a ([[Double]] -> [[Double]]
forall a. [a] -> [a]
reverse [[Double]]
details))
        go t
k [Double]
a [[Double]]
details = let ([Double]
next, [Double]
d) = Wavelet -> [Double] -> ([Double], [Double])
waveletStep Wavelet
family [Double]
a in [Double] -> Either MathError [Double]
vector ([Double]
next [Double] -> [Double] -> [Double]
forall a. [a] -> [a] -> [a]
++ [Double]
d) Either MathError [Double]
-> Either MathError WaveletTree -> Either MathError WaveletTree
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> t -> [Double] -> [[Double]] -> Either MathError WaveletTree
go (t
k t -> t -> t
forall a. Num a => a -> a -> a
- t
1) [Double]
next ([Double]
d [Double] -> [[Double]] -> [[Double]]
forall a. a -> [a] -> [a]
: [[Double]]
details)
    go levels xs []

-- | Invert multilevel analysis using orthogonal synthesis or inverse CDF 5/3 lifting.
idwt :: WaveletTree -> Either MathError [Double]
idwt :: WaveletTree -> Either MathError [Double]
idwt (WaveletTree Wavelet
family Int
_ [Double]
approximation [[Double]]
details) = ([Double] -> [Double] -> Either MathError [Double])
-> [Double] -> [[Double]] -> Either MathError [Double]
forall (t :: * -> *) (m :: * -> *) b a.
(Foldable t, Monad m) =>
(b -> a -> m b) -> b -> t a -> m b
foldM (\[Double]
a [Double]
d -> [Double] -> Either MathError [Double]
vector (Wavelet -> [Double] -> [Double] -> [Double]
inverseStep Wavelet
family [Double]
a [Double]
d)) [Double]
approximation ([[Double]] -> [[Double]]
forall a. [a] -> [a]
reverse [[Double]]
details)

-- | Admissible corrected complex Morlet, real Mexican hat, or real first Gaussian derivative.
data MotherWavelet = Morlet Double | MexicanHat | GaussianDerivative deriving (MotherWavelet -> MotherWavelet -> Bool
(MotherWavelet -> MotherWavelet -> Bool)
-> (MotherWavelet -> MotherWavelet -> Bool) -> Eq MotherWavelet
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: MotherWavelet -> MotherWavelet -> Bool
== :: MotherWavelet -> MotherWavelet -> Bool
$c/= :: MotherWavelet -> MotherWavelet -> Bool
/= :: MotherWavelet -> MotherWavelet -> Bool
Eq, Int -> MotherWavelet -> ShowS
[MotherWavelet] -> ShowS
MotherWavelet -> String
(Int -> MotherWavelet -> ShowS)
-> (MotherWavelet -> String)
-> ([MotherWavelet] -> ShowS)
-> Show MotherWavelet
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> MotherWavelet -> ShowS
showsPrec :: Int -> MotherWavelet -> ShowS
$cshow :: MotherWavelet -> String
show :: MotherWavelet -> String
$cshowList :: [MotherWavelet] -> ShowS
showList :: [MotherWavelet] -> ShowS
Show)

{- | Direct zero-extended CWT. Scales are in seconds; coefficients use dt/sqrt(s)
quadrature. Morlet subtracts its nonzero-mean correction exp(-omega²/2).
-}
cwt :: MotherWavelet -> Double -> [Double] -> [Double] -> Either MathError [[Complex Double]]
cwt :: MotherWavelet
-> Double
-> [Double]
-> [Double]
-> Either MathError [[Complex Double]]
cwt MotherWavelet
mother Double
spacing [Double]
scales [Double]
xs = do
    Bool -> String -> Either MathError ()
require
        ( Double -> Bool
finite Double
spacing
            Bool -> Bool -> Bool
&& Double
spacing Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0
            Bool -> Bool -> Bool
&& Bool -> Bool
not ([Double] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Double]
scales)
            Bool -> Bool -> Bool
&& (Double -> Bool) -> [Double] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (\Double
s -> Double -> Bool
finite Double
s Bool -> Bool -> Bool
&& Double
s Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0) [Double]
scales
            Bool -> Bool -> Bool
&& Int -> Integer
forall a. Integral a => a -> Integer
toInteger ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
scales) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
forall a. Integral a => a -> Integer
toInteger ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs) Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int) Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
20000000
        )
        String
"invalid CWT spacing/scales/work budget"
    _ <- [Double] -> Either MathError [Double]
vector [Double]
xs
    case mother of
        Morlet Double
omega -> Bool -> String -> Either MathError ()
require (Double -> Bool
finite Double
omega Bool -> Bool -> Bool
&& Double
omega Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& Double
omega Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
50) String
"Morlet omega must be 0..50"
        MotherWavelet
_ -> () -> Either MathError ()
forall a b. b -> Either a b
Right ()
    let psi Double
t = case MotherWavelet
mother of
            Morlet Double
omega -> ((Double -> Double
forall a. Floating a => a -> a
exp (-Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
t Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (-Double
0.25)) Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
0) Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
* (Double -> Complex Double
forall a. Floating a => a -> Complex a
cis (Double
omega Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
t) Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
- (Double -> Double
forall a. Floating a => a -> a
exp (-Double
omega Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
omega Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2) Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
0))
            MotherWavelet
MexicanHat -> (Double
2 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double -> Double
forall a. Floating a => a -> a
sqrt Double
3 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Floating a => a -> a -> a
** Double
0.25) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
t) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
exp (-Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
t Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2)) Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
0
            MotherWavelet
GaussianDerivative -> (-Double -> Double
forall a. Floating a => a -> a
sqrt Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (-Double
0.25) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
exp (-Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
t Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2)) Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
0
        rows =
            [ [ [Complex Double] -> Complex Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum
                    [ (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
spacing Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sqrt Double
scale Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
0) Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
* Complex Double -> Complex Double
forall a. Num a => Complex a -> Complex a
conjugate (Double -> Complex Double
psi (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
b) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
spacing Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
scale))
                    | (Int
j, Double
x) <- [Int] -> [Double] -> [(Int, Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] [Double]
xs
                    ]
              | Int
b <- [Int
0 .. [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
xs Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
              ]
            | Double
scale <- [Double]
scales
            ]
    _ <- vector (concatMap (\[Complex Double]
row -> (Complex Double -> [Double]) -> [Complex Double] -> [Double]
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap (\Complex Double
z -> [Complex Double -> Double
forall a. Complex a -> a
realPart Complex Double
z, Complex Double -> Double
forall a. Complex a -> a
imagPart Complex Double
z]) [Complex Double]
row) rows)
    Right rows

-- | Opaque complete binary low/high analysis tree in natural binary order, rather than frequency-sorted Gray order.
data WaveletPacket = PacketLeaf [Double] | PacketNode Wavelet WaveletPacket WaveletPacket
    deriving (WaveletPacket -> WaveletPacket -> Bool
(WaveletPacket -> WaveletPacket -> Bool)
-> (WaveletPacket -> WaveletPacket -> Bool) -> Eq WaveletPacket
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: WaveletPacket -> WaveletPacket -> Bool
== :: WaveletPacket -> WaveletPacket -> Bool
$c/= :: WaveletPacket -> WaveletPacket -> Bool
/= :: WaveletPacket -> WaveletPacket -> Bool
Eq, Int -> WaveletPacket -> ShowS
[WaveletPacket] -> ShowS
WaveletPacket -> String
(Int -> WaveletPacket -> ShowS)
-> (WaveletPacket -> String)
-> ([WaveletPacket] -> ShowS)
-> Show WaveletPacket
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> WaveletPacket -> ShowS
showsPrec :: Int -> WaveletPacket -> ShowS
$cshow :: WaveletPacket -> String
show :: WaveletPacket -> String
$cshowList :: [WaveletPacket] -> ShowS
showList :: [WaveletPacket] -> ShowS
Show)

-- | Recursively analyze both approximation and detail branches to the supplied depth.
waveletPackets :: Wavelet -> Int -> [Double] -> Either MathError WaveletPacket
waveletPackets :: Wavelet -> Int -> [Double] -> Either MathError WaveletPacket
waveletPackets Wavelet
family Int
levels [Double]
xs = do
    _ <- Wavelet -> Int -> [Double] -> Either MathError WaveletTree
dwt Wavelet
family Int
levels [Double]
xs
    let go t
0 [Double]
values = [Double] -> WaveletPacket
PacketLeaf [Double]
values
        go t
k [Double]
values = let ([Double]
a, [Double]
d) = Wavelet -> [Double] -> ([Double], [Double])
waveletStep Wavelet
family [Double]
values in Wavelet -> WaveletPacket -> WaveletPacket -> WaveletPacket
PacketNode Wavelet
family (t -> [Double] -> WaveletPacket
go (t
k t -> t -> t
forall a. Num a => a -> a -> a
- t
1) [Double]
a) (t -> [Double] -> WaveletPacket
go (t
k t -> t -> t
forall a. Num a => a -> a -> a
- t
1) [Double]
d)
    Right (go levels xs)

-- | Reconstruct a complete packet tree exactly up to floating roundoff.
inverseWaveletPackets :: WaveletPacket -> Either MathError [Double]
inverseWaveletPackets :: WaveletPacket -> Either MathError [Double]
inverseWaveletPackets (PacketLeaf [Double]
xs) = [Double] -> Either MathError [Double]
vector [Double]
xs
inverseWaveletPackets (PacketNode Wavelet
family WaveletPacket
a WaveletPacket
d) = do aa <- WaveletPacket -> Either MathError [Double]
inverseWaveletPackets WaveletPacket
a; dd <- inverseWaveletPackets d; vector (inverseStep family aa dd)

-- | Leaf coefficient bands in natural binary low/high tree order.
packetBands :: WaveletPacket -> [[Double]]
packetBands :: WaveletPacket -> [[Double]]
packetBands (PacketLeaf [Double]
xs) = [[Double]
xs]
packetBands (PacketNode Wavelet
_ WaveletPacket
a WaveletPacket
d) = WaveletPacket -> [[Double]]
packetBands WaveletPacket
a [[Double]] -> [[Double]] -> [[Double]]
forall a. [a] -> [a] -> [a]
++ WaveletPacket -> [[Double]]
packetBands WaveletPacket
d

-- | Soft-threshold detail coefficients; optional explicit threshold, otherwise finest-band MAD universal threshold.
waveletDenoise :: Wavelet -> Int -> Maybe Double -> [Double] -> Either MathError [Double]
waveletDenoise :: Wavelet
-> Int -> Maybe Double -> [Double] -> Either MathError [Double]
waveletDenoise Wavelet
family Int
levels Maybe Double
threshold [Double]
xs = do
    tree@(WaveletTree _ n a details) <- Wavelet -> Int -> [Double] -> Either MathError WaveletTree
dwt Wavelet
family Int
levels [Double]
xs
    let sigma =
            if [[Double]] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [[Double]]
details
                then Double
0
                else
                    let finest :: [Double]
finest = (Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Double -> Double
forall a. Num a => a -> a
abs ([[Double]]
details [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
0) in [Double] -> [Double]
forall a. Ord a => [a] -> [a]
sort [Double]
finest [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
finest Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
0.6744897501960817
        t = Double -> (Double -> Double) -> Maybe Double -> Double
forall b a. b -> (a -> b) -> Maybe a -> b
maybe (Double
sigma Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sqrt (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
log (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n))) Double -> Double
forall a. a -> a
id Maybe Double
threshold
    require (finite t && t >= 0) "invalid wavelet threshold"
    if null details
        then idwt tree
        else idwt (WaveletTree family n a [[signum x * max 0 (abs x - t) | x <- d] | d <- details])