```{-# OPTIONS -fglasgow-exts -fno-implicit-prelude #-}
module Filter.TwoWay where

import Filter.Basic

import Synthesizer.Plain.Interpolation (minLength,)
import qualified Synthesizer.Plain.Interpolation as Ip
import qualified Synthesizer.Plain.Interpolation as Interpolation

import Algebra.Module(linearComb,(*>))

import qualified Algebra.Module    as Module
import qualified Algebra.RealField as RealField
import qualified Algebra.Ring      as Ring

import Number.Complex(cis)
import Synthesizer.Utility(nest)
import qualified Data.List as List

import PreludeBase hiding (take)
import NumericPrelude

{-| A TwoWay.Signal stores values of the past and the future -}
data Signal v = Signal {past, future :: [v]}
deriving (Show, Eq)

{-| Take n values starting from time zero.
If you want clips from elsewhere,
call 'take' after 'delay'. -}
take :: Int -> Signal v -> [v]
take n (Signal _ x) = List.take n x

zipSignalWith :: (a -> b -> c) -> Signal a -> Signal b -> Signal c
zipSignalWith f (Signal xPast xFuture) (Signal yPast yFuture) =
(Signal (zipWith f xPast yPast) (zipWith f xFuture yFuture))

{-| Take the value at time zero. -}
origin :: Ring.C a => Signal a -> a
origin (Signal _ (x:_)) = x
origin _ = 0

{-| A signal that consists entirely of ones -}
ones :: Ring.C a => Signal a
ones = Signal (repeat 1) (repeat 1)

{-| shift signal in time,
keep all values but if required pad with zeros -}
Int -> Signal v -> Signal v
delay = delayGen delayOnce

delayPad :: v -> Int -> Signal v -> Signal v

{-| shift signal in time,
zero values at either ends will be flushed -}
delayOpt :: (Eq v, Additive.C v) =>
Int -> Signal v -> Signal v
delayOpt = delayGen delayOptOnce

{-| Delay by one sample. -}
Signal v -> Signal v
--delayOnce (Signal []     []) = ([],[])
delayOnce (Signal []     ys) = Signal [] (zero:ys)
delayOnce (Signal (x:xs) ys) = Signal xs (x:ys)

delayPadOnce :: v -> Signal v -> Signal v
--delayPadOnce _ (Signal []     []) = ([],[])
delayPadOnce z (Signal []     ys) = Signal [] (z:ys)
delayPadOnce _ (Signal (x:xs) ys) = Signal xs (x:ys)

delayOptOnce :: (Eq v, Additive.C v) =>
Signal v -> Signal v
--delayOptOnce (Signal []     []) = Signal [] []
delayOptOnce (Signal []     ys) = Signal [] (zero:ys)
delayOptOnce (Signal (x:xs) []) = Signal xs (if x==zero then [] else x:[])
delayOptOnce (Signal (x:xs) ys) = Signal xs (x:ys)

{-| General routine that supports delaying and prefetching
using a general one-sample delaying routine. -}
delayGen :: (Signal v -> Signal v) ->
Int -> Signal v -> Signal v
{- Using this optimization applications of recursive filters
with zero initial conditions represented by an empty list will fail.
This is because in this case the value of the first item of the future list
depends on the first item of the input future list,
whereas normally the first future value depends on no input future value,
at all.
delayGen _ _ (Signal [] []) = Signal [] []
cf. the next example -}
delayGen delOnce t =
if t < 0
then reverseTwoWay . nest (negate t) delOnce . reverseTwoWay
else nest t delOnce

reverseTwoWay :: Signal v -> Signal v
reverseTwoWay (Signal x y) = Signal y x

zero                                = Signal zero zero
(+)   (Signal y0 y1) (Signal x0 x1) = Signal (y0 + x0)   (y1 + x1)
(-)   (Signal y0 y1) (Signal x0 x1) = Signal (y0 - x0)   (y1 - x1)
negate               (Signal x0 x1) = Signal (negate x0) (negate x1)

instance (Module.C a v) => Module.C a (Signal v) where
(*>)  s      (Signal x0 x1) = Signal (s *> x0) (s *> x1)

{-| for a Signal this means a reversion of the elements -}
flipPair :: (a,b) -> (b,a)
flipPair (x,y) = (y,x)

{- This example simulates what happens if you call
apply (Feedback (Serial []) (Prim (Delay 1)) []) ([],[1])
depending on the implementation of delayGen it may work or
loop infinitely when yFuture is computed.
It's even before the first element of yFuture is computed.
Note, that the equivalent
apply (Feedback (Serial []) (Prim (Delay 1)) [0]) ([],[1])
works! (i.e. set yPast = [0] -}
testDelayGen :: Signal Double
testDelayGen =
let yPast = []
x = Signal [] [1]
y = Signal yPast yFuture
Signal _ yFuture = delayOnce (x + y)
-- Signal _ yFuture = delayOptOnce (add x y)
-- Signal _ yFuture = delayGen delayOnce 1 (add x y)
in  Signal yPast (List.take 10 yFuture)

{-| Unmodulated non-recursive filter -}
nonRecursiveFilter :: Module.C a v =>
[a] -> Signal v -> Signal v
nonRecursiveFilter m x =
linearComb m (iterate delayOnce x)

{-| Modulated non-recursive filter.
The number of values before time 0 (past) or
the filter mask lengths must be at most finite. -}
nonRecursiveFilterMod :: Module.C a v =>
Signal [a] -> Signal v -> Signal v
nonRecursiveFilterMod (Signal mpre msuf) x =
let (pre, suf) = unzip (map (\(Signal a b) -> (a,b)) (iterate delayOnce x))
in  Signal (zipWith linearComb mpre pre) (zipWith linearComb msuf suf)

{-| Interpolation allowing negative frequencies,
but requires storage of all past values. -}
interpolatePaddedZero :: (Ord a, RealField.C a) =>
b -> Interpolation.T a b
-> a -> Signal a -> Signal b -> Signal b
interpolatePaddedZero z ip phase fs (Signal xPast xFuture) =
let (phInt, phFrac) = splitFraction phase
xPadded = Signal (xPast ++ repeat z) (xFuture ++ repeat z)
in  interpolateCore ip phFrac fs

interpolatePaddedCyclic :: (Ord a, RealField.C a) =>
Interpolation.T a b
-> a -> Signal a -> Signal b -> Signal b
interpolatePaddedCyclic ip phase fs (Signal xPast xFuture) =
let (phInt, phFrac) = splitFraction phase
xCyclic = xFuture ++ reverse xPast
in  interpolateCore ip phFrac fs
-- mod is for efficiency, only
(mod (Ip.offset ip - phInt) (length xCyclic))
(Signal (cycle (reverse xCyclic)) (cycle xCyclic)))

-- note that the extrapolation may miss some of the first and some of the last points
interpolatePaddedExtrapolation :: (Ord a, RealField.C a) =>
Interpolation.T a b
-> a -> Signal a -> Signal b -> Signal b
interpolatePaddedExtrapolation ip phase fs x =
interpolateCore ip (phase - fromIntegral (Ip.offset ip)) fs x

interpolateCore :: (Ord a, Ring.C a) =>
Interpolation.T a b -> a -> Signal a -> Signal b -> Signal b
interpolateCore ip phase (Signal freqPast freqFuture) x =
Signal (interpolateHalfWay ip (1-phase) freqPast
(reverseTwoWay x)))
(interpolateHalfWay ip phase freqFuture x)

interpolateHalfWay :: (Ord a, Ring.C a) =>
Interpolation.T a b -> a -> [a] -> Signal b -> [b]
interpolateHalfWay ip phase freqs (Signal xPast xFuture) =
if phase >= 1 && minLength (1+Ip.number ip) xFuture
then interpolateHalfWay ip (phase-1) freqs
(Signal (head xFuture : xPast) (tail xFuture))
else if phase < 0 && minLength 1 xPast
then interpolateHalfWay ip (phase + 1) freqs
(Signal (tail xPast) (head xPast : xFuture))
else Ip.func ip phase xFuture :
interpolateHalfWay ip (phase + head freqs) (tail freqs)
(Signal xPast xFuture)

{-| Description of a basic filter that can be used in larger networks. -}
data T t a v =
{-^ A static filter described by its mask -}
{-^ A modulated filter described by a list of masks -}
| FracDelay (Interpolation.T t v) t
{-^ Delay the signal by a fractional amount of samples.
This is achieved by interpolation. -}
| ModFracDelay (Interpolation.T t v) (Signal t)
{-^ Delay with varying delay times. -}
| Delay Int
{-^ Delay the signal by given amount of samples. -}
| Past [v]
{-^ Replace the past by the given one.
This must be present in each recursive filter cycle
to let the magic work! -}

instance Filter Signal T where

apply (Mask         m)     = nonRecursiveFilter m
apply (ModMask      m)     = nonRecursiveFilterMod m
apply (FracDelay    ip t)  = interpolatePaddedZero zero
ip (-t) ones
apply (ModFracDelay ip ts) = interpolatePaddedZero zero
ip (- origin ts) (ts - delay (-1) ts + ones)
apply (Delay        t)     = delay t
apply (Past         x)     = Signal x . future

{- This is in principle the same as for one way filters.
How can one merge them? -}
transferFunction (Mask        m) w = linearComb m (screw (negate w))
transferFunction (FracDelay _ t) w = cis (negate w * t)
transferFunction (Delay       t) w = cis (negate w * fromIntegral t)
transferFunction (Past        _) _ = 1
transferFunction _ _ =
error "transfer function can't be computed for modulated filters"
```