{-# LANGUAGE NoImplicitPrelude #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE FlexibleInstances #-}
{- |
Copyright   :  (c) Henning Thielemann 2008
License     :  GPL

Maintainer  :  synthesizer@henning-thielemann.de
Stability   :  provisional
Portability :  requires multi-parameter type classes

Butterworth lowpass and highpass
-}
module Synthesizer.Plain.Filter.Recursive.Butterworth (
   Parameter,
   causal,
   causalPole,
   highpassCausalPole, highpassPole,
   lowpassCausalPole,  lowpassPole,
   modifier,
   parameter,
   partialParameter,
   runPole,

   -- used in Dimensional.Causal.FilterParameter
   checkedHalf,
   -- used in LLVM.Filter.Butterworth
   partialRatio,
   makeSines,
   ) where

import Synthesizer.Plain.Filter.Recursive (Passband(Lowpass,Highpass), Pole(Pole))
import qualified Synthesizer.Plain.Filter.Recursive.SecondOrderCascade as Cascade
import qualified Synthesizer.Plain.Filter.Recursive.SecondOrder as Filt2
import qualified Synthesizer.Plain.Signal   as Sig
import qualified Synthesizer.Plain.Modifier as Modifier
import qualified Synthesizer.Causal.Process as Causal
import Control.Arrow ((>>>), )

import qualified Algebra.Module                as Module
import qualified Algebra.Transcendental        as Trans
import qualified Algebra.Ring                  as Ring

import qualified Data.StorableVector as SV
import Foreign.Storable (Storable)

import NumericPrelude.Numeric
import NumericPrelude.Base



sineList, _sineListSlow, sineListFast :: (Trans.C a) => a -> [a]
sineList :: forall a. C a => a -> [a]
sineList = a -> [a]
forall a. C a => a -> [a]
sineListFast

_sineListSlow :: forall a. C a => a -> [a]
_sineListSlow a
x =
   (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map a -> a
forall a. C a => a -> a
sin ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$ (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a
xa -> a -> a
forall a. C a => a -> a -> a
*) ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$ (a -> a) -> a -> [a]
forall a. (a -> a) -> a -> [a]
iterate (a
2a -> a -> a
forall a. C a => a -> a -> a
+) a
1

sineListFast :: forall a. C a => a -> [a]
sineListFast a
x =
   let sinx :: a
sinx  = a -> a
forall a. C a => a -> a
sin a
x
       cos2x :: a
cos2x = a
2 a -> a -> a
forall a. C a => a -> a -> a
- a
4a -> a -> a
forall a. C a => a -> a -> a
*a
sinxa -> a -> a
forall a. C a => a -> a -> a
*a
sinx
       --cos2x = 2 * cos (2*x)
       sines :: [a]
sines = (-a
sinx) a -> [a] -> [a]
forall a. a -> [a] -> [a]
: a
sinx a -> [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
y1 a
y0 -> a
cos2x a -> a -> a
forall a. C a => a -> a -> a
* a
y0 a -> a -> a
forall a. C a => a -> a -> a
- a
y1) [a]
sines ([a] -> [a]
forall a. HasCallStack => [a] -> [a]
tail [a]
sines)
   in  [a] -> [a]
forall a. HasCallStack => [a] -> [a]
tail [a]
sines

makeSines :: (Trans.C a) =>
   Int -> [a]
makeSines :: forall a. C a => Int -> [a]
makeSines Int
order =
   Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take (String -> Int -> Int
checkedHalf String
"makeSines" Int
order) (a -> [a]
forall a. C a => a -> [a]
sineList (a
forall a. C a => a
pi a -> a -> a
forall a. C a => a -> a -> a
/ Int -> a
forall a b. (C a, C b) => a -> b
fromIntegral (Int
2Int -> Int -> Int
forall a. C a => a -> a -> a
*Int
order)))

partialRatio :: (Trans.C a) =>
   Int -> a -> a
partialRatio :: forall a. C a => Int -> a -> a
partialRatio Int
order a
ratio =
   (a
1a -> a -> a
forall a. C a => a -> a -> a
/a
ratioa -> Integer -> a
forall a. C a => a -> Integer -> a
^Integer
2a -> a -> a
forall a. C a => a -> a -> a
-a
1) a -> a -> a
forall a. C a => a -> a -> a
** (-a
1 a -> a -> a
forall a. C a => a -> a -> a
/ Int -> a
forall a b. (C a, C b) => a -> b
fromIntegral (Int
2Int -> Int -> Int
forall a. C a => a -> a -> a
*Int
order))



_partialLowpassParameterInstable, partialLowpassParameter :: (Trans.C a) =>
   a -> a -> a -> Filt2.Parameter a

{- must handle infinite values when 'freq' approaches 0.5 -}
_partialLowpassParameterInstable :: forall a. C a => a -> a -> a -> Parameter a
_partialLowpassParameterInstable a
ratio a
freq a
sinw =
   let wc :: a
wc    = a
ratio a -> a -> a
forall a. C a => a -> a -> a
* a -> a
forall a. C a => a -> a
tan (a
forall a. C a => a
pia -> a -> a
forall a. C a => a -> a -> a
*a
freq)
       sinw2 :: a
sinw2 = a
2 a -> a -> a
forall a. C a => a -> a -> a
* a
wc a -> a -> a
forall a. C a => a -> a -> a
* a
sinw
       wc2 :: a
wc2   = a
wc a -> a -> a
forall a. C a => a -> a -> a
* a
wc
       denom :: a
denom = a
wc2 a -> a -> a
forall a. C a => a -> a -> a
+ a
sinw2 a -> a -> a
forall a. C a => a -> a -> a
+ a
1
       c0 :: a
c0    = a
wc2 a -> a -> a
forall a. C a => a -> a -> a
/ a
denom
   in  a -> a -> a -> a -> a -> Parameter a
forall a. a -> a -> a -> a -> a -> Parameter a
Filt2.Parameter a
c0 (a
2a -> a -> a
forall a. C a => a -> a -> a
*a
c0) a
c0
          (a
2a -> a -> a
forall a. C a => a -> a -> a
*(a
1a -> a -> a
forall a. C a => a -> a -> a
-a
wc2)a -> a -> a
forall a. C a => a -> a -> a
/a
denom) ((-a
wc2a -> a -> a
forall a. C a => a -> a -> a
+a
sinw2a -> a -> a
forall a. C a => a -> a -> a
-a
1)a -> a -> a
forall a. C a => a -> a -> a
/a
denom)

-- using ratio disallows simplification by trigonometric Pythagoras' theorem
partialLowpassParameter :: forall a. C a => a -> a -> a -> Parameter a
partialLowpassParameter a
ratio a
freq =
   let phi :: a
phi      = a
forall a. C a => a
pia -> a -> a
forall a. C a => a -> a -> a
*a
freq
       rsin2phi :: a
rsin2phi = a
ratio a -> a -> a
forall a. C a => a -> a -> a
* a -> a
forall a. C a => a -> a
sin (a
2a -> a -> a
forall a. C a => a -> a -> a
*a
phi)
       cosphi2 :: a
cosphi2  = a -> a
forall a. C a => a -> a
cos a
phi a -> Integer -> a
forall a. C a => a -> Integer -> a
^ Integer
2
       c0d :: a
c0d      = (a
ratio a -> a -> a
forall a. C a => a -> a -> a
* a -> a
forall a. C a => a -> a
sin a
phi) a -> Integer -> a
forall a. C a => a -> Integer -> a
^ Integer
2
       d1d :: a
d1d      = (a
cosphi2 a -> a -> a
forall a. C a => a -> a -> a
- a
c0d) a -> a -> a
forall a. C a => a -> a -> a
* a
2
       ratsin :: a
ratsin   = a
cosphi2 a -> a -> a
forall a. C a => a -> a -> a
+ a
c0d
   in  \a
sinw ->
          let rsinsin :: a
rsinsin = a
rsin2phi a -> a -> a
forall a. C a => a -> a -> a
* a
sinw
              denom :: a
denom   = a
rsinsin a -> a -> a
forall a. C a => a -> a -> a
+ a
ratsin
              d2d :: a
d2d     = a
rsinsin a -> a -> a
forall a. C a => a -> a -> a
- a
ratsin
              c0 :: a
c0      = a
c0d a -> a -> a
forall a. C a => a -> a -> a
/ a
denom
              d1 :: a
d1      = a
d1d a -> a -> a
forall a. C a => a -> a -> a
/ a
denom
              d2 :: a
d2      = a
d2d a -> a -> a
forall a. C a => a -> a -> a
/ a
denom
          in  a -> a -> a -> a -> a -> Parameter a
forall a. a -> a -> a -> a -> a -> Parameter a
Filt2.Parameter a
c0 (a
2a -> a -> a
forall a. C a => a -> a -> a
*a
c0) a
c0 a
d1 a
d2


-- * use second order filter parameters for control

type Parameter a = Cascade.Parameter a

{-# INLINE parameter #-}
parameter ::
   (Trans.C a, Storable a) =>
   Passband -> Int -> Pole a -> Parameter a
parameter :: forall a.
(C a, Storable a) =>
Passband -> Int -> Pole a -> Parameter a
parameter Passband
kind Int
order =
   -- I hope that the 'let' is floated out of a 'map'
   let sinesVec :: Vector a
sinesVec = [a] -> Vector a
forall a. Storable a => [a] -> Vector a
SV.pack (Int -> [a]
forall a. C a => Int -> [a]
makeSines Int
order)
   in  \ (Pole a
ratio a
freq) ->
           Vector (Parameter a) -> Parameter a
forall a. Vector (Parameter a) -> Parameter a
Cascade.Parameter (Vector (Parameter a) -> Parameter a)
-> Vector (Parameter a) -> Parameter a
forall a b. (a -> b) -> a -> b
$
           (a -> Parameter a) -> Vector a -> Vector (Parameter a)
forall a b.
(Storable a, Storable b) =>
(a -> b) -> Vector a -> Vector b
SV.map
              (\a
sinw ->
                 Passband -> a -> a -> a -> Parameter a
forall a. C a => Passband -> a -> a -> a -> Parameter a
partialParameter Passband
kind (Int -> a -> a
forall a. C a => Int -> a -> a
partialRatio Int
order a
ratio) a
sinw a
freq) (Vector a -> Vector (Parameter a))
-> Vector a -> Vector (Parameter a)
forall a b. (a -> b) -> a -> b
$
           Vector a
sinesVec

partialParameter ::
   Trans.C a =>
   Passband -> a -> a -> a -> Filt2.Parameter a
partialParameter :: forall a. C a => Passband -> a -> a -> a -> Parameter a
partialParameter Passband
kind a
partRatio a
sinw a
freq =
   Passband -> (a -> Parameter a) -> a -> Parameter a
forall a. C a => Passband -> (a -> Parameter a) -> a -> Parameter a
Filt2.adjustPassband Passband
kind
      ((a -> a -> Parameter a) -> a -> a -> Parameter a
forall a b c. (a -> b -> c) -> b -> a -> c
flip (a -> a -> a -> Parameter a
forall a. C a => a -> a -> a -> Parameter a
partialLowpassParameter a
partRatio) a
sinw)
      a
freq

{-# INLINE modifier #-}
modifier ::
   (Ring.C a, Module.C a v, Storable a, Storable v) =>
   Int ->
   Modifier.Simple (Cascade.State v) (Parameter a) v v
modifier :: forall a v.
(C a, C a v, Storable a, Storable v) =>
Int -> Simple (State v) (Parameter a) v v
modifier =
   Int -> Simple (State v) (Parameter a) v v
forall a v.
(C a, C a v, Storable a, Storable v) =>
Int -> Simple (State v) (Parameter a) v v
Cascade.modifier

{-# INLINE causal #-}
causal :: (Ring.C a, Module.C a v, Storable a, Storable v) =>
   Int ->
   Causal.T (Parameter a, v) v
causal :: forall a v.
(C a, C a v, Storable a, Storable v) =>
Int -> T (Parameter a, v) v
causal Int
order =
   Int -> T (Parameter a, v) v
forall a v.
(C a, C a v, Storable a, Storable v) =>
Int -> T (Parameter a, v) v
Cascade.causal (String -> Int -> Int
checkedHalf String
"causal" Int
order)


{-# INLINE checkedHalf #-}
checkedHalf :: String -> Int -> Int
checkedHalf :: String -> Int -> Int
checkedHalf String
funcName Int
order =
   let (Int
order2,Int
r) = Int -> Int -> (Int, Int)
forall a. C a => a -> a -> (a, a)
divMod Int
order Int
2
   in  if Int
rInt -> Int -> Bool
forall a. Eq a => a -> a -> Bool
==Int
0
         then Int
order2
         else String -> Int
forall a. HasCallStack => String -> a
error (String -> Int) -> String -> Int
forall a b. (a -> b) -> a -> b
$ String
"Butterworth."String -> String -> String
forall a. [a] -> [a] -> [a]
++String
funcNameString -> String -> String
forall a. [a] -> [a] -> [a]
++String
": order must be even"

{-
lowpassCausal, highpassCausal :: (Trans.C a, Module.C a v) =>
   Int -> Causal.T (Parameter a, v) v
lowpassCausal  = causal Lowpass
highpassCausal = causal Highpass

lowpass, highpass :: (Trans.C a, Module.C a v) =>
   Int -> Sig.T (Parameter a) -> Sig.T v -> Sig.T v
lowpass  = run Lowpass
highpass = run Highpass
-}


-- * directly use frequency as control parameter

{- |
When called as @runPole kind order ratio freqs@,
the filter amplifies frequency 0 with factor 1
and frequency @freq@ with factor @ratio@.

It uses the frequency and ratio information directly
and thus cannot benefit from efficient parameter interpolation
(asynchronous run of a ControlledProcess).
-}
runPole :: (Trans.C a, Module.C a v) =>
   Passband -> Int -> Sig.T a -> Sig.T a -> Sig.T v -> Sig.T v
runPole :: forall a v.
(C a, C a v) =>
Passband -> Int -> T a -> T a -> T v -> T v
runPole Passband
kind Int
order T a
ratios T a
freqs =
   let makePartialFilter :: a -> T v -> T v
makePartialFilter a
s =
          T (Parameter a) -> T v -> T v
forall a v. (C a, C a v) => T (Parameter a) -> T v -> T v
Filt2.run (T (Parameter a) -> T v -> T v) -> T (Parameter a) -> T v -> T v
forall a b. (a -> b) -> a -> b
$
          (a -> a -> Parameter a) -> T a -> T a -> T (Parameter a)
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith
             (\a
ratio a
freq ->
                Passband -> a -> a -> a -> Parameter a
forall a. C a => Passband -> a -> a -> a -> Parameter a
partialParameter Passband
kind (Int -> a -> a
forall a. C a => Int -> a -> a
partialRatio Int
order a
ratio) a
s a
freq)
             T a
ratios T a
freqs
   in  ((T v -> T v) -> (T v -> T v) -> T v -> T v)
-> (T v -> T v) -> [T v -> T v] -> T v -> T v
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl (T v -> T v) -> (T v -> T v) -> T v -> T v
forall b c a. (b -> c) -> (a -> b) -> a -> c
(.) T v -> T v
forall a. a -> a
id ((a -> T v -> T v) -> T a -> [T v -> T v]
forall a b. (a -> b) -> [a] -> [b]
map a -> T v -> T v
forall {v}. C a v => a -> T v -> T v
makePartialFilter (Int -> T a
forall a. C a => Int -> [a]
makeSines Int
order))

causalPole :: (Trans.C a, Module.C a v) =>
   Passband -> Int -> Causal.T (Pole a, v) v
causalPole :: forall a v. (C a, C a v) => Passband -> Int -> T (Pole a, v) v
causalPole Passband
kind Int
order =
   let {-# INLINE makePartialFilter #-}
       makePartialFilter :: a -> T (Pole a, c) c
makePartialFilter a
s =
          T (Pole a) (Parameter a) -> T (Pole a, c) (Parameter a, c)
forall b c d. T b c -> T (b, d) (c, d)
forall (a :: * -> * -> *) b c d.
Arrow a =>
a b c -> a (b, d) (c, d)
Causal.first
             ((Pole a -> Parameter a) -> T (Pole a) (Parameter a)
forall a b. (a -> b) -> T a b
Causal.map (\(Pole a
ratio a
freq) ->
                Passband -> a -> a -> a -> Parameter a
forall a. C a => Passband -> a -> a -> a -> Parameter a
partialParameter Passband
kind (Int -> a -> a
forall a. C a => Int -> a -> a
partialRatio Int
order a
ratio) a
s a
freq))
          T (Pole a, c) (Parameter a, c)
-> T (Parameter a, c) c -> T (Pole a, c) c
forall {k} (cat :: k -> k -> *) (a :: k) (b :: k) (c :: k).
Category cat =>
cat a b -> cat b c -> cat a c
>>>
          T (Parameter a, c) c
forall a v. (C a, C a v) => T (Parameter a, v) v
Filt2.causal
   in  [T (Pole a, v) v] -> T (Pole a, v) v
forall c x. [T (c, x) x] -> T (c, x) x
Causal.chainControlled ([T (Pole a, v) v] -> T (Pole a, v) v)
-> [T (Pole a, v) v] -> T (Pole a, v) v
forall a b. (a -> b) -> a -> b
$ (a -> T (Pole a, v) v) -> [a] -> [T (Pole a, v) v]
forall a b. (a -> b) -> [a] -> [b]
map a -> T (Pole a, v) v
forall {a} {c}. (C a, C a c) => a -> T (Pole a, c) c
makePartialFilter ([a] -> [T (Pole a, v) v]) -> [a] -> [T (Pole a, v) v]
forall a b. (a -> b) -> a -> b
$ Int -> [a]
forall a. C a => Int -> [a]
makeSines Int
order


lowpassCausalPole, highpassCausalPole :: (Trans.C a, Module.C a v) =>
   Int -> Causal.T (Pole a, v) v
lowpassCausalPole :: forall a v. (C a, C a v) => Int -> T (Pole a, v) v
lowpassCausalPole  = Passband -> Int -> T (Pole a, v) v
forall a v. (C a, C a v) => Passband -> Int -> T (Pole a, v) v
causalPole Passband
Lowpass
highpassCausalPole :: forall a v. (C a, C a v) => Int -> T (Pole a, v) v
highpassCausalPole = Passband -> Int -> T (Pole a, v) v
forall a v. (C a, C a v) => Passband -> Int -> T (Pole a, v) v
causalPole Passband
Highpass

lowpassPole, highpassPole :: (Trans.C a, Module.C a v) =>
   Int -> Sig.T a -> Sig.T a -> Sig.T v -> Sig.T v
lowpassPole :: forall a v. (C a, C a v) => Int -> T a -> T a -> T v -> T v
lowpassPole  = Passband -> Int -> T a -> T a -> T v -> T v
forall a v.
(C a, C a v) =>
Passband -> Int -> T a -> T a -> T v -> T v
runPole Passband
Lowpass
highpassPole :: forall a v. (C a, C a v) => Int -> T a -> T a -> T v -> T v
highpassPole = Passband -> Int -> T a -> T a -> T v -> T v
forall a v.
(C a, C a v) =>
Passband -> Int -> T a -> T a -> T v -> T v
runPole Passband
Highpass