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 Show
shiftParam :: (Scalar a, P.Fractional a, Trans.C a) =>
Int -> a -> a -> Parameter a
shiftParam order weight phase =
let
normalVector = map negate
(map (scalarProdScrewExp weight order phase 0) [1..order])
normalMatrix = map (\j ->
map (scalarProdScrewExp weight order phase j) [1..order]) [1..order]
in Parameter (one_ket_solution normalMatrix normalVector)
makePhase :: (RealTrans.C a, ZeroTestable.C a) => Parameter a -> a -> a
makePhase (Parameter ks) frequency =
let omega = 2*pi * frequency
omegas = iterate (omega+) omega
denom = 1+sum (zipWith (\k w -> k*cos w +: k*sin w) ks omegas)
in 2 * Complex.phase denom omega*(fromIntegral (length ks))
scalarProdScrewExp :: Trans.C a => a -> Int -> a -> Int -> Int -> a
scalarProdScrewExp r order phase k j =
let (intCos,intSin) = integrateScrewExp r (k+jorder)
in 2 * (fst (integrateScrewExp r (kj))
(cos phase * intCos + sin phase * intSin))
screwProd :: Trans.C a => Int -> a -> Int -> Int -> a -> a
screwProd order phase k j omega =
let z0 = cis (fromIntegral k * omega)
cis phase * cis (fromIntegral (orderk) * omega)
z1 = cis (fromIntegral j * omega)
cis phase * cis (fromIntegral (orderj) * omega)
in real z0 * real z1 + imag z0 * imag z1
integrateScrewExp :: Trans.C a => a -> Int -> (a,a)
integrateScrewExp r kInt =
let k = fromIntegral kInt
q = (exp (2*pi*r) 1) / (r^2 + k^2)
in (r*q, k*q)
integrateNum :: (Field.C a, Module.C a v) => Int -> (a,a) -> (a->v) -> v
integrateNum n (lo,hi) f =
let xs = map (\k -> lo + (hilo) * fromIntegral k / fromIntegral n)
[1..(n1)]
in ((hilo) / fromIntegral n) *>
(foldl (+) ((1/2 `asTypeOf` lo) *> (f lo + f hi))
(map f xs))