{-# LANGUAGE Arrows, TemplateHaskell, BangPatterns, 
             ExistentialQuantification, FlexibleContexts, 
             FunctionalDependencies, ScopedTypeVariables #-}

-- | Haskell adaptation of some unit generators in CSound.
-- Convention: Optional arguments in some CSound unit generators
-- sometimes carry different semantics depending on the way the
-- generator is called.  They are encoded in algebraic datatypes
-- instead (see 'pluck' for example).  Single optional argument
-- is normally encoded in the Haskell Maybe type.

-- CSound's i-type is updated only once on every note's initialization
-- pass.  They are represented as unlifted arguments here
-- (i.e. non-signal).

-- Most unit generators in CSound take in a signal 'amp' as input,
-- essentially amplifies its result by 'amp'.  Since this feature is
-- easily expressed using arrow syntax, we will leave it to the
-- library user instead of requiring an 'amp' signal for every
-- unit generator.  
module Euterpea.Audio.CSoundGenerators (
  returnA, 
  Table,
  oscil, oscili, oscils, oscil1, oscil1i,
  table, tablei, tableIx, tableiIx,
  buzz,
  pluck,
  PluckDecayMethod(..),
  delayt, delay, vdelay,
  comb,
  reson, areson,
  tone, atone,
  rand, randi, randh,
  balance,
  butterlp, butterhp, butterbp, butterbr,
  line, expon, linseg, expseg, linen, 
  envlpx,
  integral,
  -- GEN routines
  gen05, gen05', exponential1,
  gen07, gen07', lineSeg1,
  gen09, gen09', compSine2,
  gen10, gen10', compSine1,
  gen12, gen12'
) where

import Prelude hiding (init)

import Euterpea.Audio.Basics
import Euterpea.Audio.Types
import Control.Arrow hiding (returnA)
import Control.CCA.ArrowP
import Control.CCA.Types
import Control.CCA.CCNF
import Data.Array.Base (unsafeAt)
import Data.Array.Unboxed

import Debug.Trace

import Language.Haskell.TH
import Language.Haskell.TH.Instances
import Language.Haskell.TH.Syntax

import Foreign.Marshal
import Foreign.Marshal.Alloc
import Foreign.Ptr
import Foreign.Storable

import GHC.IO

import System.Random

-- type SigFun clk a b = ArrowP ASyn clk a b


wrap val bound = if val > bound then wrap val (val-bound) else val

clip val lower upper 
    | val <= lower = lower
    | val >= upper = upper
    | otherwise    = val

-- Raises 'a' to the power of 'b' using logarithms.
pow a b = exp (log a * b)

-- Returns the fractional part of 'x'.
frac x = if x > 1 then x - fromIntegral (truncate x) else x

data Table = Table 
    !Int                   -- size
    !(UArray Int Double)   -- table implementation
    ExpQ                   -- TH expression to construct 
                           -- the table at compile time
    !Bool                   -- Whether the table is normalized

instance Show Table where
    show (Table sz a _ n) = "Table with " ++ show sz ++ " entries: " ++ 
                            show a

instance Language.Haskell.TH.Syntax.Lift Table where
    lift (Table sz uarr fexp norm) =
         [| funToTable ($(fexp)) fexp norm sz |]

funToTable :: (Double->Double) -> ExpQ -> Bool -> Int -> Table
funToTable f f' normalize size = 
    let delta = 1 / fromIntegral size
        ys = take size (map f [0, delta.. ]) ++ [head ys]
             -- make table one size larger as an extended guard point
        zs = if normalize then map (/ maxabs ys) ys else ys
        maxabs = maximum . map abs
    in Table size (listArray (0, size) zs) f' normalize

readFromTable :: Table -> Double -> Double
readFromTable (Table sz array _ _) pos = 
    let idx = truncate (fromIntegral sz * pos)  -- range must be [0,size]
    in array `unsafeAt` idx
{-# INLINE [0] readFromTable #-}

readFromTableA t = arr' [| readFromTable t |] (readFromTable t)

readFromTableRaw :: Table -> Int -> Double
readFromTableRaw (Table _ a _ _) idx = a `unsafeAt` idx


-- Like readFromTable, but with linear interpolation
readFromTablei :: Table -> Double -> Double
readFromTablei (Table sz array _ _) pos = 
    let idx  = fromIntegral sz * pos  -- fractional "index" in table ([0,sz])
        idx0 = (truncate idx) `mod` sz       :: Int
        idx1 = idx0 + 1                      :: Int
        val0 = array `unsafeAt` idx0
        val1 = array `unsafeAt` idx1
    in val0 + (val1 - val0) * (idx - fromIntegral idx0)
{-# INLINE [0] readFromTablei #-}

readFromTableiA t = arr' [| readFromTablei t |] (readFromTablei t)

-- | Accesses table values by direct indexing with linear interpolation.
-- The index 'pos' is expected to be normalized (between 0 and 1).  Values
-- out of bounds are either clipped or wrapped.
tablei :: (Clock p, ArrowInit a) => 
          Table   -- ^ Table to read from.
       -> Bool    -- ^ Whether to wrap around index; 
                  --   if not, index is clipped within bounds
       -> ArrowP a p Double Double
tablei tab True =
    proc pos -> do 
      returnA -< readFromTablei tab (wrap pos 1)
tablei tab False =
    proc pos -> do
      returnA -< readFromTablei tab (clip pos 0 1)

-- | Accesses table values by direct indexing.
-- The index is normalized (between 0 and 1).
table :: (Clock p, ArrowInit a) => Table -> Bool -> ArrowP a p Double Double
table tab True =
    proc pos -> do 
      returnA -< readFromTable tab (wrap pos 1)
table tab False =
    proc pos -> do
      returnA -< readFromTable tab (clip pos 0 1)

-- | Like tablei, but the index is taken in as raw 
-- (between 0 and (size of table - 1), inclusive)
tableiIx :: (Clock p, ArrowInit a) => Table -> Bool -> ArrowP a p Double Double
tableiIx tab@(Table sz array _ _) True =
    proc idx -> do
      let idx0 = (truncate idx) `mod` sz
          val0 = readFromTableRaw tab idx0
          val1 = readFromTableRaw tab (idx0 + 1)
      returnA -< val0 + (val1 - val0) * (idx - fromIntegral idx0)
tableiIx tab@(Table sz _ _ _) False =
    proc idx -> do
      let pos = idx / fromIntegral (sz-1)
      returnA -< readFromTablei tab (clip pos 0 1)

-- | Like table, but index taken as raw.
tableIx :: (Clock p, ArrowInit a) => Table -> Bool -> ArrowP a p Double Double
tableIx tab@(Table sz array _ _) True =
    proc idx -> do
      returnA -< readFromTableRaw tab (truncate idx `mod` (sz-1))
tableIx tab@(Table sz array _ _) False =
    proc idx -> do
      returnA -< readFromTableRaw tab (clip (truncate idx) 0 (sz-1))

-- | 'oscil' generates periodic signals consisting of the values
-- returned from sampling a stored function table. The internal phase
-- is simultaneously advanced in accordance with the input signal
-- 'freq'.
oscil :: (Clock p, ArrowInit a) =>
         Table 
      -> Double  -- ^ Initial phase of sampling, expressed as a
                 -- fraction of a cycle (0 to 1).
      -> ArrowP a p Double Double
oscil table iphs = oscil_ iphs >>> readFromTableA table

-- | 'oscili' is like 'oscil', but with linear interpolation.
oscili :: (Clock p, ArrowInit a) => 
          Table 
       -> Double 
       -> ArrowP a p Double Double
oscili table iphs = oscil_ iphs >>> readFromTableiA table

-- helper function for oscil and oscili
oscil_ :: (Clock p, ArrowInit a) => Double -> ArrowP a p Double Double
oscil_ phs = 
    let sr = rate (undefined :: p)
    in proc freq -> do
      rec 
        let delta = 1 / sr * freq
            phase = if next > 1 then frac next else next
        next <- init phs -< frac (phase + delta)
      returnA -< phase

-- | Simple, fast sine oscillator, that uses only one multiply, and two
-- add operations to generate one sample of output, and does not
-- require a function table.
oscils :: (Clock p, ArrowInit a) =>
          Double -> ArrowP a p Double Double
oscils freq =
  let omh = 2 * pi * freq / sr
      d   = sin omh
      c   = 2 * cos omh
      sr  = rate (undefined :: p)
      sf  = proc amp -> do
               rec
                 let r = c * d2 - d1
                 d1 <- init 0         -< d2
                 d2 <- init d         -< r
               returnA -< r * amp
  in sf


-- | 'oscil1' accesses values by sampling once through the function
-- table at a rate determined by 'dur'. For the first 'del' seconds,
-- the point of scan will reside at the first location of the table;
-- it will then begin moving through the table at a constant rate,
-- reaching the end in another 'dur' seconds; from that time on
-- (i.e. after 'del' + 'dur' seconds) it will remain pointing at the
-- last location. Each value obtained from sampling is then multiplied
-- by an amplitude factor kamp before being written into the result.
oscil1 :: (Clock p, ArrowChoice a, ArrowInit a) =>
          Table
       -> Double
       -- ^ delay in seconds before 'oscil1' incremental sampling begins
       -> Double
       -- ^ duration in seconds to sample through the table just once.
       -> ArrowP a p Double Double
oscil1 = oscil1_ oscil

-- | Like 'oscil1', but with linear interpolation.
oscil1i :: (Clock p, ArrowChoice a, ArrowInit a) => 
           Table
        -> Double                  
           -- ^ delay in seconds before 'oscil1' incremental sampling begins.
        -> Double                  
           -- ^ duration in seconds to sample through the table just once.
        -> ArrowP a p Double Double
oscil1i = oscil1_ oscili

-- helper function for oscil1 and oscil1i
oscil1_ :: (Clock p, ArrowChoice a, ArrowInit a) => (Table -> Double -> ArrowP a p Double Double)
           -> Table -> Double -> Double -> ArrowP a p Double Double
oscil1_ osc table@(Table sz _ _ _) del dur =
  let sr = rate (undefined :: p)
      t1 = del * sr
      t2 = t1 + dur * sr
      v0 = readFromTableRaw table 0
      v2 = readFromTableRaw table (sz-1)
  in proc amp -> do
       let state i = if i < t1 then 0 else if i < t2 then 1 else 2
       i <- countUp -< ()
       y <- case state (fromIntegral i) of
              0 -> returnA        -< v0
              1 -> osc table 0  -< 1 / dur
              2 -> returnA        -< v2
       returnA -< y * amp

foscil, foscili :: (Clock p, ArrowInit a) => Table -> ArrowP a p (Double,Double,Double,Double) Double
foscil table =
    proc (freq,carfreq,modfreq,modindex) -> do
      returnA -< 0

foscili table =
    proc (freq,carfreq,modfreq,modindex) -> do
      returnA -< 0

loscil :: (Clock p, ArrowInit a) => Table -> ArrowP a p Double Double
loscil table = 
    proc freq -> do
      returnA -< 0


-- | Output a set of harmonically related sine partials. 
buzz :: Clock p => 
        Table   -- ^ table number of a stored function containing a
                -- sine wave. A large table of at least 8192 points is
                -- recommended.
     -> Double  -- ^ initial phase of the fundamental frequency,
                -- expressed as a fraction of a cycle (0 to 1).
     -> Signal p (Double,Int) Double -- ^ 'freq' is the base frequency
                                     -- in cycles per second, and
                                     -- 'nharms' is the total number
                                     -- of harmonics requested.
buzz table initialPhase =
    let sr = rate (undefined :: p)
    in proc (freq, nharms) -> do
      rec
        let delta = 1 / sr * freq
            phase = if next > 1 then frac next else next
        next <- init initialPhase -< frac (phase + delta)
      returnA -< sum [ readFromTable table (frac (phase * fromIntegral pn)) | pn <- [1..nharms] ] / fromIntegral nharms

instance Lift PluckDecayMethod where
    lift SimpleAveraging = [| SimpleAveraging |]
    lift (WeightedAveraging a b) = [| WeightedAveraging a b |]
    -- TODO: rest of the methods

data PluckDecayMethod
    = SimpleAveraging
      -- ^ A simple smoothing process.
    | StretchedAveraging Double  
      -- ^ with smoothing time stretched by a factor.
    | SimpleDrum Double
      -- ^ The range from pitch to noise is controlled by a 'roughness
      -- factor' (0 to 1). Zero gives the plucked string effect, while
      -- 1 reverses the polarity of every sample (octave down, odd
      -- harmonics). The setting .5 gives an optimum snare drum.
    | StretchedDrum Double Double  
      -- ^ Combines both roughness and stretch factors. parm1 is
      -- roughness (0 to 1), and parm2 the stretch factor (=1).
    | WeightedAveraging Double Double 
      -- ^ As SimpleAveraging, with parm1 weighting the current sample
      -- (the status quo) and iparm2 weighting the previous adjacent
      -- one. iparm1 + iparm2must be <= 1.
    | RecursiveFilter
      -- ^ 1st order recursive filter, with coefs .5. Unaffected by
      -- parameter values.
pluck :: Clock p => Table -> Double -> PluckDecayMethod -> Signal p Double Double
pluck table pitch method = 
    let sr = rate (undefined :: p) 
    in proc cps -> do
      rec 
        z <- delayt (max 64 (truncate (sr / pitch))) table -< y
        z' <- init 0 -< z
        let y = case method of 
                  SimpleAveraging -> 0.5 * (z + z') -- or is this "RecursiveFilter?"
                  WeightedAveraging a b -> z * a + z' * b
                  _ -> error "pluck: method not implemented"
      returnA -< y

grain :: Table 
         -- ^ Grain waveform. This can be just a sine wave or a sampled sound. 
      -> Table
         -- ^ Amplitude envelope used for the grains.
      -> Double
         -- ^ Maximum grain duration in seconds. This is the biggest
         -- value to be assigned to 'gdur'.
      -> Bool
         -- ^ If 'True', all grains will begin reading from the
         -- beginning of the 'gfn' table. If 'False', grains
         -- will start reading from random 'gfn' table positions.
      -> Signal p (Double,Double,Double,Double,Double) Double
grain gfn wfn mgdur grnd = 
    proc (pitch,dens,ampoff,pitchoff,gdur) -> do
        returnA -< 0


-- delayr and delayw are obsolete. 
-- One can use a fixed-time delay with native recursive arrow syntax to achieve
-- modified feedback loops.

updateBuf :: Buf -> Int -> Double -> IO Double
updateBuf (Buf _ a) i u = a `seq` i `seq` u `seq` do
    let p = a `advancePtr` i
    x' <- peek p
    poke p u
    return x'

peekBuf (Buf sz a) i = peek (a `advancePtr` (min (sz-1) i))

-- TODO: deal with pre-initialized buffers
instance Lift Buf where
    lift (Buf sz _) = [| mkArr sz |]

mkArr :: Int -> Buf
mkArr n = n `seq` Buf n (unsafePerformIO $ 
            Foreign.Marshal.newArray (replicate n 0))

mkArrWithTable size t = Buf size (unsafePerformIO $
    Foreign.Marshal.newArray (map (readFromTable t) [0, (1/sz)..((sz-1)/sz)]))
      where sz = fromIntegral size

data Buf = Buf !Int !(Ptr Double)

-- A fixed-length delay line, initialized using a table
delayt :: Clock p => Int -> Table -> Signal p Double Double
delayt size table =
    let sr = rate (undefined :: p)
        buf = mkArrWithTable size table
    in proc x -> do
        rec
          let i' = if i == size-1 then 0 else i+1
          i <- init 0 -< i'
          y <- init 0 -< x  
        -- TODO: this proc can't be strict on x, but how can we deal with strictness better without this hack?
        returnA -< unsafePerformIO $ updateBuf buf i y

-- | A fixed-length delay line.
delay :: Clock p => Double -> Signal p Double Double
delay maxdel =
    let sr = rate (undefined :: p)
        sz = truncate (sr * maxdel)
        buf = mkArr sz
    in proc x -> do
        rec
          let i' = if i == sz-1 then 0 else i+1
          i <- init 0 -< i'
          y <- init 0 -< x  
        returnA -< unsafePerformIO $ updateBuf buf i y

delay1 :: Clock p => Double -> Signal p (Double, Double) Double
delay1 = vdelay

-- | delay line with one tap.
vdelay :: Clock p => Double -> Signal p (Double, Double) Double
vdelay maxdel =
    let sr = rate (undefined :: p)
        sz = truncate (sr * maxdel)
        buf = mkArr sz
    in proc (sig,dlt) -> do
      rec
        let i' = if i == sz-1 then 0 else i+1
            dl = min maxdel dlt
            tap = i - truncate (sr * dl)
            tapidx = if tap < 0 then sz + tap else tap
        i <- init 0 -< i'
        y <- init 0 -< sig
      returnA -< unsafePerformIO $ do
        s <- peekBuf buf tapidx
        updateBuf buf i y
        return s

-- delay line with two taps.
delay2 :: Double -> Signal p (Double, Double, Double) Double
delay2 maxdel = 
    proc (sig, dlt1, dlt2) -> do
      returnA -< 0

-- delay line with three taps.
delay3 :: Double -> Signal p (Double, Double, Double, Double) Double
delay3 maxdel = 
    proc (sig, dlt1, dlt2, dlt3) -> do
      returnA -< 0

-- delay line with four taps.
delay4 :: Double -> Signal p (Double, Double, Double, Double, Double) Double
delay4 maxdel = 
    proc (sig, dlt1, dlt2, dlt3, dlt4) -> do
      returnA -< 0

instance Language.Haskell.TH.Syntax.Lift StdGen where
    lift g = [| g |]

-- | Generate uniform white noise with an R.M.S value of 'amp' / sqrt 2,
-- where 'seed' is the random seed.
rand :: Int -> Signal p Double Double
rand seed =
    let gen = mkStdGen seed
    in proc amp -> do
      rec
        let (a,g') = random g
        g <- init gen -< g'
      returnA -< (a * 2 - 1) * amp

-- | Generates a controlled random number series with interpolation
-- between each new number, where the internal pseudo-random formula
-- generate uniform white noise with an R.M.S value of 'amp' / sqrt 2,
-- 'seed' is the random seed.
randi :: Clock p => Int -> Signal p (Double, Double) Double
randi seed =
    let sr = rate (undefined :: p)
        gen = mkStdGen seed
        (i_n1, i_g1) = random gen
        (i_n2, i_g2) = random i_g1
        i_pr = (i_n1, i_n2, i_g2)
    in proc (amp, cps) -> do
      let bound = sr / cps
      rec
        state <- init (0, i_pr) -< state'
        let (cnt, pr@(n1, n2, g)) = state
            n = n1 + (n2 - n1) * cnt / bound
            state' = if cnt + 1 < bound then (cnt + 1, pr)
                                        else let (n3, g') = random g
                                             in (0, (n2, n3, g'))
      returnA -< (n * 2 - 1) * amp

-- | Generates a controlled random number series with interpolation
-- between each new number, where the internal pseudo-random formula
-- generate uniform white noise with an R.M.S value of 'amp' / sqrt 2,
-- 'seed' is the random seed.
randh :: Clock p => Int -> Signal p (Double, Double) Double
randh seed =
    let sr = rate (undefined :: p)
        gen = mkStdGen seed
        (i_n1, i_g) = random gen
        i_pr = (i_n1, i_g)
    in proc (amp, cps) -> do
      let bound = sr / cps
      rec
        state <- init (0, i_pr) -< state'
        let (cnt, pr@(n, g)) = state
            state' = if cnt + 1 < bound then (cnt + 1, pr)
                                        else let (n', g') = random g
                                             in (0, (n', g'))
      returnA -< (n * 2 - 1) * amp

balance :: Clock p =>
           Int
        -> Signal p (Double, Double) Double
balance ihp =
    proc (sig, ref) -> do
      rec
        (sqrsum, refsum) <- init (0, 0) -< (sqrsum', refsum')
        let sqrsum' = c1 * sig * sig + c2 * sqrsum
            refsum' = c1 * ref * ref + c2 * refsum
            ratio   = if sqrsum == 0 then sqrt $ refsum
                                     else sqrt $ refsum / sqrsum
      returnA -< sig * ratio

  where sr = rate (undefined :: p)
        tpidsr = 2 * pi / sr      -- tpidsr = two-pi over sr

        b  = 2 - cos (fromIntegral ihp * tpidsr)
        c1 = 1 - c2
        c2 = b - sqrt (b * b - 1)
        dup a = (a, a)

data ResonData = ResonData{ rsnKcf     :: !Double
                          , rsnKbw     :: !Double
                          , rsnCosf    :: !Double
                          , rsnC1      :: !Double
                          , rsnC2      :: !Double
                          , rsnC3      :: !Double
                          , rsnYt1     :: !Double
                          , rsnYt2     :: !Double
                          }
rsnDefault = ResonData (-1) (-1) 0 0 0 0 0 0

-- | A second-order resonant filter.
reson :: Clock p =>
         Int -- ^ 'scale': 1 signifies a peak response factor of 1, i.e. all
             -- frequencies other than kcf are attenuated in accordance with
             -- the (normalized) response curve; 2 raises the response
             -- factor so that its overall RMS value equals 1; 0 ignifies
             -- no scaling of the signal, leaving that to some later
             -- adjustment (like balance).
      -> Signal p (Double, Double, Double) Double
             -- ^ 'sig' is the signal to be filtered,
             -- 'kcf' is the center frequency of the filter,
             -- and 'kbw' is the bandwidth of it.
reson scale =
    proc (sig, kcf, kbw) -> do
      rec
        rsnData  <- init rsnDefault -< rsnData'
        currData <- if kcf == rsnKcf rsnData && kbw == rsnKbw rsnData
                         then returnA -<  rsnData
                         else update  -< (rsnData, kcf, kbw)
        let ResonData{ rsnC1 = c1, rsnC2 = c2, rsnC3 = c3,
                       rsnYt1 = yt1, rsnYt2 = yt2 } = currData
            a = c1 * sig + c2 * yt1 - c3 * yt2
            rsnData' = currData{ rsnYt1 = a, rsnYt2 = yt1 }
      returnA -< a

  where sr = rate (undefined :: p)
        tpidsr = 2 * pi / sr      -- tpidsr = two-pi over sr
        update = proc (rsnData, kcf, kbw) -> do
          -- kcf or kbw changed, recalc consts
          let cosf = cos $ kcf * tpidsr   -- cos (2pi * freq / rate)
              c3   = exp $ - kbw * tpidsr -- exp (-2pi * bwidth / rate)
                              -- (note on csound code) mtpdsr = -tpidsr

              -- c1   Gain for input signal.
              -- c2   (Minused) gain for output of delay 1.
              -- c3   Gain for output of delay 2.

              c3p1 = c3 + 1
              c3t4 = c3 * 4
              c2   = c3t4 * cosf / c3p1

              omc3  = 1 - c3
              c2sqr = c2 * c2
              c1 = case scale of
                1 -> omc3 * sqrt (1 - c2sqr / c3t4)
                2 -> sqrt $ (c3p1 * c3p1 - c2sqr) * omc3 / c3p1
                _ -> 1.0
          returnA -< rsnData{ rsnKcf = kcf, rsnKbw = kbw, rsnCosf = cosf,
                              rsnC1  = c1,  rsnC2  = c2,  rsnC3   = c3   }

-- | A notch filter whose transfer functions are the complements
-- of the reson opcode.
areson scale = proc (sig, kcf, kbw) -> do
  r <- reson scale -< (sig, kcf, kbw)
  returnA -< sig - r

data ButterData = ButterData !Double !Double !Double !Double !Double

sqrt2 = sqrt 2

blpset freq sr = ButterData a1 a2 a3 a4 a5
  where c = 1 / tan (pidsr * freq)
        csq = c * c; pidsr = pi / sr
        a1 = 1 / (1 + sqrt2 * c + csq)
        a2 = 2 * a1
        a3 = a1
        a4 = 2 * (1 - csq) * a1
        a5 = (1 - sqrt2 * c + csq) * a1
        
bhpset freq sr = ButterData a1 a2 a3 a4 a5
  where c = tan (pidsr * freq)
        csq = c * c; pidsr = pi / sr
        a1 = 1 / (1 + sqrt2 * c + csq)
        a2 = (-2) * a1
        a3 = a1
        a4 = 2 * (csq - 1) * a1
        a5 = (1 - sqrt2 * c + csq) * a1
        
bbpset freq band sr = ButterData a1 a2 a3 a4 a5
  where c = 1 / tan (pidsr * band)
        d = 2 * cos (2 * pidsr * freq)
        pidsr = pi / sr
        a1 = 1 / (1 + c)
        a2 = 0
        a3 = negate a1
        a4 = negate (c * d * a1)
        a5 = (c - 1) * a1
        
bbrset freq band sr = ButterData a1 a2 a3 a4 a5
  where c = tan (pidsr * band)
        d = 2 * cos (2 * pidsr * freq)
        pidsr = pi / sr
        a1 = 1 / (1 + c)
        a2 = negate d * a1
        a3 = a1
        a4 = a2
        a5 = (1 - c) * a1

-- | Implementation of a second-order low-pass Butterworth filter,
-- where 'sig' is the input signal to be filtered, and
-- 'freq' is the cutoff center frequency.
butterlp :: Clock p => Signal p (Double, Double) Double
butterlp = 
  let sr = rate (undefined :: p) 
  in proc (sig, freq) -> do
    butter -< (sig, blpset freq sr)

-- | A high-pass Butterworth filter.
butterhp :: Clock p => Signal p (Double, Double) Double
butterhp = 
  let sr = rate (undefined :: p)
  in proc (sig, freq) -> do
    butter -< (sig, bhpset freq sr)

-- | A band-pass Butterworth filter where 'band' is the bandwidth.
-- 'butterbr -< (s, 2000, 100)' will pass only 1950 to 2050 Hz in 's'.
butterbp :: Clock p => Signal p (Double, Double, Double) Double
butterbp = 
  let sr = rate (undefined :: p)
  in proc (sig, freq, band) -> do
    butter -< (sig, bbpset freq band sr)

-- | A band-reject Butterworth filter where 'band' is the bandwidth.
-- 'butterbr -< (s, 4000, 1000)' will filter 's' such that frequencies 
-- between 3500 to 4500 Hz are rejected.
butterbr :: Clock p => Signal p (Double, Double, Double) Double
butterbr = 
  let sr = rate (undefined :: p)
  in proc (sig, freq, band) -> do
    butter -< (sig, bbrset freq band sr)

-- Helper function for various Butterworth filters.
butter :: Clock p => Signal p (Double,ButterData) Double
butter = proc (sig, ButterData a1 a2 a3 a4 a5) -> do
    rec let t = sig - a4 * y' - a5 * y''
            y = t * a1 + a2 * y' + a3 * y''
        y'  <- init 0 -< t
        y'' <- init 0 -< y'
    returnA -< y


-- | This filter reiterates input with an echo density determined by
-- loop time 'looptime'. The attenuation rate is independent and is
-- determined by 'rvt', the reverberation time (defined as the time in
-- seconds for a signal to decay to 1/1000, or 60dB down from its
-- original amplitude). Output from a 'comb' filter will appear only
-- after 'looptime' seconds.
comb :: Clock p => 
        Double -- ^ loop time in seconds, which determines the “echo
               -- density” of the reverberation. This in turn
               -- characterizes the “color” of the comb filter whose
               -- frequency response curve will contain 'looptime' *
               -- sr/2 peaks spaced evenly between 0 and sr/2 (the
               -- Nyquist frequency). Loop time can be as large as
               -- available memory will permit.
     -> Signal p (Double, Double) Double
comb looptime = 
    let log001 = -6.9078
        del = delay looptime
    in proc (sig, rvt) -> do
      let gain = exp (log001 * looptime / rvt)
      rec
        r <- del -< sig + r * gain
      returnA -< r

-- | A first-order recursive low-pass filter with variable frequency
-- response. 'hp' is the response curve's half-power point, in
-- Hertz. Half power is defined as peak power / sqrt 2.
tone :: Clock p => Signal p (Double,Double) Double
tone = 
    let sr = rate (undefined :: p)
    in proc (sig, hp) -> do
        rec
           let y' = c1 * sig + c2 * y
               b = 2 - cos (2 * pi * hp / sr)
               c2 = b - sqrt (b * b - 1.0)
               c1 = 1 - c2
           y <- init 0 -< y'
        returnA -< y

-- | A hi-pass filter whose transfer functions are the complements of
-- the 'tone' opcode.  'atone' is thus a form of high-pass filter
-- whose transfer functions represent the “filtered out” aspects of
-- their complements. However, power scaling is not normalized in
-- 'atone' but remains the true complement of the corresponding
-- unit. Thus an audio signal, filtered by parallel matching 'tone'
-- and 'atone' units, would under addition simply reconstruct the
-- original spectrum.
atone :: Clock p => Signal p (Double,Double) Double
atone = proc (sig, hp) -> do
       y <- tone -< (sig, hp)
       returnA -< sig - y

--  | 'line' generates control or audio signals whose values move
--  linearly from an initial value to a final one.  A common error
--  with this opcode is to assume that the value of 'b' is held after
--  the time 'dur'. 'line' does not automatically end or stop at the
--  end of the duration given. If your note length is longer than
--  'dur' seconds, the resulting value will not come to rest at 'b',
--  but will instead continue to rise or fall with the same rate. If a
--  rise (or fall) and then hold is required that the 'linseg' opcode
--  should be considered instead.
line :: Clock p => 
        Double  -- ^ Starting value.
     -> Double  -- ^ Duration in seconds.
     -> Double  -- ^ Value after 'dur' seconds.
     -> Signal p Double Double
line a dur b =
    let sr = rate (undefined :: p)
    in proc sig -> do
      rec
        y <- init a -< y + (b-a) * (1 / sr / dur)
      returnA -< y * sig



-- | Trace an exponential curve between specified points. 
expon :: Clock p => 
        Double  -- ^ Starting value.  Zero is illegal for exponentials. 
     -> Double  -- ^ Duration in seconds.                
     -> Double  -- ^ Value after 'dur' seconds.  For exponentials,
                -- must be non-zero and must agree in sign with 'a'.
     -> Signal p Double Double
expon a dur b =
    let sr = rate (undefined :: p)
    in proc sig -> do
      rec
        y <- init a -< y * pow (b/a) (1 / sr / dur)
      returnA -< y * sig

-- Unfortunately, line and expon cannot be abstracted to a common
-- function because Template Haskell doesn't like higher-order functions.

data Tab = Tab [Double] !Int !(UArray Int Double)

instance Language.Haskell.TH.Syntax.Lift Tab where
    lift (Tab xs sz uarr) =
        [| Tab xs sz (listArray (0, sz-1) xs) |]

aAt (Tab _ sz a) i = unsafeAt a (min (sz-1) i)



-- Helper function for linseg and expseg.
seghlp :: Clock p =>
           [Double]  -- ^ List of points to trace through.
        -> [Double]  -- ^ List of durations for each line segment.
                     -- Needs to be one element fewer than 'iamps'.
        -> Signal p () (Double,Double,Double,Double)
seghlp iamps idurs =
    let sr = rate (undefined :: p)
        sz = length iamps
        amps = Tab iamps sz (listArray (0, sz-1) iamps)
        durs = Tab idurs (sz-1) (listArray (0, sz-2) (map (*sr) idurs))
    in proc _ -> do
      -- TODO: this is better defined using 'integral', but which is faster?
      rec
        let (t', i') = if t >= durs `aAt` i 
                       then if i == sz-2 then (t+1, i) else (0, i+1)
                       else (t+1, i)
        i <- init 0 -< i'
        t <- init 0 -< t'
      let a1 = aAt amps i
          a2 = aAt amps (i+1)
          d  = aAt durs i
      returnA -< (a1,a2,t,d)


-- | Trace a series of line segments between specified points.
linseg :: Clock p => 
           [Double]  -- ^ List of points to trace through.
        -> [Double]  -- ^ List of durations for each line segment.
                     -- Needs to be one element fewer than 'amps'.
        -> Signal p () Double
linseg amps durs = 
    let sf = seghlp amps durs
    in proc () -> do
      (a1,a2,t,d) <- sf -< ()
      returnA -< a1 + (a2-a1) * (t / d)


-- | Trace a series of exponential segments between specified points.
expseg :: Clock p => 
           [Double]  -- ^ List of points to trace through.
        -> [Double]  -- ^ List of durations for each line segment.
                     -- Needs to be one element fewer than 'amps'.
        -> Signal p () Double
expseg (a:amps) durs = 
    let amps' = max 0.001 a : amps 
        sf = seghlp amps' durs
    in proc () -> do
      (a1,a2,t,d) <- sf -< ()
      returnA -< a1 * pow (a2/a1) (t / d)

-- | Applies a straight line rise and decay pattern to an input amp
-- signal.  Rise modifications are applied for the first 'rise'
-- seconds, and decay from time 'dur' - 'dec'. If these periods are
-- separated in time there will be a steady state during which the
-- input signal will be unmodified. If the overall duration idur is
-- exceeded in performance, the final decay will continue on in the
-- same direction, going negative.
linen :: (Clock p) =>
         Double  -- ^ rise time in seconds.
      -> Double  -- ^ overall duration in seconds.
      -> Double  -- ^ decay time in seconds.
      -> Signal p Double Double
linen rise dur dec = 
    let sf = linseg [0,1,1,0] [rise, dur-rise-dec, dec]
    in proc sig -> do
      env <- sf -< ()
      returnA -< sig * env

-- | Apply an envelope consisting of 3 segments:
--   1. stored function rise shape
--   2. modified exponential pseudo steady state
--   3. exponential decay
--
--   Rise modifications are applied for the first 'rise' seconds, and
--   decay from time 'dur' - 'dec'. If these periods are separated in 
--   time there will be a steady state during which amp will be 
--   modified by the first exponential pattern. If rise and decay 
--   periods overlap then both modifications will be in effect for 
--   that time. If the overall duration 'dur' is exceeded in 
--   performance, the final decay will continue on in the same 
--   direction, tending asymptotically to zero.
envlpx :: Clock p =>
          Double  -- ^ rise time in seconds.
       -> Double  -- ^ overall duration in seconds.
       -> Double  -- ^ decay time in seconds.
       -> Table   -- ^ table of stored rise shape.
       -> Double  
          -- ^ attenuation factor, by which the last value of the
          --   'envlpx' rise is modified during the note's pseudo
          --   steady state. A factor greater than 1 causes an
          --   exponential growth and a factor less than 1 creates an
          --   exponential decay. A factor of 1 will maintain a true
          --   steady state at the last rise value. Note that this
          --   attenuation is not by fixed rate (as in a piano), but
          --   is sensitive to a note's duration. However, if 'atss'
          --   is negative (or if steady state < 4 k-periods) a fixed
          --   attenuation rate of 'abs' 'atss' per second will be
          --   used. 0 is illegal.
       -> Double  
          -- ^ attenuation factor by which the closing steady state
          --   value is reduced exponentially over the decay
          --   period. This value must be positive and is normally of
          --   the order of .01. A large or excessively small value is
          --   apt to produce a cutoff which is audible. A zero or
          --   negative value is illegal.
       -> Signal p Double Double
envlpx rise dur dec tab atss atdec = 
    let sr     = rate (undefined :: p)
        cnt1   = (dur - rise - dec) * sr + 0.5 
                 -- num of samples in steady state
        mlt1   = pow atss  (1 / cnt1)
        mlt2   = pow atdec (1 / sr / dec)
    in proc amp -> do
      rec 
        i <- countUp -< ()
        let state i
              | i < rise * sr      = 0
              | i < (dur-dec) * sr = 1
              | otherwise          = 2
        y  <- init (readFromTableRaw tab 0) -< y'
        y' <- case state (fromIntegral i) of 
                 0 -> table tab False -< fromIntegral i / (rise*sr+0.5)
                 1 -> returnA -< y * mlt1
                 2 -> returnA -< y * mlt2
      returnA -< amp * y'

-------------------------------
-- GEN routines
-------------------------------
type TableSize = Int
type PartialNum = Double
type PartialStrength = Double
type PhaseOffset = Double

type StartPt = Double
type SegLength = Double
type EndPt = Double

type DoubleSegFun = (Double, StartPt) -> [(SegLength, EndPt)] -> Double -> Double

gen05 :: TableSize
         -- ^ The size of the table to be produced. 
     ->  StartPt
         -- ^ The y-coordinate for the start point, (0,y). 
     -> [(SegLength, EndPt)]
         -- ^ Pairs of segment lengths and y-coordinates. The segment
         -- lengths are the projection along the x-axis. The first
         -- pair will define the line from (0, startPt) to (segLength,
         -- endPt).
     -> Table
gen05  size sp segs = gen05_ sp segs True size
gen05' size sp segs = gen05_ sp segs False size
gen05_ sp segs = funToTable (interpLine sp segs interpExpLine) 
                            [| interpLine sp segs interpExpLine |]
exponential1 = gen05

gen07 :: TableSize
         -- ^ The size of the table to be produced. 
     ->  StartPt
         -- ^ The y-coordinate for the start point, (0,y). 
     -> [(SegLength, EndPt)]
         -- ^ Pairs of segment lengths and y-coordinates. The segment
         -- lengths are the projection along the x-axis. The first
         -- pair will define the line from (0, startPt) to (segLength,
         -- endPt).
     -> Table
gen07  size sp segs = gen07_ sp segs True size
gen07' size sp segs = gen07_ sp segs False size
gen07_ sp segs = funToTable (interpLine sp segs interpStraightLine) 
                            [| interpLine sp segs interpStraightLine |]
lineSeg1 = gen07

-- | Make a table from a collection of sine waves at different offsets
-- and strengths.
gen09 :: TableSize
         -- ^ The size of the table to be produced.
      -> [(PartialNum, PartialStrength, PhaseOffset)]
         -- ^ List of triples of the partial (0,1,...), partial
         -- strength on [0,1], and phase offset on [0,360].
      -> Table
gen09  size ps = gen09_ ps True size
gen09' size ps = gen09_ ps False size
gen09_ ps = funToTable (makeCompositeSineFun ps) 
                       [| makeCompositeSineFun ps |]
compSine2 = gen09

gen10f pss x = let phase = 2 * pi * x 
               in sum (zipWith (*) [ sin (phase * pn) | pn <- [1..] ] pss)

gen10 :: TableSize -> [PartialStrength] -> Table
gen10  size pss = gen10_ pss True size
gen10' size pss = gen10_ pss False size
gen10_ pss = funToTable (gen10f pss) [| gen10f pss |]

compSine1 = gen10

-- | Generates the log of a modified Bessel function of the second
-- kind, order 0, suitable for use in amplitude-modulated FM.
gen12 :: TableSize 
      -> Double  -- ^ specifies the x interval [0 to +xint] over which
                 -- the function is defined.
      -> Table
gen12  size xint = gen12_ xint True size
gen12' size xint = gen12_ xint False size
gen12_ xint = funToTable (gen12f xint) [| gen12f xint |]
gen12f xint x =
    log $ 1 +
        let tsquare = x * x * xint * xint / 3.75 / 3.75
        in sum $ zipWith (*) [ 3.5156229, 3.0899424, 1.2067492,
                               0.2659732, 0.0360768, 0.0045813 ]
               $ iterate (*tsquare) tsquare


-- Begin utility functions for gen05 & gen07--

normalizeSegs :: [(SegLength, entPt)] -> [(SegLength, entPt)]
normalizeSegs segs =
    let s = sum (map fst segs)
        fact = if (s > 1) then (1/s) else 1 -- don't force max<1 up to max=1
    in  map (\(x,y) -> (x*fact, y)) segs


interpLine :: StartPt
              -- ^ The y-coordinate for the start point (0,y).
           -> [(SegLength, EndPt)]
              -- ^ Pairs of segment lengths (projected on the x-axis)
              -- and y-coordinates (end points).
           -> DoubleSegFun
              -- ^ The function to use for interpolation
           -> Double
              -- ^ The x-coordinate for which to find the
              -- corresponding f(x)=y.
           -> Double
interpLine sp [] d f = 0 -- catchall case
interpLine sp points f d = f (0,sp) (normalizeSegs points) d 


{-
The exponential interpolation function stretches e^x between two 
endpoints for each pair of points.
-}
interpExpLine :: (Double, StartPt)
                 -- ^ The startpoing as (x,y)
              -> [(SegLength, EndPt)]
                 -- ^ A list of line segments with (x',y) where x' is
                 -- a length projected on the x-axis
              -> Double
                 -- ^ The target x-coordinate to find a corresponding
                 -- y value for
              -> Double
interpExpLine (s1, e1) [] d = e1 -- termination case, end of list
interpExpLine (s1, e1) ((s2, e2):t) d = 
    if d > s2 then interpExpLine (s2, e2) t (d-s2) else
    let  h = e2 - e1 
         x = if h<0 then s2-d else d
    in   if s2<=0 then e2 else -- accomodate discontinuities
         (abs h)*((exp (x/s2))-1)/((exp 1)-1) + (min e1 e2)


interpStraightLine :: (Double, StartPt)
                 -- ^ The startpoing as (x,y)
              -> [(SegLength, EndPt)]
                 -- ^ A list of line segments with (x',y) where x' is
                 -- a length projected on the x-axis
              -> Double
                 -- ^ The target x-coordinate to find a corresponding
                 -- y value for
              -> Double
interpStraightLine (s1, e1) [] d = e1 -- termination case, end of list
interpStraightLine (s1, e1) ((s2, e2):t) d = 
    if d > s2 then interpStraightLine (s2, e2) t (d-s2) else
    let  h = e2 - e1 -- height of triangle
         s = h/s2 -- slope of triangle
    in   if s2<=0 then e2 else 
         e1 + (s*d) -- start point plus slope times distance



-- Function to find a particular point at a particular strength

makeSineFun :: (PartialNum, PartialStrength, PhaseOffset)
                -- ^ Triple of the partial (0,1,...), partial strength
                -- on [0,1], and phase offset on [0,360].
             -> Double
                -- ^ The x coordinate for which to find f(x)=y
             -> Double
makeSineFun (pNum, pStrength, pOffset) x = 
    let x' = x * 2 * pi -- convert [0,1] to [0,pi] radians
        po = (pOffset/360) * 2 * pi -- convert [0,360] to [0,pi] radians
    in  pStrength * sin (x' * pNum + po)



-- For a particular point, sum all partials

makeCompositeSineFun :: [(PartialNum, PartialStrength, PhaseOffset)]
                        -- ^ List of triples of the partial (0,1,...),
                        -- partial strength on [0,1], and phase offset
                        -- on [0,360].
                     -> Double
                        -- ^ The x coordinate for which to find f(x)=y
                     -> Double
makeCompositeSineFun [] x = 0
makeCompositeSineFun (p:ps) x = makeSineFun p x + makeCompositeSineFun ps x