{-# LANGUAGE NoImplicitPrelude #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE FlexibleInstances #-}
module Synthesizer.Filter.TwoWay where

import Synthesizer.Filter.Basic

import qualified Synthesizer.Plain.Signal as Sig
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 qualified Algebra.Additive  as Additive

import Number.Complex(cis, )
import Data.Function.HT (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 -}
delay :: (Additive.C v) =>
   Int -> Signal v -> Signal v
delay = delayGen delayOnce

delayPad :: v -> Int -> Signal v -> Signal v
delayPad z = delayGen (delayPadOnce z)

{-| 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. -}
delayOnce :: (Additive.C v) =>
   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


instance (Additive.C v) => Additive.C (Signal v) where
   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
          (delayPad z (Ip.offset ip - phInt) xPadded)

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
          (delayPad (error "interpolate: infinite signal needs no zero padding")
                    (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
             (delayPadOnce (error "interpolateCore: infinite signal needs no zero padding")
                           (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 && Sig.lengthAtLeast (1+Ip.number ip) xFuture
    then interpolateHalfWay ip (phase-1) freqs
            (Signal (head xFuture : xPast) (tail xFuture))
    else if phase < 0 && Sig.lengthAtLeast 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 =
     Mask [a]
        {-^ A static filter described by its mask -}
   | ModMask (Signal [a])
        {-^ 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"