{-# 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
-}
module Synthesizer.Plain.Filter.Recursive.AllpassPoly where

import qualified Algebra.Module                as Module
import qualified Algebra.RealTranscendental    as RealTrans
import qualified Algebra.Transcendental        as Trans
import qualified Algebra.Field                 as Field
import qualified Algebra.ZeroTestable          as ZeroTestable

import Number.Complex (cis,(+:),real,imag)
import qualified Number.Complex as Complex
import Orthogonals(Scalar,one_ket_solution)

import qualified Prelude as P
import NumericPrelude.Numeric
import NumericPrelude.Base



newtype Parameter a = Parameter [a]
   deriving Int -> Parameter a -> ShowS
[Parameter a] -> ShowS
Parameter a -> String
(Int -> Parameter a -> ShowS)
-> (Parameter a -> String)
-> ([Parameter a] -> ShowS)
-> Show (Parameter a)
forall a. Show a => Int -> Parameter a -> ShowS
forall a. Show a => [Parameter a] -> ShowS
forall a. Show a => Parameter a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: forall a. Show a => Int -> Parameter a -> ShowS
showsPrec :: Int -> Parameter a -> ShowS
$cshow :: forall a. Show a => Parameter a -> String
show :: Parameter a -> String
$cshowList :: forall a. Show a => [Parameter a] -> ShowS
showList :: [Parameter a] -> ShowS
Show

{- | Compute coefficients for an allpass that shifts low frequencies
     by approximately the shift you want.
     To achieve this we solve a linear least squares problem,
     where low frequencies are more weighted than high ones.
     The output is a list of coefficients for an arbitrary order allpass. -}
shiftParam :: (Scalar a, P.Fractional a, Trans.C a) =>
   Int -> a -> a -> Parameter a
shiftParam :: forall a.
(Scalar a, Fractional a, C a) =>
Int -> a -> a -> Parameter a
shiftParam Int
order a
weight a
phase =
    let {- construct matrix for normal equations -}
        normalVector :: [a]
normalVector = (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map a -> a
forall a. C a => a -> a
negate
           ((Int -> a) -> [Int] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a -> Int -> a -> Int -> Int -> a
forall a. C a => a -> Int -> a -> Int -> Int -> a
scalarProdScrewExp a
weight Int
order a
phase Int
0) [Int
1..Int
order])
        normalMatrix :: [[a]]
normalMatrix = (Int -> [a]) -> [Int] -> [[a]]
forall a b. (a -> b) -> [a] -> [b]
map (\Int
j ->
            (Int -> a) -> [Int] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a -> Int -> a -> Int -> Int -> a
forall a. C a => a -> Int -> a -> Int -> Int -> a
scalarProdScrewExp a
weight Int
order a
phase Int
j) [Int
1..Int
order]) [Int
1..Int
order]
    in  [a] -> Parameter a
forall a. [a] -> Parameter a
Parameter ([[a]] -> [a] -> [a]
forall a. (Scalar a, Fractional a) => [[a]] -> [a] -> [a]
one_ket_solution [[a]]
normalMatrix [a]
normalVector)

{-
  GNUPlot.plotFunc (GNUPlot.linearScale 500 (0,1)) ((fwrap (-pi,pi)).(makePhase (shiftParam 6 (-6) (-pi/2::Double))))
-}
makePhase :: (RealTrans.C a, ZeroTestable.C a) => Parameter a -> a -> a
makePhase :: forall a. (C a, C a) => Parameter a -> a -> a
makePhase (Parameter [a]
ks) a
frequency =
    let omega :: a
omega  = a
2a -> a -> a
forall a. C a => a -> a -> a
*a
forall a. C a => a
pi a -> a -> a
forall a. C a => a -> a -> a
* a
frequency
        omegas :: [a]
omegas = (a -> a) -> a -> [a]
forall a. (a -> a) -> a -> [a]
iterate (a
omegaa -> a -> a
forall a. C a => a -> a -> a
+) a
omega
        denom :: T a
denom = T a
1T a -> T a -> T a
forall a. C a => a -> a -> a
+[T a] -> T a
forall a. C a => [a] -> a
sum ((a -> a -> T a) -> [a] -> [a] -> [T a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\a
k a
w -> a
ka -> a -> a
forall a. C a => a -> a -> a
*a -> a
forall a. C a => a -> a
cos a
w a -> a -> T a
forall a. a -> a -> T a
+: a
ka -> a -> a
forall a. C a => a -> a -> a
*a -> a
forall a. C a => a -> a
sin a
w) [a]
ks [a]
omegas)
    in  a
2 a -> a -> a
forall a. C a => a -> a -> a
* T a -> a
forall a. (C a, C a) => T a -> a
Complex.phase T a
denom a -> a -> a
forall a. C a => a -> a -> a
- a
omegaa -> a -> a
forall a. C a => a -> a -> a
*(Int -> a
forall a b. (C a, C b) => a -> b
fromIntegral ([a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
ks))

{- integrate (0,2*pi) (\omega -> exp (r*omega) * screwProd order phase k j omega) -}
scalarProdScrewExp :: Trans.C a => a -> Int -> a -> Int -> Int -> a
scalarProdScrewExp :: forall a. C a => a -> Int -> a -> Int -> Int -> a
scalarProdScrewExp a
r Int
order a
phase Int
k Int
j =
    let (a
intCos,a
intSin) = a -> Int -> (a, a)
forall a. C a => a -> Int -> (a, a)
integrateScrewExp a
r (Int
kInt -> Int -> Int
forall a. C a => a -> a -> a
+Int
jInt -> Int -> Int
forall a. C a => a -> a -> a
-Int
order)
    in  a
2 a -> a -> a
forall a. C a => a -> a -> a
* ((a, a) -> a
forall a b. (a, b) -> a
fst (a -> Int -> (a, a)
forall a. C a => a -> Int -> (a, a)
integrateScrewExp a
r (Int
kInt -> Int -> Int
forall a. C a => a -> a -> a
-Int
j)) a -> a -> a
forall a. C a => a -> a -> a
-
              (a -> a
forall a. C a => a -> a
cos a
phase a -> a -> a
forall a. C a => a -> a -> a
* a
intCos a -> a -> a
forall a. C a => a -> a -> a
+ a -> a
forall a. C a => a -> a
sin a
phase a -> a -> a
forall a. C a => a -> a -> a
* a
intSin))

screwProd :: Trans.C a => Int -> a -> Int -> Int -> a -> a
screwProd :: forall a. C a => Int -> a -> Int -> Int -> a -> a
screwProd Int
order a
phase Int
k Int
j a
omega =
    let z0 :: T a
z0 = a -> T a
forall a. C a => a -> T a
cis (Int -> a
forall a b. (C a, C b) => a -> b
fromIntegral Int
k a -> a -> a
forall a. C a => a -> a -> a
* a
omega) T a -> T a -> T a
forall a. C a => a -> a -> a
-
                       a -> T a
forall a. C a => a -> T a
cis a
phase T a -> T a -> T a
forall a. C a => a -> a -> a
* a -> T a
forall a. C a => a -> T a
cis (Int -> a
forall a b. (C a, C b) => a -> b
fromIntegral (Int
orderInt -> Int -> Int
forall a. C a => a -> a -> a
-Int
k) a -> a -> a
forall a. C a => a -> a -> a
* a
omega)
        z1 :: T a
z1 = a -> T a
forall a. C a => a -> T a
cis (Int -> a
forall a b. (C a, C b) => a -> b
fromIntegral Int
j a -> a -> a
forall a. C a => a -> a -> a
* a
omega) T a -> T a -> T a
forall a. C a => a -> a -> a
-
                       a -> T a
forall a. C a => a -> T a
cis a
phase T a -> T a -> T a
forall a. C a => a -> a -> a
* a -> T a
forall a. C a => a -> T a
cis (Int -> a
forall a b. (C a, C b) => a -> b
fromIntegral (Int
orderInt -> Int -> Int
forall a. C a => a -> a -> a
-Int
j) a -> a -> a
forall a. C a => a -> a -> a
* a
omega)
    in  T a -> a
forall a. T a -> a
real T a
z0 a -> a -> a
forall a. C a => a -> a -> a
* T a -> a
forall a. T a -> a
real T a
z1 a -> a -> a
forall a. C a => a -> a -> a
+ T a -> a
forall a. T a -> a
imag T a
z0 a -> a -> a
forall a. C a => a -> a -> a
* T a -> a
forall a. T a -> a
imag T a
z1

{- integrate (0,2*pi) (\omega -> (exp (r*omega) +: 0) * cis (k*omega)) -}
integrateScrewExp :: Trans.C a => a -> Int -> (a,a)
integrateScrewExp :: forall a. C a => a -> Int -> (a, a)
integrateScrewExp a
r Int
kInt =
    let k :: a
k = Int -> a
forall a b. (C a, C b) => a -> b
fromIntegral Int
kInt
        q :: a
q = (a -> a
forall a. C a => a -> a
exp (a
2a -> a -> a
forall a. C a => a -> a -> a
*a
forall a. C a => a
pia -> a -> a
forall a. C a => a -> a -> a
*a
r) a -> a -> a
forall a. C a => a -> a -> a
- a
1) a -> a -> a
forall a. C a => a -> a -> a
/ (a
ra -> Integer -> a
forall a. C a => a -> Integer -> a
^Integer
2 a -> a -> a
forall a. C a => a -> a -> a
+ a
ka -> Integer -> a
forall a. C a => a -> Integer -> a
^Integer
2)
    in  (a
ra -> a -> a
forall a. C a => a -> a -> a
*a
q, -a
ka -> a -> a
forall a. C a => a -> a -> a
*a
q)

{- Should be moved to NumericPrelude.Numeric -}
integrateNum :: (Field.C a, Module.C a v) => Int -> (a,a) -> (a->v) -> v
integrateNum :: forall a v. (C a, C a v) => Int -> (a, a) -> (a -> v) -> v
integrateNum Int
n (a
lo,a
hi) a -> v
f =
    let xs :: [a]
xs = (Int -> a) -> [Int] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (\Int
k -> a
lo a -> a -> a
forall a. C a => a -> a -> a
+ (a
hia -> a -> a
forall a. C a => a -> a -> a
-a
lo) a -> a -> a
forall a. C a => a -> a -> a
* Int -> a
forall a b. (C a, C b) => a -> b
fromIntegral Int
k a -> a -> a
forall a. C a => a -> a -> a
/ Int -> a
forall a b. (C a, C b) => a -> b
fromIntegral Int
n)
                 [Int
1..(Int
nInt -> Int -> Int
forall a. C a => a -> a -> a
-Int
1)]
    in  ((a
hia -> a -> a
forall a. C a => a -> a -> a
-a
lo) a -> a -> a
forall a. C a => a -> a -> a
/ Int -> a
forall a b. (C a, C b) => a -> b
fromIntegral Int
n) a -> v -> v
forall a v. C a v => a -> v -> v
*>
        ((v -> v -> v) -> v -> [v] -> v
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl v -> v -> v
forall a. C a => a -> a -> a
(+) ((a
1a -> a -> a
forall a. C a => a -> a -> a
/a
2 a -> a -> a
forall a. a -> a -> a
`asTypeOf` a
lo) a -> v -> v
forall a v. C a v => a -> v -> v
*> (a -> v
f a
lo v -> v -> v
forall a. C a => a -> a -> a
+ a -> v
f a
hi))
               ((a -> v) -> [a] -> [v]
forall a b. (a -> b) -> [a] -> [b]
map a -> v
f [a]
xs))