signal.clj

(ns clojure2d.extra.signal
  "Signal processing and generation.
  
  ## Signal

  Signal is array of doubles from range `[-1.0.1.0]` packed into `Signal` type.

  Signal can be:

  * obtained from `Pixels` with [[pixels->signal]] function.
  * generated by calling [[signal-from-wave]].
  * loaded from file with [[load-signal]] function. Signal should be encoded as 16-bit signed integer big endian.

  ### Pixels as Signal

  `Pixels` (`Image`) can be treated as `Signal`. Conversion to Signal is based on strategies of converting image to RAW and then converting to audio. It includes channel data layout and packing into integer, encoding, endianess, etc.

  To convert `Pixels` to `Signal` use [[pixels->signal]] function. To convert back use [[signal->pixels]]. [[signal->pixels]] requires target `Pixels` object to store result of conversion. Target is mutated then.

  To filter `Pixels` directly (without explicit conversion to and from Signals) you can use [[filter-channels]] with [[effects-filter]].
  
  ### Wave as Signal

  You can create [[wave]] function from oscillator. You can also sum waves with [[sum-waves]].

  To sample wave to signal, call [[wave->signal]] with following parameters:

  * `f` - wave function
  * `samplerate` - sample rate (samples per second)
  * `seconds` - how many seconds generate
  
  ### File operations

  You can [[save-signal]] or [[load-signal]]. Representation is 16 bit signed, big endian. Use Audacity or SoX to convert to/from audio files.

  ## Signal processing

  To process Signal use [[apply-effects]] function on it.

  Effect is signal filter, created with [[effect]] multimethod. Effects can be composed with [[compose-effects]]. Effect can be treated as function and can be called for given sample.

  Each effect has it's own parametrization which should be passed during creation.
  
  List of all available effects is under [[effects-list]] value.

  ### Effects parametrization

  Each effect has its own parametrization

  #### :simple-lowpass, :simple-highpass

  * `:rate` - sample rate (default 44100.0)
  * `:cutoff` - cutoff frequency (default 2000.0)

  #### :biquad-eq

  Biquad equalizer
  
  * `:fc` - center frequency
  * `:gain` - gain
  * `:bw` - bandwidth (default: 1.0)
  * `:fs` - sampling rate (defatult: 44100.0)

  #### :biquad-hs, :biquad-ls

  Biquad highpass and lowpass shelf filters
  
  * `:fc` - center frequency
  * `:gain` - gain
  * `:slope` - shelf slope (default 1.5)
  * `:fs` - sampling rate (default 44100.0)

  #### :biquad-lp, :biquad-hp, :biquad-bp

  Biquad lowpass, highpass and bandpass filters
  
  * `:fc` - cutoff/center frequency
  * `:bw` - bandwidth (default 1.0)
  * `:fs` - sampling rate (default 44100.0)

  #### :dj-eq

  * `:high` - high frequency gain (10000Hz)
  * `:mid` - mid frequency gain (1000Hz)
  * `:low` - low frequency gain (100Hz)
  * `:shelf-slope` - shelf slope for high frequency (default 1.5)
  * `:peak-bw` - peak bandwidth for mid and low frequencies (default 1.0)
  * `:rate` - sampling rate (default 44100.0)

  #### :phaser-allpass

  * `:delay` - delay factor (default: 0.5)

  #### :divider

  * `:denom` (long, default 2.0)

  #### :fm

  Modulate and demodulate signal using frequency

   * `:quant` - quantization value (0.0 - if no quantization, default 10)
   * `:omega` - carrier factor (default 0.014)
   * `:phase` - deviation factor (default 0.00822)

  #### :bandwidth-limit

  https://searchcode.com/file/18573523/cmt/src/lofi.cpp#

  * `:rate` - sample rate (default 44100.0)
  * `:freq` - cutoff frequency (default 1000.0)

  #### :distort

  * `:factor` - distortion factor (default 1.0)

  #### :foverdrive

  Fast overdrive
  
  * `:drive` - drive (default 2.0)

  #### :decimator

  * `:bits` - bit depth (default 2)
  * `:fs` - decimator sample rate (default 4410.0)
  * `:rate` - input sample rate (default 44100.0)

  #### :basstreble

  * `:bass` - bass gain (default 1.0)
  * `:treble` - treble gain (default 1.0)
  * `:gain` - gain (default 0.0)
  * `:rate` - sample rate (default 44100.0)
  * `:slope` - slope for both (default 0.4)
  * `:bass-freq` - bass freq (default 250.0)
  * `:treble-freq` - treble freq (default 4000.0)

  #### :echo

  * `:delay` - delay time in seconds (default 0.5)
  * `:decay` - decay (amount echo in signal, default 0.5)
  * `:rate` - sample rate (default 44100.0)
  
  _Warning! Echo filter uses mutable array as a internal state, don't use the same filter in paraller processing._

  #### :vcf303

  * `:rate` - sample rate (default 44100.0)
  * `:trigger` - boolean, trigger some action (default `false`), set true when you reset filter every line
  * `:cutoff` - cutoff frequency (values 0-1, default 0.8)
  * `:resonance` - resonance (values 0-1, default 0.8)
  * `:env-mod` - envelope modulation (values 0-1, default 0.5)
  * `:decay` - decay (values 0-1, default 1.0)
  * `:gain` - gain output signal (default: 1.0)
  
  #### :slew-limit

  http://git.drobilla.net/cgit.cgi/omins.lv2.git/tree/src/slew_limiter.c

  * `:rate` - sample rate
  * `:maxrise` - maximum change for rising signal (in terms of 1/rate steps, default 500)
  * `:maxfall` - maximum change for falling singal (default 500)

  #### :mda-thru-zero

  * `:rate` - sample rate
  * `:speed` - effect rate
  * `:depth`
  * `:mix`
  * `:depth-mod`
  * `:feedback`
  
  _Warning: like `:echo` internal state is kept in doubles array._"
  {:metadoc/categories {:eff "Effects"
                        :wave "Wave"
                        :sig "Signal"}}
  (:require [fastmath.core :as m]
            [clojure2d.pixels :as p]
            [clojure2d.core :refer :all]
            [clojure.java.io :refer :all]
            [fastmath.vector :as v]
            [fastmath.random :as r])
  (:import [clojure2d.pixels Pixels]
           [fastmath.vector Vec2 Vec3]
           [clojure.lang Counted IFn Seqable Sequential]
           [clojure2d.java.signal Converter]))

(set! *warn-on-reflection* true)
(set! *unchecked-math* :warn-on-boxed)
(m/use-primitive-operators)

(deftype Signal [^doubles signal]
  Object
  (toString [_] (str "size=" (alength signal)))
  Counted
  (count [_] (alength signal))
  Sequential
  Seqable
  (seq [_] (seq signal)) )

(defn signal
  "Signal creator"
  {:metadoc/categories #{:sig}}
  [xs]
  (Signal. (double-array xs)))

(defn- coding-value
  "Encode `coding` as Converter class value."
  [coding]
  (condp = coding
    :alaw Converter/CODING_ALAW
    :ulaw Converter/CODING_ULAW
    :alaw-rev Converter/CODING_ALAW_REV
    :ulaw-rev Converter/CODING_ULAW_REV
    0))

(defn pixels->signal
  "Convert `Pixels` to `Signal` for given specifications.

  Specification for conversion is based on Audacity importing options and image to RAW conversions practices. Spec contains:

  * `planar?` - layout of channels data, planar (true, default) or interleaved (false)
  * `signed?` - channel data are encoded as signed (true, default) or unsigned (false) integer.
  * `little-endian?` - for 16 and 24 bits pack channel data with little-endian order (true, default) or big-endian (false)
  * `bits` - pack channel data in 8 (default), 16 or 24 bits.
  * `channels` - list of channels to convert to signal, default `[0 1 2]` (for all you can use `:all` keyword).
  * `coding` - encode signal with one of the encoders: `:alaw`, `:alaw-rev`, `:ulaw`, `:ulaw-rev` or `:none` (default).

  Pixels can be restored from Signal with [[signal->pixels]] function."
  {:metadoc/categories #{:sig}}
  ([^Pixels p {:keys [planar? signed? little-endian? ^int bits channels coding]
               :or {planar? true little-endian? true bits 8 signed? false channels [0 1 2] coding :none}}]
   (let [channels (int-array (if (= :all channels) [0 1 2 3] channels))]
     (Signal. (Converter/toSignal
               (.p p) channels
               (dec (>> bits 3))
               little-endian? signed? planar?
               (coding-value coding)))))
  ([p]
   (pixels->signal p {})))

(defn signal->pixels
  "Convert `Signal` to `Pixels` storing result into `target` (mutating it!).

  Specification is the same as in [[pixels-signal]] functions."
  {:metadoc/categories #{:sig}}
  ([^Signal sig ^Pixels target {:keys [planar? signed? little-endian? ^int bits channels coding]
                                :or {planar? true little-endian? true bits 8 signed? false channels [0 1 2] coding :none}}]
   (let [channels (int-array (if (= :all channels) [0 1 2 3] channels))]
     (Converter/fromSignal
      (.signal sig)
      (.p target)
      channels
      (dec (>> bits 3))
      little-endian? signed? planar?
      (coding-value coding))
     target))
  ([sig target]
   (signal->pixels sig target {})))


;; ## Signal processing
;;
;; Looks like a fastest implementation (reduce and vector types are 8-10 times slower).

(declare single-pass)

(deftype SampleAndState [^double sample state]
  Object
  (toString [_] (str sample))) ; sample and effect state, StateWithF or vector of StateWithF

;; Type representing list node consisting resulting sample, effect functions, state value and link to next node (or nil if last node)
(deftype EffectsList [effect-name ^double sample effect state next]
  Object
  (toString [_] (str (if next 
                       (str next " -> " (name effect-name))
                       (name effect-name)) " (" sample ")"))
  IFn
  (invoke [_] sample)
  (invoke [e n]
    (single-pass e n)))

(defn- effect-node
  "Create `EffectsList` node from effect function and initial state."
  [effect-name f]
  (EffectsList. effect-name 0.0 f (f) nil))

(defn- compose-two-effects
  "Join two effect lists."
  [^EffectsList e1 ^EffectsList e2]
  (EffectsList. (.effect-name e1) (.sample e1) (.effect e1) (.state e1) (if (.next e1)
                                                                          (compose-two-effects (.next e1) e2)
                                                                          e2)))

(defn compose-effects
  "Compose effects."
  {:metadoc/categories #{:eff}}
  [^EffectsList e & es]
  (reduce compose-two-effects e es))

(defn- reset-effects
  "Resets effects state to initial one."
  [^EffectsList e]
  (EffectsList. (.effect-name e) 0.0 (.effect e) ((.effect e)) (when (.next e)
                                                                 (reset-effects (.next e)))))

(defn- single-pass
  "Process on sample using effects, returns `EffectsList` with result and new effect states."
  [^EffectsList e ^double sample]
  (if-not (.next e)
    (let [^SampleAndState r ((.effect e) sample (.state e))]
      (EffectsList. (.effect-name e) (.sample r) (.effect e) (.state r) nil))
    (let [^EffectsList prev (single-pass (.next e) sample)]
      (let [^SampleAndState r ((.effect e) (.sample prev) (.state e))]
        (EffectsList. (.effect-name e) (.sample r) (.effect e) (.state r) prev)))))

(defn apply-effects
  "Apply effects to signal.

  If `reset` is positive, reinit state each `reset` number of samples.

  Returns new signal."
  {:metadoc/categories #{:eff}}
  ([^Signal s effects ^long reset]
   (let [len (count s)
         in (.signal s)
         out (double-array len)]
     (loop [idx (int 0)
            effects-and-state effects]
       (when (< idx len)
         (let [sample (aget ^doubles in idx)
               ^EffectsList res (single-pass effects-and-state sample)
               idx+ (inc idx)] 
           (aset ^doubles out idx ^double (.sample res))
           (recur idx+
                  (if (and (pos? reset) (zero? ^long (mod idx+ reset)))
                    (reset-effects effects)
                    res)))))
     (Signal. out)))
  ([sig effects] (apply-effects sig effects 0)))

(defn- process-pixels
  "Process pixels as signal based on config for given pixels p. Store result in given target t."
  [p t effects config config-back reset]
  (signal->pixels (apply-effects (pixels->signal p config) effects reset) t config-back))

(defn effects-filter
  "Creates filter to process `Pixels` as `Signal` using [[filter-channels]].

  Provide configurations to properly convert Pixels to Signal and back.
  Optionally set `reset` value (reinit effects' state after `reset` samples).

  Filter operates on one channel at time and is defined to be used with [[filter-channels]].
  If you want to operate on `Pixels` directly, use [[apply-effects-to-pixels]]."
  {:metadoc/categories #{:eff}}
  ([effects config config-back reset]
   (fn [ch target p]
     (let [c {:channels [ch]}]
       (process-pixels p target effects (merge config c) (merge config-back c) reset))))
  ([effects]
   (effects-filter effects {} {} 0))
  ([effects reset]
   (effects-filter effects {} {} reset))
  ([effects config config-back]
   (effects-filter effects config config-back 0)))

(defn apply-effects-to-pixels
  "Apply effects directly to `Pixels`.

  Provide configurations to properly convert Pixels to Signal and back.
  Optionally set `reset` value (reinit effects' state after `reset` samples).
    
  If you prefer operating through [[filter-channels]] use [[effect-filter]]."
  {:metadoc/categories #{:eff}}
  ([effects config config-back reset pixels]
   (process-pixels pixels (p/clone-pixels pixels) effects config config-back reset))
  ([effects pixels]
   (apply-effects-to-pixels effects {} {} 0 pixels))
  ([effects reset pixels]
   (apply-effects-to-pixels effects {} {} reset pixels))
  ([effects config config-back pixels]
   (apply-effects-to-pixels effects config config-back 0 pixels)))

;; ## Helper functions

(defn db->linear
  "DB to Linear"
  ^double [^double x]
  (m/pow 10.0 (/ x 20.0)))

(defn linear->db
  "Linear to DB"
  ^double [x]
  (* 20.0 (m/log10 x)))

;; ## Effects / Filters

(defmulti effect
  "Create effect for given name (as keyword) and optional parametrization.

  List of all possible effects is under [[effects-list]].

  Effect is a custom type which contains: name, sample (result of last call), effect function and current state.

  Effect can be considered as function: call with sample to effect with next state or call without parameter to obtain latest result. Effects are also composable with [[compose-effects]]."
  {:metadoc/categories #{:eff}}
  (fn [m & conf] m))

;; ### Simple Low/High pass filters
(defn- calc-filter-alpha
  "Calculate alpha factor"
  ^double [^double rate ^double cutoff]
  (let [tinterval (/ rate)
        tau (/ (* cutoff m/TWO_PI))]
    (/ tinterval (+ tau tinterval))))

(defmethod effect :simple-lowpass
  ([m] (effect m {}))
  ([m {:keys [rate cutoff]
       :or {rate 44100.0 cutoff 2000.0}}]
   (let [alpha (calc-filter-alpha rate cutoff)] 
     (effect-node m (fn
                      ([^double sample ^double prev]
                       (let [s1 (* sample alpha)
                             s2 (- prev (* prev alpha))
                             nprev (+ s1 s2)]
                         (SampleAndState. nprev nprev)))
                      ([] 0.0))))))

(defmethod effect :simple-highpass
  ([m] (effect m {}))
  ([m conf]
   (let [lpfilter (effect :simple-lowpass conf)]
     (effect-node m (fn
                      ([^double sample lp]
                       (let [^EffectsList res (single-pass lp sample)]
                         (SampleAndState. (- sample (.sample res)) res)))
                      ([] lpfilter))))))

;; ### Biquad filters

;; Store biquad effect configuration in `BiquadConf` type
(deftype BiquadConf [^double b0 ^double b1 ^double b2 ^double a1 ^double a2])

(defn- biquad-eq-params
  "Calculate configuration for biquad equalizer
   fc - center frequency
   gain
   bw - bandwidth
   fs - sample rate"
  [^double fc ^double gain ^double bw ^double fs]
  (let [w (/ (* m/TWO_PI (m/constrain fc 1.0 (* 0.5 fs))) fs)
        cw (m/cos w)
        sw (m/sin w)
        J (m/pow 10.0 (* gain 0.025))
        g (-> bw
              (m/constrain 0.0001 4.0)
              (* m/LN2_2 w)
              (/ sw)
              (m/sinh)
              (* sw))
        a0r (/ (+ 1.0 (/ g J)))

        b0 (* a0r (+ 1.0 (* g J)))
        b1 (* a0r -2.0 cw)
        b2 (* a0r (- 1.0 (* g J)))
        a1 (- b1)
        a2 (* a0r (- (/ g J) 1.0))]
    (BiquadConf. b0 b1 b2 a1 a2)))

(defn- biquad-hs-params 
  "Calculate configuration for biquad high shelf
   fc - center frequency
   gain
   slope - shelf slope
   fs - sample rate"
  [^double fc ^double gain ^double slope ^double fs]
  (let [w (/ (* m/TWO_PI (m/constrain fc 1.0 (* 0.5 fs))) fs)
        cw (m/cos w)
        sw (m/sin w)
        A (m/pow 10.0 (* gain 0.025))
        iA (inc A)
        dA (dec A)
        b (m/sqrt (- (/ (inc (* A A)) (m/constrain slope 0.0001 1.0)) (* dA dA)))
        apc (* cw iA)
        amc (* cw dA)
        bs (* b sw)
        a0r (->> amc
                 (- iA)
                 (+ bs)
                 (/ 1.0))

        b0 (* a0r A (+ iA amc bs))
        b1 (* a0r A -2.0 (+ dA apc))
        b2 (* a0r A (- (+ iA amc) bs))
        a1 (* a0r -2.0 (- dA apc))
        a2 (* a0r (+ (dec (- A)) amc bs))]
    (BiquadConf. b0 b1 b2 a1 a2)))

(defn- biquad-ls-params 
  "Calculate configuration for biquad low shelf
   fc - center frequency
   gain
   slope - shelf slope
   fs - sample rate"
  [^double fc ^double gain ^double slope ^double fs]
  (let [w (/ (* m/TWO_PI (m/constrain fc 1.0 (* 0.5 fs))) fs)
        cw (m/cos w)
        sw (m/sin w)
        A (m/pow 10.0 (* gain 0.025))
        iA (inc A)
        dA (dec A)
        b (m/sqrt (- (/ (inc (* A A)) (m/constrain slope 0.0001 1.0)) (* dA dA)))
        apc (* cw iA)
        amc (* cw dA)
        bs (* b sw)
        a0r (->> amc
                 (+ iA)
                 (+ bs)
                 (/ 1.0))

        b0 (* a0r A (- (+ iA bs) amc))
        b1 (* a0r A 2.0 (- dA apc))
        b2 (* a0r A (- iA amc bs))
        a1 (* a0r 2.0 (+ dA apc))
        a2 (* a0r (+ bs (- (dec (- A)) amc)))]
    (BiquadConf. b0 b1 b2 a1 a2)))

(defn- biquad-lp-params
  "Calculate configuration for biquad low pass"
  [^double fc ^double bw ^double fs]
  (let [omega (* m/TWO_PI (/ fc fs))
        sn (m/sin omega)
        cs (m/cos omega)
        alpha (* sn (m/sinh (* m/LN2_2 bw (/ omega sn))))
        a0r (/ (inc alpha))
        cs- (- 1.0 cs)
        
        b0 (* a0r 0.5 cs-)
        b1 (* a0r cs-)
        b2 b0
        a1 (* a0r 2.0 cs)
        a2 (* a0r (dec alpha))]
    (BiquadConf. b0 b1 b2 a1 a2)))

(defn- biquad-hp-params
  "Calculate configuration for biquad high pass"
  [^double fc ^double bw ^double fs]
  (let [omega (* m/TWO_PI (/ fc fs))
        sn (m/sin omega)
        cs (m/cos omega)
        alpha (* sn (m/sinh (* m/LN2_2 bw (/ omega sn))))
        a0r (/ (inc alpha))
        cs+ (inc cs)
        
        b0 (* a0r 0.5 cs+)
        b1 (* a0r (- cs+))
        b2 b0
        a1 (* a0r 2.0 cs)
        a2 (* a0r (dec alpha))]
    (BiquadConf. b0 b1 b2 a1 a2)))

(defn- biquad-bp-params
  "Calculate configuration for biquad band pass"
  [^double fc ^double bw ^double fs]
  (let [omega (* m/TWO_PI (/ fc fs))
        sn (m/sin omega)
        cs (m/cos omega)
        alpha (* sn (m/sinh (* m/LN2_2 bw (/ omega sn))))
        a0r (/ (inc alpha))
        
        b0 (* a0r alpha)
        b1 0.0
        b2 (* a0r (- alpha))
        a1 (* a0r 2.0 cs)
        a2 (* a0r (dec alpha))]
    (BiquadConf. b0 b1 b2 a1 a2)))

;; Store state in `StateBiquad` type.
(deftype StateBiquad [^double x2 ^double x1 ^double y2 ^double y1])

(defn- make-biquad-filter
  "Create biquad effect based on passed configuration"
  [m ^BiquadConf c]
  (effect-node m (fn
                   ([^double sample ^StateBiquad state]
                    (let [y (-> (* (.b0 c) sample)
                                (+ (* (.b1 c) (.x1 state)))
                                (+ (* (.b2 c) (.x2 state)))
                                (+ (* (.a1 c) (.y1 state)))
                                (+ (* (.a2 c) (.y2 state))))]
                      (SampleAndState. y (StateBiquad. (.x1 state) sample (.y1 state) y))))
                   ([] (StateBiquad. 0.0 0.0 0.0 0.0)))))

;; ### Biquad equalizer

(defmethod effect :biquad-eq
  ([m] (effect m {}))
  ([m {:keys [fc gain bw fs]
       :or {fc 0.0 gain 0.0 bw 1.0 fs 44100.0}}]
   (make-biquad-filter m (biquad-eq-params fc gain bw fs))))

;; ### Biquad high/low shelf

(defmethod effect :biquad-hs
  ([m] (effect m {}))
  ([m {:keys [fc gain slope fs]
       :or {fc 0.0 gain 0.0 slope 1.5 fs 44100.0}}]
   (make-biquad-filter m (biquad-hs-params fc gain slope fs))))

(defmethod effect :biquad-ls
  ([m] (effect m {}))
  ([m {:keys [fc gain slope fs]
       :or {fc 0.0 gain 0.0 slope 1.5 fs 44100.0}}]
   (make-biquad-filter m (biquad-ls-params fc gain slope fs))))

;; ### Biquad lowpass/highpass/bandpass

(defn- lhb-params
  "Create parameters for lp hp and bp biquad filters."
  [f {:keys [fc bw fs]
      :or {fc 0.0 bw 1.0 fs 44100.0}}]
  (f fc bw fs))

(defmethod effect :biquad-lp ([m] (effect m {})) ([m conf] (make-biquad-filter m (lhb-params biquad-lp-params conf))))
(defmethod effect :biquad-hp ([m] (effect m {})) ([m conf] (make-biquad-filter m (lhb-params biquad-hp-params conf))))
(defmethod effect :biquad-bp ([m] (effect m {})) ([m conf] (make-biquad-filter m (lhb-params biquad-bp-params conf))))

;; ### DJ Equalizer

(defmethod effect :dj-eq
  ([m] (effect m {}))
  ([m {:keys [hi mid low shelf-slope peak-bw ^double rate]
       :or {hi 0.0 mid 0.0 low 0.0 shelf-slope 1.5 peak-bw 1.0 rate 44100.0}}]
   (let [b (compose-effects
            (effect :biquad-hs {:fc (* rate (/ 10000.0 44100.0)) :gain hi :slope shelf-slope :fs rate})
            (effect :biquad-eq {:fc (* rate (/ 1000.0 44100.0)) :gain mid :bw peak-bw :fs rate})
            (effect :biquad-eq {:fc (* rate (/ 100.0 44100.0)) :gain low :bw peak-bw :fs rate}))]
     (effect-node m (fn
                      ([sample state]
                       (let [^EffectsList res (single-pass state sample)]
                         (SampleAndState. (.sample res) res)))
                      ([] b))))))

;; ### Phaser

(defmethod effect :phaser-allpass
  ([m] (effect m {}))
  ([m {:keys [^double delay]
       :or {delay 0.5}}]
   (let [a1 (/ (- 1.0 delay) (inc delay))]
     (effect-node m (fn
                      ([^double sample ^double zm1]
                       (let [y (+ zm1 (* sample (- a1)))
                             new-zm1 (+ sample (* y a1))]
                         (SampleAndState. y new-zm1)))
                      ([] 0.0))))))

;; ### Divider

(deftype StateDivider [^double out ^double amp ^double count ^double lamp ^double last ^int zeroxs])

(defmethod effect :divider
  ([m] (effect m {}))
  ([m {:keys [^long denom]
       :or {denom 2.0}}]
   (effect-node m
                (fn
                  ([^double sample ^StateDivider state]
                   (let [count (inc (.count state))
                         ^StateDivider s1 (if (bool-or (bool-and (> sample 0.0) (<= (.last state) 0.0))
                                                       (bool-and (neg? sample) (>= (.last state) 0.0)))
                                            (if (== denom 1)
                                              (StateDivider. (if (pos? (.out state)) -1.0 1.0) 0.0 0.0 (/ (.amp state) count) (.last state) 0)
                                              (StateDivider. (.out state) (.amp state) count (.lamp state) (.last state) (inc (.zeroxs state))))
                                            (StateDivider. (.out state) (.amp state) count (.lamp state) (.last state) (.zeroxs state)))
                         amp (+ (.amp s1) (m/abs sample))
                         ^StateDivider s2 (if (bool-and (> denom 1)
                                                        (== ^long (rem (.zeroxs s1) denom) (dec denom)))
                                            (StateDivider. (if (pos? (.out s1)) -1.0 1.0) 0.0 0 (/ amp (.count s1)) (.last s1) 0)
                                            (StateDivider. (.out s1) amp (.count s1) (.lamp s1) (.last s1) (.zeroxs s1)))]
                     (SampleAndState. (* (.out s2) (.lamp s2)) (StateDivider. (.out s2) (.amp s2) (.count s2) (.lamp s2) sample (.zeroxs s2)))))
                  ([] (StateDivider. 1.0 0.0 0.0 0.0 0.0 0.0))))))

;; ### FM filter

(deftype StateFm [^double pre ^double integral ^double t lp])

(defmethod effect :fm
  ([m] (effect m {}))
  ([m {:keys [^double quant ^double omega ^double phase]
       :or {quant 10.0 omega 0.014 phase 0.00822}}]
   (let [lp-chain (compose-effects (effect :simple-lowpass {:rate 100000 :cutoff 25000})
                                   (effect :simple-lowpass {:rate 100000 :cutoff 10000})
                                   (effect :simple-lowpass {:rate 100000 :cutoff 1000}))]
     (effect-node m
                  (fn
                    ([^double sample ^StateFm state]
                     (let [sig (* sample phase)
                           new-integral (+ (.integral state) sig)
                           m (m/cos (+ new-integral (* omega (.t state))))
                           m (if (pos? quant)
                               (m/norm (unchecked-int (m/norm m -1.0 1.0 0.0 quant)) 0.0 quant -1.0 1.0)
                               m)
                           dem (m/abs (- m (.pre state)))
                           ^EffectsList res (single-pass (.lp state) dem)
                           demf (/ (* 2.0 (- (.sample res) omega)) phase)]
                       (SampleAndState. (m/constrain demf -1.0 1.0) (StateFm. m new-integral (inc (.t state)) res))))
                    ([] (StateFm. 0.0 0.0 0.0 lp-chain)))))))


;; ### Bandwidth limit

(defmethod effect :bandwidth-limit
  ([m] (effect m {}))
  ([m {:keys [^double freq ^double rate]
       :or {freq 1000.0 rate 44100.0}}]
   (let [dx (/ freq rate)]
     (effect-node m (fn
                      ([^double sample ^double state]
                       (let [res (if (>= sample state)
                                   (min (+ state dx) sample)
                                   (max (- state dx) sample))]
                         (SampleAndState. res res)))
                      ([] 0.0))))))

;; ### Distortion

(defmethod effect :distort
  ([m] (effect m {}))
  ([m {:keys [^double factor]
       :or {factor 1.0}}]
   (let [nfact (inc factor)]
     (effect-node m (fn 
                      ([^double sample state]
                       (let [div (+ factor (m/abs sample))
                             res (* nfact (/ sample div))]
                         (SampleAndState. res state)))
                      ([]))))))

;; ### Fast overdrive

(defmethod effect :foverdrive
  ([m] (effect m {}))
  ([m {:keys [^double drive]
       :or {drive 2.0}}]
   (let [drivem1 (dec drive)]
     (effect-node m (fn
                      ([^double sample state]
                       (let [fx (m/abs sample)
                             res (/ (* sample (+ fx drive)) (inc (+ (* sample sample) (* fx drivem1))))]
                         (SampleAndState. res state)))
                      ([]))))))

;; ### Decimator

(deftype StateDecimator [^double count ^double last])

(defmethod effect :decimator
  ([m] (effect m {}))
  ([m {:keys [^double bits ^double fs ^double rate]
       :or {bits 2.0 fs 4410.0 rate 44100.0}}]
   (let [step (m/pow 0.5 (- bits 0.9999))
         stepr (/ step)
         ratio (/ fs rate)]
     (effect-node m (fn
                      ([^double sample ^StateDecimator state]
                       (let [ncount (+ (.count state) ratio)]
                         (if (>= ncount 1.0)
                           (let [delta (* step ^double (m/remainder (->> sample
                                                                         m/sgn
                                                                         (* step 0.5)
                                                                         (+ sample)
                                                                         (* stepr)) 1.0))
                                 last (- sample delta)]
                             (SampleAndState. last (StateDecimator. (dec ncount) last)))
                           (SampleAndState. (.last state) (StateDecimator. ncount (.last state))))))
                      ([] (StateDecimator. 0.0 0.0)))))))

;; ### BassTreble

(deftype StateBassTreble [^double xn1Bass ^double xn2Bass ^double yn1Bass ^double yn2Bass
                          ^double xn1Treble ^double xn2Treble ^double yn1Treble ^double yn2Treble])

(defmethod effect :basstreble
  ([m] (effect m {}))
  ([m {:keys [^double bass ^double treble ^double gain ^double rate ^double slope ^double bass-freq ^double treble-freq]
       :or {bass 1.0 treble 1.0 gain 0.0 rate 44100.0 slope 0.4 bass-freq 250.0 treble-freq 4000.0}}]
   (let [data-gain (db->linear gain)
         wb (/ (* m/TWO_PI bass-freq) rate)
         wt (/ (* m/TWO_PI treble-freq) rate) 
         cwb (m/cos wb)
         cwt (m/cos wt)
         ab (m/exp (/ (* 2.302585092994046 bass) 40.0))
         ab+ (inc ab)
         ab- (dec ab)
         at (m/exp (/ (* 2.302585092994046 treble) 40.0))
         at+ (inc at)
         at- (dec at)
         bb (m/sqrt (- (/ (inc (m/sq ab)) slope) (m/sq (dec ab))))
         bt (m/sqrt (- (/ (inc (m/sq at)) slope) (m/sq (dec at))))
         bswb (* bb (m/sin wb))
         bswt (* bt (m/sin wt))

         b0b (* ab (+ (- ab+ (* ab- cwb)) bswb))
         b1b (* 2.0 ab (- ab- (* ab+ cwb)))
         b2b (* ab (- (- ab+ (* ab- cwb)) bswb))
         a0b (+ (+ ab+ (* ab- cwb)) bswb)
         a1b (* -2.0 (+ ab- (* ab+ cwb)))
         a2b (- (+ ab+ (* ab- cwb)) bswb)

         b0t (* at (+ (+ at+ (* at- cwt)) bswt))
         b1t (* -2.0 at (+ at- (* at+ cwt)))
         b2t (* at (- (+ at+ (* at- cwt)) bswt))
         a0t (+ (- at+ (* at- cwt)) bswt)
         a1t (* 2.0 (- at- (* at+ cwt)))
         a2t (- (- at+ (* at- cwt)) bswt)]
     (effect-node m (fn
                      ([^double sample ^StateBassTreble state]
                       (let [outb (/ (-> (* b0b sample)
                                         (+ (* b1b (.xn1Bass state)))
                                         (+ (* b2b (.xn2Bass state)))
                                         (- (* a1b (.yn1Bass state)))
                                         (- (* a2b (.yn2Bass state)))) a0b)
                             outt (/ (-> (* b0t outb)
                                         (+ (* b1t (.xn1Treble state)))
                                         (+ (* b2t (.xn2Treble state)))
                                         (- (* a1t (.yn1Treble state)))
                                         (- (* a2t (.yn2Treble state)))) a0t)]
                         (SampleAndState. (* outt data-gain)
                                          (StateBassTreble. sample (.xn1Bass state) outb (.yn1Bass state)
                                                            outb (.xn1Treble state) outt (.yn1Treble state)))))
                      ([] (StateBassTreble. 0.0 0.0 0.0 0.0
                                            0.0 0.0 0.0 0.0)))))))

;; ### Echo (audacity)

(deftype StateEcho [^doubles buffer ^int position])

(defmethod effect :echo
  ([m] (effect m {}))
  ([m {:keys [^double delay ^double decay ^double rate]
       :or {delay 0.5 decay 0.5 rate 44100.0}}]
   (let [buffer-len (int (min 10000000 (* delay rate)))]
     (effect-node m (fn
                      ([^double sample ^StateEcho state]
                       (let [result (+ sample (* decay (aget ^doubles (.buffer state) (.position state))))]
                         (aset ^doubles (.buffer state) (.position state) result)
                         (SampleAndState. result (StateEcho. (.buffer state) (rem (inc (.position state)) buffer-len)))))
                      ([]
                       (StateEcho. (double-array buffer-len 0.0) 0)))))))

;; ### Vcf303
;;
(deftype StateVcf303 [^double d1 ^double d2 ^double c0 ^int env-pos ^Vec3 abc])

(defmethod effect :vcf303
  ([m] (effect m {}))
  ([m {:keys [^double rate trigger ^double cutoff ^double resonance ^double env-mod ^double decay ^double gain]
       :or {rate 44100.0 trigger false cutoff 0.8 resonance 0.8 env-mod 0.5 decay 1.0 gain 1.0}}]
   (let [scale (/ m/PI rate)
         e0 (* scale
               (m/exp (-> (- 5.613 (* 0.8 env-mod))
                          (+ (* 2.1553 cutoff))
                          (- (* 0.7696 (- 1.0 resonance))))))        
         d (m/pow (->> decay
                       (* 2.3)
                       (+ 0.2)
                       (* rate)
                       (/ 1.0)
                       (m/pow 0.1)) 64.0)
         r (m/exp (- (* 3.455 resonance) 1.20))
         recalc-abc (fn [^double vc0]
                      (let [whopping (+ e0 vc0)
                            k (m/exp (/ (- whopping) r))
                            a (* (+ k k) (m/cos (+ whopping whopping)))
                            b (* (- k) k)
                            c (* 0.2 (- (- 1.0 a) b))]
                        (Vec3. a b c)))
         init-c0       (if trigger
                         (- (* scale
                               (m/exp (-> (+ 6.109 (* 1.5876 env-mod))
                                          (+ (* 2.1553 cutoff))
                                          (- (* 1.2 (- 1.0 resonance)))))) e0)
                         0.0)]    
     (effect-node m (fn
                      ([^double sample ^StateVcf303 state]
                       (let [^Vec3 abc (.abc state)
                             result (-> (* (.x abc) (.d1 state))
                                        (+ (* (.y abc) (.d2 state)))
                                        (+ (* (.z abc) sample)))
                             d2 (.d1 state)
                             d1 result
                             env-pos (inc (.env-pos state))]
                         (if (>= env-pos 64)
                           (let [c0 (* d (.c0 state))]
                             (SampleAndState. (* gain result)
                                              (StateVcf303. d1 d2 c0 0 (recalc-abc c0))))
                           (SampleAndState. (* gain result)
                                            (StateVcf303. d1 d2 (.c0 state) env-pos abc)))))
                      ([] (StateVcf303. 0.0 0.0 init-c0 0 (recalc-abc init-c0))))))))

;; ### Slew limiter

(defmethod effect :slew-limit
  ([m] (effect m {}))
  ([m {:keys [^double rate ^double maxrise ^double maxfall]
       :or {rate 44100.0 maxrise 500.0 maxfall 500.0}}]
   (let [maxinc (/ maxrise rate)
         maxdec (- (/ maxfall rate))]
     (effect-node m (fn
                      ([^double sample ^double prev]
                       (let [increment (- sample prev) 
                             nsample (+ prev (m/constrain increment maxdec maxinc))]
                         (SampleAndState. nsample nsample)))
                      ([] 0.0))))))


;; 
(deftype StateMdaThruZero [^doubles buffer ^double ph ^long bp ^double f])

(defmethod effect :mda-thru-zero
  ([m] (effect m {}))
  ([m {:keys [^double rate ^double speed ^double depth ^double mix ^double depth-mod ^double feedback]
       :or {rate 44100.0 speed 0.3 depth 0.43 mix 0.47 feedback 0.3 depth-mod 1.0}}]
   (let [rat (/ (* (m/pow 10.0 (- 2.0 (* 3.0 speed))) 2.0) rate)
         dep (* 2000.0 (m/sq depth))
         dem (- dep (* dep depth-mod))
         dep (- dep dem)
         wet mix
         dry (- 1.0 wet)
         fb (- (* 1.9 feedback) 0.95)]
     (effect-node m (fn
                      ([^double sample ^StateMdaThruZero state]
                       (let [ph (+ (.ph state) rat)
                             ph (if (> ph 1.0) (- ph 2.0) ph)
                             bp (bit-and (dec (.bp state)) 0x7ff)]
                         (aset ^doubles (.buffer state) bp (+ sample (* fb (.f state))))
                         (let [tmpf (+ dem (* dep (- 1.0 (m/sq ph))))
                               tmp (unchecked-int tmpf)
                               tmpf (- tmpf tmp)
                               tmp (bit-and (+ tmp bp) 0x7ff)
                               tmpi (bit-and (inc tmp) 0x7ff)
                               f (aget ^doubles (.buffer state) tmp)
                               f (+ (* tmpf (- (aget ^doubles (.buffer state) tmpi) f)) f)
                               result (+ (* sample dry) (* f wet))]
                           (SampleAndState. result (StateMdaThruZero. (.buffer state) ph bp f)))))
                      ([] (StateMdaThruZero. (double-array 2048) 0.0 0 0.0)))))))

(def ^{:metadoc/categories #{:eff}
       :doc "List of effects."}
  effects-list (sort (keys (methods effect))))

;; ## File operations

(defn save-signal
  "Save signal to file.

  Representation is: 16 bit signed, big endian file
  You can use Audacity/SOX utilities to convert files to audio."
  {:metadoc/categories #{:sig}}
  [^Signal s filename]
  (make-parents filename)
  (let [^java.io.DataOutputStream out (java.io.DataOutputStream. (output-stream filename))]
    (try
      (dotimes [i (alength ^doubles (.signal s))]
        (.writeShort out (short (m/cnorm (aget ^doubles (.signal s) i) -1.0 1.0 Short/MIN_VALUE Short/MAX_VALUE))))
      (.flush out)
      (finally (. out clojure.core/close)))
    s))

(defn load-signal
  "Read signal from file

  Expected representation is 16 bit signed, big endian file."
  {:metadoc/categories #{:sig}}
  [filename]
  (let [^java.io.File f (file filename)
        len (/ (.length f) 2)
        ^java.io.DataInputStream in (java.io.DataInputStream. (input-stream filename))
        ^doubles buffer (double-array len)]
    (try
      (dotimes [i len]
        (aset ^doubles buffer (int i) (double (m/cnorm (.readShort in) Short/MIN_VALUE Short/MAX_VALUE -1.0 1.0))))
      (finally (. in clojure.core/close)))
    (Signal. buffer)))

;; ## Signal generators
;;
;; Here you have defined multimethods to create waves from various oscilators
;;
;; Parameters are:
;;
;; * oscilator name (see `oscillators` variable)
;; * frequency
;; * amplitude
;; * phase (0-1)
;;
;; Multimethod creates oscillator function accepting `double` (time) and resulting `double` from [-1.0 1.0] range.

(defmulti wave
  "Create waves from various oscilators

  Parameters are:

  * oscilator name (see `oscillators` variable)
  * frequency
  * amplitude
  * phase (0-1)
  
  Multimethod creates oscillator function accepting `double` (time) and resulting `double` from [-1.0 1.0] range.

  To convert `wave` to Signal, call [[signal-from-wave]].

  To add waves, call [[sum-waves]]."
  {:metadoc/categories #{:wave}}
  (fn [f _ _ _] f))

(defmethod wave :sin [_ ^double f ^double a ^double p]
  (fn ^double [^double x]
    (* a
       (m/sin (+ (* p m/TWO_PI) (* x m/TWO_PI f))))))

(def ^:private snoise (r/fbm-noise {:octaves 2}))

(defmethod wave :noise [_ ^double f ^double a ^double p]
  (fn ^double [^double x]
    (* a 2.0
       (- ^double (snoise (* (+ p x) f) 1.23456789) 0.5))))

(defmethod wave :saw [_ ^double f ^double a ^double p] 
  (fn ^double [^double x]
    (let [rp (* 2.0 a)
          p2 (* f ^double (mod (+ (* a p) a x) 1.0))]
      (* rp (- p2 (m/floor p2) 0.5)))))

(defmethod wave :square [_ ^double f ^double a ^double p]
  (fn ^double [^double x]
    (if (< ^double (mod (+ p (* x f)) 1.0) 0.5)
      a
      (- a))))

(defmethod wave :triangle [_ ^double f ^double a ^double p]
  (let [saw (wave :saw f a p)]
    (fn ^double [^double x]
      (- (* 2.0 (m/abs (saw x))) a))))

(defmethod wave :cut-triangle [_ ^double f ^double a ^double p]
  (let [tri (wave :triangle f a p)]
    (fn ^double [^double x]
      (let [namp (* 0.5 a)]
        (* 2.0 (m/constrain ^double (tri x) (- namp) namp))))))

(defmethod wave :constant [_ _ ^double a _] (constantly a))

(def ^{:doc "List of oscillator names used with [[wave]]"
       :metadoc/categories #{:wave}}
  oscillators (sort (keys (methods wave))))

(defn sum-waves
  "Create wave which is sum of all waves."
  {:metadoc/categories #{:wave}}
  [& fs]
  (reduce #(fn [^double x] (+ ^double (%1 x) ^double (%2 x))) fs))

(defn wave->signal
  "Create Signal from wave. Parameters are: f - oscillator, samplerate and number of seconds."
  {:metadoc/categories #{:wave}}
  [f ^double samplerate ^double seconds]
  (let [len (* samplerate seconds)
        ^doubles buffer (double-array len)]
    (dotimes [i len]
      (aset ^doubles buffer i (m/constrain ^double (f (m/norm i 0 len 0 seconds)) -1.0 1.0)))
    (Signal. buffer)))