-----------------------------------------------------------------------------
-- |
-- Module  :  ForSyDe.Shallow.Utility.FilterLib
-- Copyright   :  (c) ForSyDe Group, KTH 2007-2008
-- License     :  BSD-style (see the file LICENSE)
-- 
-- Maintainer  :  forsyde-dev@ict.kth.se
-- Stability   :  experimental
-- Portability :  portable
--
-- This is the filter library for ForSyDe heterogeneous MoCs - CT-MoC, SR-MoC,
-- and Untimed-MoC.
--
-- The filters at CT-MoC are based on the filters implemented at SR-MoC and Untimed-MoC, 
-- which means a signal in CT-MoC is always digitalized by a A\/D converter, processed by 
-- digital filters at SR or Untimed domain, and converted back into a CT domain signal by 
-- a D\/A converter. A CT-filter is composed of one A\/D converter, one digital filter, and 
-- one D\/A converter.
--
-- The implementation of the filters at untimed and synchronous domains follows the
-- way in a paper /An introduction to Haskell with applications to digital signal processing, 
-- David M. Goblirsch, in Proceedings of the 1994 ACM symposium on Applied computing./,
-- where the details of the FIR\/IIR, AR, and ARMA filters are introduced. The ARMA filter
-- is the kernel part of our general linear filter 'zLinearFilter' in Z-domain at SR\/Untimed
-- MoC, and is also the kernel digital filter for the linear filter 'sLinearFilter' in 
-- S-domain at CT-MoC.
-----------------------------------------------------------------------------
module ForSyDe.Shallow.Utility.FilterLib (
      -- *FIR filter
      firFilter,
      -- *AR and ARMA filter trim
      arFilterTrim, armaFilterTrim,
      -- *The solver mode
      SolverMode(..),
      -- *The general linear filter in S-domain
      sLinearFilter,
      -- *The general linear filter in Z-domain
      zLinearFilter,
      -- *s2z domain coefficient tranformation
      s2zCoef,
      -- *The Z-domain to ARMA coefficient tranformation
      h2ARMACoef
     )
    where 

import ForSyDe.Shallow.MoC
--import ForSyDe.Shallow.MoC.CT
import ForSyDe.Shallow.Core
import ForSyDe.Shallow.Utility.PolyArith
import Data.List (zipWith5)

-- |The FIR filter. Let '[x_n]' denote the input signal, '[y_n]' denote the ouput
-- signal, and '[h_n]' the impulse response of the filter. Suppose the length of
-- the impulse responses is M samples. The formula for '[y_n]' is 
-- $sum_{k=0}^{M-1} h_k*x_{n-k}$.
firFilter :: (Num a) => [a]  -- ^Coefficients of the FIR filter
         -> Signal a -- ^Input signal
         -> Signal a -- ^Output signal
firFilter :: [a] -> Signal a -> Signal a
firFilter [a]
hs Signal a
xs = ([a] -> a -> [a]) -> ([a] -> a -> a) -> [a] -> Signal a -> Signal a
forall a b c.
(a -> b -> a) -> (a -> b -> c) -> a -> Signal b -> Signal c
mealySY [a] -> a -> [a]
forall a. [a] -> a -> [a]
stateF ([a] -> [a] -> a -> a
forall b. Num b => [b] -> [b] -> b -> b
outF [a]
hs) (Int -> a -> [a]
forall a. Int -> a -> [a]
repeatN ([a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
hs) a
0) Signal a
xs
  where
    stateF :: [a] -> a -> [a]
stateF [a]
xs0 a
x = [a] -> a -> [a]
forall a. [a] -> a -> [a]
fixedList [a]
xs0 a
x
    outF :: [b] -> [b] -> b -> b
outF [b]
hs [b]
xs0 b
x = [b] -> [b] -> b
forall b. Num b => [b] -> [b] -> b
iprod [b]
hs ([b] -> b) -> [b] -> b
forall a b. (a -> b) -> a -> b
$ [b] -> b -> [b]
forall a. [a] -> a -> [a]
fixedList [b]
xs0 b
x

-- |The autoregressive filter is the simplest IIR filter. It is characterized 
-- by a list of numbers '[a_1,a_2,...,a_M]' called the autoregression 
-- coefficients and a single number 'b' called the gain. 'M' is the order of 
-- the filter. Let '[x_n]' denote the input signal, '[y_n]' denote the ouput
-- signal. The formula for '[y_n]' is  $\sum_{k=1}^M {a_k*y_{n-k}+b*x_n}$. 
-- Although it is an IIR filter, here we only get the same length of ouput 
-- signal as the input signal.
arFilterTrim :: (Num a, Fractional a) => 
           [a]   -- ^Coefficients of the AR filter.
         -> a    -- ^The gain
         -> Signal a -- ^Input signal
         -> Signal a -- ^Output signal
arFilterTrim :: [a] -> a -> Signal a -> Signal a
arFilterTrim [a]
as a
b Signal a
xs = 
    ([a] -> a -> [a]) -> ([a] -> a -> a) -> [a] -> Signal a -> Signal a
forall a b c.
(a -> b -> a) -> (a -> b -> c) -> a -> Signal b -> Signal c
mealySY ([a] -> a -> [a] -> a -> [a]
forall a. Num a => [a] -> a -> [a] -> a -> [a]
stateF [a]
as a
b) ([a] -> a -> [a] -> a -> a
forall a. Num a => [a] -> a -> [a] -> a -> a
outF [a]
as a
b) (Int -> a -> [a]
forall a. Int -> a -> [a]
repeatN ([a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
as) a
0) Signal a
xs
  where
    stateF :: [a] -> a -> [a] -> a -> [a]
stateF [a]
as a
b [a]
xs0 a
x = [a] -> a -> [a]
forall a. [a] -> a -> [a]
fixedList [a]
xs0 (a -> [a]) -> a -> [a]
forall a b. (a -> b) -> a -> b
$ [a] -> a -> [a] -> a -> a
forall a. Num a => [a] -> a -> [a] -> a -> a
outF [a]
as a
b [a]
xs0 a
x 
    outF :: [a] -> a -> [a] -> a -> a
outF [a]
as a
b [a]
xs0 a
x = a
ba -> a -> a
forall a. Num a => a -> a -> a
*a
x a -> a -> a
forall a. Num a => a -> a -> a
+ ([a] -> [a] -> a
forall b. Num b => [b] -> [b] -> b
iprod [a]
as [a]
xs0)

-- |The ARMA filter combines the FIR and AR filters. Let '[x_n]' denote the 
-- input signal, '[y_n]' denote the ouput signal. The formula for '[y_n]' is
--  $\sum_{k=1}^M {a_k*y_{n-k}+b*x_n} + sum_{i=0}^{N-1} b_i*x_{n-i}$. The ARMA
-- filter can be defined as the composition of an FIR filter having the impulse
-- reponse '[b_0,b_1,...,b_N-1]' and an AR filter having the regression 
-- coefficients '[a_1,a_2,...,a_M]' and a gain of '1'. Although it is an IIR 
-- filter, here we only get the same length of ouput signal as the input signal.
armaFilterTrim :: (Num a, Fractional a) => 
             [a]  -- ^Coefficients of the FIR filter
          -> [a]  -- ^Coefficients of the AR filter.
          -> Signal a -- ^Input signal
          -> Signal a -- ^Output signal
armaFilterTrim :: [a] -> [a] -> Signal a -> Signal a
armaFilterTrim [a]
bs [a]
as = [a] -> a -> Signal a -> Signal a
forall a. (Num a, Fractional a) => [a] -> a -> Signal a -> Signal a
arFilterTrim [a]
as a
1 (Signal a -> Signal a)
-> (Signal a -> Signal a) -> Signal a -> Signal a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [a] -> Signal a -> Signal a
forall a. Num a => [a] -> Signal a -> Signal a
firFilter [a]
bs


-- |The solver mode.
data SolverMode = S2Z   -- ^Tustin tranfer from s-domain to z-domain
        | RK4   -- ^Runge Kutta 4 with fixed simulation steps
  deriving (Int -> SolverMode -> ShowS
[SolverMode] -> ShowS
SolverMode -> String
(Int -> SolverMode -> ShowS)
-> (SolverMode -> String)
-> ([SolverMode] -> ShowS)
-> Show SolverMode
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [SolverMode] -> ShowS
$cshowList :: [SolverMode] -> ShowS
show :: SolverMode -> String
$cshow :: SolverMode -> String
showsPrec :: Int -> SolverMode -> ShowS
$cshowsPrec :: Int -> SolverMode -> ShowS
Show, SolverMode -> SolverMode -> Bool
(SolverMode -> SolverMode -> Bool)
-> (SolverMode -> SolverMode -> Bool) -> Eq SolverMode
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: SolverMode -> SolverMode -> Bool
$c/= :: SolverMode -> SolverMode -> Bool
== :: SolverMode -> SolverMode -> Bool
$c== :: SolverMode -> SolverMode -> Bool
Eq)

-- |The general linear filter in S-domain at CT-MoC. As the kernel
-- implementation is in Z-domain, the smapling rate should be specified. 
-- It is used on the S-transformation with the following forms, with 
-- coefficients for the descending powers of 's' and m < n.
--
-- >    b_0*s^m + b_1*s^m-1 + ... + b_m-1*s^1 + b_m*s^0
-- >H(s) = ------------------------------------------------     (Eq 1)
-- >    a_0*s^n + a_1*s^n-1 + ... + a_n-1*s^1 + a_n*s^0
--
-- If we multiply both the numerator and the denominator with s^-n, we get 
-- another equivelent canonical form
--
-- >    b_0*s^m-n + b_1*s^m-n-1 + ... + b_m-1*s^1-n + b_m*s^-n
-- >H(s) = -----------------------------------------------------    (Eq 2)
-- >    a_0*s^0 + a_1*s^-1 + ... + a_n-1*s^1-n + a_n*s^-n
--
-- To instantiate a filter, with sampling interval 'T ', we use
--
-- > filter1 = sLinearFilter T [1] [2,1]
-- 
-- to get a filter  with the transfer function
-- 
-- >      1
-- >H(s) = --------
-- >   2*s + 1
-- 
-- and
--
-- > filter2 = sLinearFilter T [2,1] [1,2,2]
--
-- to get another filter with the transfer function
-- 
-- >       2*s +1
-- >H(s) = ----------------
-- >    s^2 + 2*s + 2
--
-- There are two solver modes, 'S2Z' and 'RK4'.
-- Caused by the precision problem, the time interval in CT uses Rational data
-- type and the digital data types in all the domains are Double.
sLinearFilter :: (Num a, Fractional a, Show a, Eq a) =>
    SolverMode   -- ^Specify the solver mode
    ->  Rational     -- ^Fixed simulation interval
    -> [a]  -- ^Coefficients of the polynomial numerator in s-domain
    -> [a]  -- ^Coefficients of the polynomial denominator in s-domain
    -> Signal (SubsigCT a)-- ^Input CT-signal
    -> Signal (SubsigCT a)-- ^Output CT-signal
sLinearFilter :: SolverMode
-> Rational
-> [a]
-> [a]
-> Signal (SubsigCT a)
-> Signal (SubsigCT a)
sLinearFilter SolverMode
filterMode Rational
step [a]
bs [a]
as Signal (SubsigCT a)
inS =  Signal (SubsigCT a)
outS 
  where
    -- A2D conversion
    inSDigital :: Signal a
inSDigital = Rational -> Signal (SubsigCT a) -> Signal a
forall a.
(Num a, Show a) =>
Rational -> Signal (SubsigCT a) -> Signal a
a2dConverter Rational
step Signal (SubsigCT a)
inS
    -- D2A conversion
    outS :: Signal (SubsigCT a)
outS =  DACMode -> Rational -> Signal a -> Signal (SubsigCT a)
forall a.
(Fractional a, Show a) =>
DACMode -> Rational -> Signal a -> Signal (SubsigCT a)
d2aConverter DACMode
DAhold Rational
step Signal a
outSDigital
    -- Digital filter
    outSDigital :: Signal a
outSDigital | SolverMode
filterMode SolverMode -> SolverMode -> Bool
forall a. Eq a => a -> a -> Bool
== SolverMode
S2Z = [a] -> [a] -> Signal a -> Signal a
forall a.
(Num a, Fractional a) =>
[a] -> [a] -> Signal a -> Signal a
armaFilterTrim [a]
bs' [a]
as' Signal a
inSDigital
        | Bool
otherwise =  Rational -> [a] -> [a] -> Signal a -> Signal a
forall a.
(Fractional a, Show a, Eq a) =>
Rational -> [a] -> [a] -> Signal a -> Signal a
rk4FilterDigital Rational
step [a]
as [a]
bs Signal a
inSDigital
          where ([a]
bs',[a]
as') = ([a], [a]) -> ([a], [a])
forall a. (Num a, Fractional a) => ([a], [a]) -> ([a], [a])
h2ARMACoef (([a], [a]) -> ([a], [a])) -> ([a], [a]) -> ([a], [a])
forall a b. (a -> b) -> a -> b
$ Rational -> [a] -> [a] -> ([a], [a])
forall a.
(Num a, Fractional a, Eq a) =>
Rational -> [a] -> [a] -> ([a], [a])
s2zCoef Rational
step [a]
bs [a]
as

-- |Digital filter using Runge Kutta 4 solver.
rk4FilterDigital :: (Fractional a, Show a, Eq a) => 
        Rational -> [a] -> [a] -> Signal a -> Signal a
rk4FilterDigital :: Rational -> [a] -> [a] -> Signal a -> Signal a
rk4FilterDigital Rational
step [a]
as [a]
bs Signal a
inSDigital = Signal a
outSDigital
  where
    -- Below are the skeletons of the RK-4 solver, with
    -- input signal 'inSDigital' and output signal 'outSDigital'
    -- Coefficients handling
    as'' :: [a]
as'' = (a -> Bool) -> [a] -> [a]
forall a. (a -> Bool) -> [a] -> [a]
dropWhile (\a
x -> a
xa -> a -> Bool
forall a. Eq a => a -> a -> Bool
==a
0.0) [a]
as
    a0 :: a
a0 = [a] -> a
forall a. [a] -> a
head [a]
as''
    -- Normalized the coefficients
    as' :: [a]
as' = [a] -> [a]
forall a. [a] -> [a]
reverse ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$ [a] -> [a]
forall a. [a] -> [a]
tail ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$ (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (\a
x -> -a
xa -> a -> a
forall a. Fractional a => a -> a -> a
/a
a0) [a]
as''
    bs' :: [a]
bs' = [a] -> [a]
forall a. [a] -> [a]
reverse ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$ (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (\a
x -> a
xa -> a -> a
forall a. Fractional a => a -> a -> a
/a
a0) [a]
bs
    -- Order of the filter
    orderFilter :: Int
orderFilter = [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
as'
    -- The last state function, '0' is for the time 
    fXn :: [a] -> a
fXn = [a] -> [a] -> a
forall b. Num b => [b] -> [b] -> b
iprod (a
0a -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a]
as')
    -- The functions for the observalbe state matrix 'A'
    stateFunctions :: [[a] -> a]
stateFunctions = Int -> [[a] -> a]
forall t t1. (Num t, Num t1, Eq t) => t -> [[t1] -> t1]
ffn' Int
orderFilter [[a] -> a] -> [[a] -> a] -> [[a] -> a]
forall a. [a] -> [a] -> [a]
++ [[a] -> a
fXn]
    -- Initial states
    initialStates :: [a]
initialStates = Int -> a -> [a]
forall a. Int -> a -> [a]
repeatN Int
orderFilter a
0.0
    inputSteps :: Signal a
inputSteps = [a] -> Signal a
forall a. [a] -> Signal a
signal ([a] -> Signal a) -> [a] -> Signal a
forall a b. (a -> b) -> a -> b
$ a -> [a]
forall a. a -> [a]
repeat a
step'
    -- The states signal
    statesSignal :: Signal [a]
statesSignal = a -> [a] -> [[a] -> a] -> Signal a -> Signal a -> Signal [a]
forall a.
(Num a, Fractional a) =>
a -> [a] -> [[a] -> a] -> Signal a -> Signal a -> Signal [a]
rks4InSY a
0.0 [a]
initialStates [[a] -> a]
stateFunctions 
             Signal a
inputSteps Signal a
inSDigital --xs
    -- The ouput digital signal 
    outSDigital :: Signal a
outSDigital = ([a] -> a) -> Signal [a] -> Signal a
forall a b. (a -> b) -> Signal a -> Signal b
mapSY ([a] -> [a] -> a
forall b. Num b => [b] -> [b] -> b
iprod [a]
bs') Signal [a]
statesSignal
    -- The fixed simulation step
    step' :: a
step' = Rational -> a
forall a. Fractional a => Rational -> a
fromRational Rational
step

-- The length of the function list is 'n-1' for nth order filter
ffn' :: (Num t, Num t1, Eq t) => t -> [[t1] -> t1]
ffn' :: t -> [[t1] -> t1]
ffn' t
n = Int -> t -> [[t1] -> t1]
forall t1 t. (Num t1, Num t, Eq t) => Int -> t -> [[t1] -> t1]
ffn Int
0 t
n

-- Construct the functions for the diagonal '1'
ffn :: (Num t1, Num t, Eq t) => Int -> t -> [[t1] -> t1]
ffn :: Int -> t -> [[t1] -> t1]
ffn Int
_ t
1 = []
ffn Int
m t
n = Int -> [t1] -> t1
forall t. Num t => Int -> [t] -> t
ff1 Int
m ([t1] -> t1) -> [[t1] -> t1] -> [[t1] -> t1]
forall a. a -> [a] -> [a]
: Int -> t -> [[t1] -> t1]
forall t1 t. (Num t1, Num t, Eq t) => Int -> t -> [[t1] -> t1]
ffn (Int
mInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) (t
nt -> t -> t
forall a. Num a => a -> a -> a
-t
1) 

ff1 :: Num t => Int -> [t] -> t
ff1 :: Int -> [t] -> t
ff1 Int
m = [t] -> [t] -> t
forall b. Num b => [b] -> [b] -> b
iprod ([t
0,t
0] [t] -> [t] -> [t]
forall a. [a] -> [a] -> [a]
++ (Int -> t -> [t]
forall a. Int -> a -> [a]
repeatN Int
m t
0) [t] -> [t] -> [t]
forall a. [a] -> [a] -> [a]
++ [t
1] [t] -> [t] -> [t]
forall a. [a] -> [a] -> [a]
++ (t -> [t]
forall a. a -> [a]
repeat t
0) )

-- |RK-4 to solve the 1st-order ODEs, with input signal.
rks4InSY :: (Num a, Fractional a) =>
      a     -- ^The initial time
  -> [a]       -- ^The initial state values
  -> [([a] -> a)] -- ^List of the functions of the ODEs.
  ->  Signal  a  -- ^Input signal of steps
  ->  Signal  a  -- ^Input signal
  ->  Signal [a] -- ^Next state signal
rks4InSY :: a -> [a] -> [[a] -> a] -> Signal a -> Signal a -> Signal [a]
rks4InSY a
x0 [a]
ys0 [[a] -> a]
fFs Signal a
hs Signal a
us = ([a] -> a -> a -> a -> [a])
-> [a] -> Signal a -> Signal a -> Signal a -> Signal [a]
forall a b c d.
(a -> b -> c -> d -> a)
-> a -> Signal b -> Signal c -> Signal d -> Signal a
scanl3SY [a] -> a -> a -> a -> [a]
stateF [a]
ys0 Signal a
xs Signal a
hs Signal a
us
  where
    stateF :: [a] -> a -> a -> a -> [a]
stateF [a]
ysn a
xn a
h a
ut = (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Num a => a -> a -> a
(+) (Int -> a -> [a]
forall a. Int -> a -> [a]
repeatN Int
orderODE' a
0.0 [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [a
uta -> a -> a
forall a. Num a => a -> a -> a
*a
h]) ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$ --Input value
              a -> a -> [[a] -> a] -> [a] -> [a]
forall a.
(Num a, Fractional a) =>
a -> a -> [[a] -> a] -> [a] -> [a]
rks4 a
h a
xn [[a] -> a]
fFs [a]
ysn 
    xs :: Signal a
xs = (a -> a -> a) -> a -> Signal a -> Signal a
forall a b. (a -> b -> a) -> a -> Signal b -> Signal a
scanldSY a -> a -> a
forall a. Num a => a -> a -> a
(+) a
x0 Signal a
hs
    -- Order -1 of the ODEs
    orderODE' :: Int
orderODE' = [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
ys0 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1

-- |One step RK-4 for the 1st-order ordinary differential equations (ODEs).
rks4 ::  (Num a, Fractional a) =>
     a    -- ^The step
     ->  a    -- ^Initial value of time
     -> [[a] -> a] -- ^List of the funcitons of the ODEs.
     -> [a]   -- ^List of the value at the current state
     -> [a]   -- ^List of the value at the next state
rks4 :: a -> a -> [[a] -> a] -> [a] -> [a]
rks4 a
h a
x0 [[a] -> a]
fFs [a]
ys0 = [a]
ys1
  where
    h_2 :: a
h_2 = a
ha -> a -> a
forall a. Fractional a => a -> a -> a
/a
2.0
    ks1 :: [a]
ks1 = (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a
ha -> a -> a
forall a. Num a => a -> a -> a
*) ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$ [a] -> [[a] -> a] -> [a]
forall a b. a -> [a -> b] -> [b]
map' (a
x0a -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a]
ys0) [[a] -> a]
fFs  --    -- map (h *) $ applyFt x0 fFs ys0
    ks2 :: [a]
ks2 = (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a
ha -> a -> a
forall a. Num a => a -> a -> a
*) ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$ [a] -> [[a] -> a] -> [a]
forall a b. a -> [a -> b] -> [b]
map' (a
x0a -> a -> a
forall a. Num a => a -> a -> a
+a
h_2a -> [a] -> [a]
forall a. a -> [a] -> [a]
:(a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\a
y a
k-> a
ya -> a -> a
forall a. Num a => a -> a -> a
+a
ka -> a -> a
forall a. Fractional a => a -> a -> a
/a
2.0) [a]
ys0 [a]
ks1) [[a] -> a]
fFs 
    ks3 :: [a]
ks3 = (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a
ha -> a -> a
forall a. Num a => a -> a -> a
*) ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$ [a] -> [[a] -> a] -> [a]
forall a b. a -> [a -> b] -> [b]
map' (a
x0a -> a -> a
forall a. Num a => a -> a -> a
+a
h_2a -> [a] -> [a]
forall a. a -> [a] -> [a]
:(a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\a
y a
k-> a
ya -> a -> a
forall a. Num a => a -> a -> a
+a
ka -> a -> a
forall a. Fractional a => a -> a -> a
/a
2.0) [a]
ys0 [a]
ks2) [[a] -> a]
fFs 
    ks4 :: [a]
ks4 = (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a
ha -> a -> a
forall a. Num a => a -> a -> a
*) ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$ [a] -> [[a] -> a] -> [a]
forall a b. a -> [a -> b] -> [b]
map' (a
x0a -> a -> a
forall a. Num a => a -> a -> a
+a
ha -> [a] -> [a]
forall a. a -> [a] -> [a]
:(a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\a
y a
k-> a
ya -> a -> a
forall a. Num a => a -> a -> a
+a
k) [a]
ys0 [a]
ks3) [[a] -> a]
fFs 
    ys1 :: [a]
ys1 = (a -> a -> a -> a -> a -> a)
-> [a] -> [a] -> [a] -> [a] -> [a] -> [a]
forall a b c d e f.
(a -> b -> c -> d -> e -> f)
-> [a] -> [b] -> [c] -> [d] -> [e] -> [f]
zipWith5 (\a
y0 a
k1 a
k2 a
k3 a
k4 -> a
y0 a -> a -> a
forall a. Num a => a -> a -> a
+ a
k1a -> a -> a
forall a. Fractional a => a -> a -> a
/a
6 a -> a -> a
forall a. Num a => a -> a -> a
+ a
k2a -> a -> a
forall a. Fractional a => a -> a -> a
/a
3 a -> a -> a
forall a. Num a => a -> a -> a
+ a
k3a -> a -> a
forall a. Fractional a => a -> a -> a
/a
3 a -> a -> a
forall a. Num a => a -> a -> a
+ a
k4a -> a -> a
forall a. Fractional a => a -> a -> a
/a
6)
        [a]
ys0 [a]
ks1 [a]
ks2 [a]
ks3 [a]
ks4

-- |The general linear filter in Z-domain.
zLinearFilter :: Fractional a => [a] -> [a] -> Signal a -> Signal a
zLinearFilter :: [a] -> [a] -> Signal a -> Signal a
zLinearFilter [a]
bs [a]
as = [a] -> [a] -> Signal a -> Signal a
forall a.
(Num a, Fractional a) =>
[a] -> [a] -> Signal a -> Signal a
armaFilterTrim [a]
bs' [a]
as'
  where
    bs' :: [a]
bs' = (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map ((\a
x a
y-> a
ya -> a -> a
forall a. Fractional a => a -> a -> a
/a
x ) ([a] -> a
forall a. [a] -> a
head [a]
as)) [a]
bs
    as' :: [a]
as' = (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map ((\a
x a
y-> -a
ya -> a -> a
forall a. Fractional a => a -> a -> a
/a
x ) ([a] -> a
forall a. [a] -> a
head [a]
as)) ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$ [a] -> [a]
forall a. [a] -> [a]
tail [a]
as

-- |s2z domain coefficient tranformation with a specified sampling rate.
-- The Tustin transformation is used for the transfer, with
--
-- >  2(z - 1)  
-- > s = ----------                 (Eq 3)
-- >  T(z + 1)
--
-- in which, 'T' is the sampling interval.
s2zCoef :: (Num a, Fractional a, Eq a) =>
    Rational      -- ^ Sampling rate in Z-domain
    -> [a]        -- ^ Coefficients of the polynomial numerator in s-domain
    -> [a]        -- ^ Coefficients of the polynomial denominator in s-domain
    -> ([a], [a]) -- ^ Tuple of the numerator and denominator 
                  --   coefficients in Z-domain
s2zCoef :: Rational -> [a] -> [a] -> ([a], [a])
s2zCoef Rational
sampleT [a]
bs [a]
as = ([a] -> [a]
forall a. [a] -> [a]
reverse [a]
bs', [a] -> [a]
forall a. [a] -> [a]
reverse [a]
as')
  where
    ([a]
bs',[a]
as') = Poly a -> ([a], [a])
forall a. Num a => Poly a -> ([a], [a])
getCoef Poly a
hZ    
    bsInv :: [a]
bsInv = [a] -> [a]
forall a. [a] -> [a]
reverse [a]
bs
    asInv :: [a]
asInv = [a] -> [a]
forall a. [a] -> [a]
reverse [a]
as
    numerator' :: Poly a
numerator' = (Poly a -> (a, Poly a) -> Poly a)
-> Poly a -> [(a, Poly a)] -> Poly a
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl (\Poly a
x (a, Poly a)
y -> Poly a -> Poly a -> Poly a
forall a. (Num a, Eq a) => Poly a -> Poly a -> Poly a
addPoly Poly a
x (Poly a -> Poly a) -> Poly a -> Poly a
forall a b. (a -> b) -> a -> b
$ a -> Poly a -> Poly a
forall a. Num a => a -> Poly a -> Poly a
scalePoly ((a, Poly a) -> a
forall a b. (a, b) -> a
fst (a, Poly a)
y) ((a, Poly a) -> Poly a
forall a b. (a, b) -> b
snd (a, Poly a)
y)) 
                 ([a] -> Poly a
forall a. [a] -> Poly a
Poly [a
0]) ([(a, Poly a)] -> Poly a) -> [(a, Poly a)] -> Poly a
forall a b. (a -> b) -> a -> b
$ [a] -> [Poly a] -> [(a, Poly a)]
forall a b. [a] -> [b] -> [(a, b)]
zip [a]
bsInv [Poly a]
sList
    denominator' :: Poly a
denominator' = (Poly a -> (a, Poly a) -> Poly a)
-> Poly a -> [(a, Poly a)] -> Poly a
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl (\Poly a
x (a, Poly a)
y -> Poly a -> Poly a -> Poly a
forall a. (Num a, Eq a) => Poly a -> Poly a -> Poly a
addPoly Poly a
x (Poly a -> Poly a) -> Poly a -> Poly a
forall a b. (a -> b) -> a -> b
$ a -> Poly a -> Poly a
forall a. Num a => a -> Poly a -> Poly a
scalePoly ((a, Poly a) -> a
forall a b. (a, b) -> a
fst (a, Poly a)
y) ((a, Poly a) -> Poly a
forall a b. (a, b) -> b
snd (a, Poly a)
y)) 
                   ([a] -> Poly a
forall a. [a] -> Poly a
Poly [a
0]) ([(a, Poly a)] -> Poly a) -> [(a, Poly a)] -> Poly a
forall a b. (a -> b) -> a -> b
$ [a] -> [Poly a] -> [(a, Poly a)]
forall a b. [a] -> [b] -> [(a, b)]
zip [a]
asInv [Poly a]
sList
    hZ :: Poly a
hZ = Poly a -> Poly a -> Poly a
forall a. Num a => Poly a -> Poly a -> Poly a
divPoly Poly a
numerator' Poly a
denominator'
    -- Tustin transform
    s :: Poly a
s = (Poly a, Poly a) -> Poly a
forall a. (Poly a, Poly a) -> Poly a
PolyPair ([a] -> Poly a
forall a. [a] -> Poly a
Poly [-a
2,a
2],[a] -> Poly a
forall a. [a] -> Poly a
Poly [Rational -> a
forall a. Fractional a => Rational -> a
fromRational Rational
sampleT,Rational -> a
forall a. Fractional a => Rational -> a
fromRational Rational
sampleT])
    sList :: [Poly a]
sList = (Int -> Poly a) -> [Int] -> [Poly a]
forall a b. (a -> b) -> [a] -> [b]
map (Poly a -> Int -> Poly a
forall a. Num a => Poly a -> Int -> Poly a
powerPoly Poly a
s) [Int
0..]

-- |The Z-domain to ARMA coefficient tranformation. It is used on the 
-- Z-transfer function
--
-- >    b_0*z^m-n + b_1*z^m-n-1 + ... + b_m-1*z^1-n + b_m*z^-n
-- >H(z) = -----------------------------------------------------    (Eq 4)
-- >    a_0*z^0 + a_1*z^-1 + ... + a_n-1*z^1-n + a_n*z^-n
--
-- which is normalized as
--
-- >    b_0/a_0*z^m-n + b_1/a_0*z^m-n-1 + ... + b_m/a_0*z^-n
-- >H(z) = -------------------------------------------------------  (Eq 5)
-- >    1 + a_1/a_0*z^-1 + ... + a_n-1/a_0*z^1-n + a_n/a_0*z^-n
--
-- The implementation coudl be
--
-- >y(k) = b_0/a_0*x_k+m-n + b_1/a_0*x_k+m-n-1 + ... + b_m/a_0*x_k-n
-- >                        (Eq 6)
-- >           - a_1/a_0*y_k-1 - ... - a_n/a_0*y_k-n
--
-- Then, we could get the coefficients of the ARMA filter.
h2ARMACoef :: (Num a, Fractional a) =>
       ([a], [a]) -- ^Coefficients in Z-domain
    -> ([a], [a]) -- ^Coefficients of the ARMA filter
h2ARMACoef :: ([a], [a]) -> ([a], [a])
h2ARMACoef ([a]
bs,[a]
as) = (a -> [a] -> [a]
forall a. Num a => a -> [a] -> [a]
scalePolyCoef a
a0_1 [a]
bs, 
          a -> [a] -> [a]
forall a. Num a => a -> [a] -> [a]
scalePolyCoef (a
0a -> a -> a
forall a. Num a => a -> a -> a
-a
a0_1) ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$ [a] -> [a]
forall a. [a] -> [a]
tail [a]
as)
  where
    a0_1 :: a
a0_1 = a
1.0a -> a -> a
forall a. Fractional a => a -> a -> a
/ [a] -> a
forall a. [a] -> a
head [a]
as

-- Helper functions

map' :: a -> [a->b] -> [b]
map' :: a -> [a -> b] -> [b]
map' = ([a -> b] -> a -> [b]) -> a -> [a -> b] -> [b]
forall a b c. (a -> b -> c) -> b -> a -> c
flip (([a -> b] -> a -> [b]) -> a -> [a -> b] -> [b])
-> ([a -> b] -> a -> [b]) -> a -> [a -> b] -> [b]
forall a b. (a -> b) -> a -> b
$ [a -> b] -> a -> [b]
forall (t :: * -> *) (m :: * -> *) a.
(Traversable t, Monad m) =>
t (m a) -> m (t a)
sequence 


-- |Computes the inner product.
iprod :: Num b => [b] -> [b] -> b
iprod :: [b] -> [b] -> b
iprod [b]
xs [b]
ys = [b] -> b
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [b
xb -> b -> b
forall a. Num a => a -> a -> a
*b
y | (b
x, b
y) <- [b] -> [b] -> [(b, b)]
forall a b. [a] -> [b] -> [(a, b)]
zip [b]
xs [b]
ys]

-- |Repeat an element for a given times.
repeatN :: Int -> a -> [a]
repeatN :: Int -> a -> [a]
repeatN Int
n = Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take Int
n ([a] -> [a]) -> (a -> [a]) -> a -> [a]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. a -> [a]
forall a. a -> [a]
repeat

-- |Maintain a fixed length of list like Fifo, except the outputs are ignored.
fixedList :: [a] -> a -> [a]
fixedList :: [a] -> a -> [a]
fixedList [a]
xs a
y = Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take ([a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
xs) ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$ a
ya -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a]
xs