{-# 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
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
showList :: [Parameter a] -> ShowS
$cshowList :: forall a. Show a => [Parameter a] -> ShowS
show :: Parameter a -> String
$cshow :: forall a. Show a => Parameter a -> String
showsPrec :: Int -> Parameter a -> ShowS
$cshowsPrec :: forall a. Show a => Int -> 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 = forall a b. (a -> b) -> [a] -> [b]
map forall a. C a => a -> a
negate
           (forall a b. (a -> b) -> [a] -> [b]
map (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 = forall a b. (a -> b) -> [a] -> [b]
map (\Int
j ->
            forall a b. (a -> b) -> [a] -> [b]
map (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  forall a. [a] -> Parameter a
Parameter (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
2forall a. C a => a -> a -> a
*forall a. C a => a
pi forall a. C a => a -> a -> a
* a
frequency
        omegas :: [a]
omegas = forall a. (a -> a) -> a -> [a]
iterate (a
omegaforall a. C a => a -> a -> a
+) a
omega
        denom :: T a
denom = T a
1forall a. C a => a -> a -> a
+forall a. C a => [a] -> a
sum (forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\a
k a
w -> a
kforall a. C a => a -> a -> a
*forall a. C a => a -> a
cos a
w forall a. a -> a -> T a
+: a
kforall a. C a => a -> a -> a
*forall a. C a => a -> a
sin a
w) [a]
ks [a]
omegas)
    in  a
2 forall a. C a => a -> a -> a
* forall a. (C a, C a) => T a -> a
Complex.phase T a
denom forall a. C a => a -> a -> a
- a
omegaforall a. C a => a -> a -> a
*(forall a b. (C a, C b) => a -> b
fromIntegral (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) = forall a. C a => a -> Int -> (a, a)
integrateScrewExp a
r (Int
kforall a. C a => a -> a -> a
+Int
jforall a. C a => a -> a -> a
-Int
order)
    in  a
2 forall a. C a => a -> a -> a
* (forall a b. (a, b) -> a
fst (forall a. C a => a -> Int -> (a, a)
integrateScrewExp a
r (Int
kforall a. C a => a -> a -> a
-Int
j)) forall a. C a => a -> a -> a
-
              (forall a. C a => a -> a
cos a
phase forall a. C a => a -> a -> a
* a
intCos forall a. C a => a -> a -> a
+ forall a. C a => a -> a
sin a
phase 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 = forall a. C a => a -> T a
cis (forall a b. (C a, C b) => a -> b
fromIntegral Int
k forall a. C a => a -> a -> a
* a
omega) forall a. C a => a -> a -> a
-
                       forall a. C a => a -> T a
cis a
phase forall a. C a => a -> a -> a
* forall a. C a => a -> T a
cis (forall a b. (C a, C b) => a -> b
fromIntegral (Int
orderforall a. C a => a -> a -> a
-Int
k) forall a. C a => a -> a -> a
* a
omega)
        z1 :: T a
z1 = forall a. C a => a -> T a
cis (forall a b. (C a, C b) => a -> b
fromIntegral Int
j forall a. C a => a -> a -> a
* a
omega) forall a. C a => a -> a -> a
-
                       forall a. C a => a -> T a
cis a
phase forall a. C a => a -> a -> a
* forall a. C a => a -> T a
cis (forall a b. (C a, C b) => a -> b
fromIntegral (Int
orderforall a. C a => a -> a -> a
-Int
j) forall a. C a => a -> a -> a
* a
omega)
    in  forall a. T a -> a
real T a
z0 forall a. C a => a -> a -> a
* forall a. T a -> a
real T a
z1 forall a. C a => a -> a -> a
+ forall a. T a -> a
imag T a
z0 forall a. C a => a -> 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 = forall a b. (C a, C b) => a -> b
fromIntegral Int
kInt
        q :: a
q = (forall a. C a => a -> a
exp (a
2forall a. C a => a -> a -> a
*forall a. C a => a
piforall a. C a => a -> a -> a
*a
r) forall a. C a => a -> a -> a
- a
1) forall a. C a => a -> a -> a
/ (a
rforall a. C a => a -> Integer -> a
^Integer
2 forall a. C a => a -> a -> a
+ a
kforall a. C a => a -> Integer -> a
^Integer
2)
    in  (a
rforall a. C a => a -> a -> a
*a
q, -a
kforall 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 = forall a b. (a -> b) -> [a] -> [b]
map (\Int
k -> a
lo forall a. C a => a -> a -> a
+ (a
hiforall a. C a => a -> a -> a
-a
lo) forall a. C a => a -> a -> a
* forall a b. (C a, C b) => a -> b
fromIntegral Int
k forall a. C a => a -> a -> a
/ forall a b. (C a, C b) => a -> b
fromIntegral Int
n)
                 [Int
1..(Int
nforall a. C a => a -> a -> a
-Int
1)]
    in  ((a
hiforall a. C a => a -> a -> a
-a
lo) forall a. C a => a -> a -> a
/ forall a b. (C a, C b) => a -> b
fromIntegral Int
n) forall a v. C a v => a -> v -> v
*>
        (forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl forall a. C a => a -> a -> a
(+) ((a
1forall a. C a => a -> a -> a
/a
2 forall a. a -> a -> a
`asTypeOf` a
lo) forall a v. C a v => a -> v -> v
*> (a -> v
f a
lo forall a. C a => a -> a -> a
+ a -> v
f a
hi))
               (forall a b. (a -> b) -> [a] -> [b]
map a -> v
f [a]
xs))