;;; -*- syntax: common-lisp; package: clm; base: 10; mode: lisp -*-

;;; this is a translation to CLM of Perry Cook's Physical Modelling Toolkit.

;;; reedtable
(def-clm-struct reed (offset 0.6) (slope -0.8))

(defmacro reedtable (r sample) 
  `(min 1.0 (+ (reed-offset ,r) (* (reed-slope ,r) ,sample))))

;;; bowtable
(def-clm-struct bowt (offset 0.0) (slope 1.0))

(defmacro bowtable (b sample) 
  `(max 0.0 (- 1.0 (abs (* (bowt-slope ,b) (+ ,sample (bowt-offset ,b)))))))

;;;jettable
(defmacro jettable (sample0) 
  `(let ((sample ,sample0)) (max -1.0 (min 1.0 (* sample (- (* sample sample) 1.0))))))

;;; one-zero filter (slightly different from CLM's)
(def-clm-struct onez (gain 0.5) (zerocoeff 1.0) (input 0.0))

(defmacro onezero (b sample0)
  `(let ((sample ,sample0))
     (prog1
	 (* (onez-gain ,b) (+ sample (* (onez-zerocoeff ,b) (onez-input ,b))))
       (setf (onez-input ,b) sample))))

;;; one-pole filter (also slightly different)
(def-clm-struct onep (polecoeff 0.9) (gain 1.0) (sgain 0.1) (output 0.0))

(defun fixup-sgain (p)
  (if (> (onep-polecoeff p) 0.0) 
      (setf (onep-sgain p) (* (onep-gain p) (- 1.0 (onep-polecoeff p))))
    (setf (onep-sgain p) (* (onep-gain p) (+ 1.0 (onep-polecoeff p))))))

(defun set-pole (p aval) (setf (onep-polecoeff p) aval) (fixup-sgain p))

(defun set-gain (p aval) (setf (onep-gain p) aval) (fixup-sgain p))

(defmacro onepole (p sample0) 
  `(let ((sample ,sample0))
     (setf (onep-output ,p) (+ (* sample (onep-sgain ,p)) 
			       (* (onep-polecoeff ,p) (onep-output ,p))))))

;;; biquadfilter (a version of CLM's formnt generator)
(def-clm-struct biq (pc0 0.0) (pc1 0.0) (zc0 0.0) (zc1 0.0) (gain 1.0) (out0 0.0) (out1 0.0) (in0 0.0) (in1 0.0))

(defmacro biquad (b sample0)
  `(let ((sample ,sample0))
     (let ((temp (+ (* (biq-zc0 ,b) (biq-in0 ,b)) (* (biq-zc1 ,b) (biq-in1 ,b)))))
       (setf (biq-in1 ,b) (biq-in0 ,b))
       (setf (biq-in0 ,b) (* (biq-gain ,b) sample))
       (incf temp (+ (biq-in0 ,b) (* (biq-pc0 ,b) (biq-out0 ,b)) (* (biq-pc1 ,b) (biq-out1 ,b))))
       (setf (biq-out1 ,b) (biq-out0 ,b))
       (setf (biq-out0 ,b) temp))))

;;;lipfilter (a biquad filter)
(defun lip-set-freq (b freq) 
  (setf (biq-pc0 b) (* 2.0 0.999 (cos (/ (* pi 2 freq) sampling-rate))))
  (setf (biq-pc1 b) (* -0.999 0.999))
  (setf (biq-gain b) 0.02))

(defmacro lip (b mouthsample0 boresample0)
  `(let ((mouthsample ,mouthsample0)
	 (boresample ,boresample0))
     (let ((temp (biquad ,b (- mouthsample boresample))))
       (setf temp (min 1.0 (* temp temp)))
       (+ (* temp mouthsample) (* (- 1.0 temp) boresample)))))

(def-clm-struct dcb (input 0.0) (output 0.0))
(defmacro dcblock (b sample0) 
  `(let ((sample ,sample0))
     (prog1
	 (setf (dcb-output ,b) (+ sample (- (* 0.99 (dcb-output ,b)) (dcb-input ,b))))
       (setf (dcb-input ,b) sample))))
;;; in the sndtools program sndblockdc, the 0.99 is replaced by (- 1.0 (/ 7.0 adaption_time))

;;;delaylineA
(def-clm-struct dla (inpoint INTEGER) (outpoint INTEGER) (lastin 0.0) (length INTEGER) (output 0.0) (input ARRAY) alpha coeff)

(defun make-delayA (len) 
  (let ((nd (make-dla :length len :input (make-short-float-array len :initial-element 0.0) :inpoint 0 :outpoint 0)))
    (set-delayA nd (* 0.5 len))
    nd))

(defun set-delayA (d lag)
  (let ((outpointer (+ (dla-inpoint d) (- 2.0 lag))))
    (loop while (minusp outpointer) do (incf outpointer (dla-length d)))
    (setf (dla-outpoint d) (floor outpointer))
    (setf (dla-alpha d) (- outpointer (dla-outpoint d)))
    (setf (dla-coeff d) (/ (- 1.0 (dla-alpha d)) (+ 1.0 (dla-alpha d))))))

(defmacro delayA (d sample0)
  `(let ((sample ,sample0))
     (let ((temp 0.0))
       (setf (aref (dla-input ,d) (dla-inpoint ,d)) sample)
       (incf (dla-inpoint ,d))
       (if (= (dla-inpoint ,d) (dla-length ,d)) (setf (dla-inpoint ,d) 0))
       (setf temp (aref (dla-input ,d) (dla-outpoint ,d)))
       (incf (dla-outpoint ,d))
       (if (>= (dla-outpoint ,d) (dla-length ,d)) (decf (dla-outpoint ,d) (dla-length ,d)))
       (setf (dla-output ,d) (+ (* (- (dla-coeff ,d)) (dla-output ,d)) (dla-lastin ,d) (* temp (dla-coeff ,d))))
       (setf (dla-lastin ,d) temp)
       (dla-output ,d))))

;;; DelayLineL
(def-clm-struct dll (inpoint INTEGER) (outpoint INTEGER) (length INTEGER) (output 0.0) (input ARRAY) alpha omAlpha)

(defun make-delayL (len) 
  (let ((nd (make-dll :length len :input (make-short-float-array len :initial-element 0.0) :inpoint 0 :outpoint 0)))
    (set-delayL nd (* 0.5 len))
    nd))

(defun set-delayL (d lag)
  (let ((outpointer (+ (dll-inpoint d) (- 1 lag))))
    (loop while (minusp outpointer) do (incf outpointer (dll-length d)))
    (setf (dll-outpoint d) (floor outpointer))
    (setf (dll-alpha d) (- outpointer (dll-outpoint d)))
    (setf (dll-omalpha d) (- 1.0 (dll-alpha d)))))

(defmacro delayL (d sample0)
  `(let ((sample ,sample0))
     (setf (aref (dll-input ,d) (dll-inpoint ,d)) sample)
     (incf (dll-inpoint ,d))
     (if (= (dll-inpoint ,d) (dll-length ,d)) (setf (dll-inpoint ,d) 0))
     (setf (dll-output ,d) (* (aref (dll-input ,d) (dll-outpoint ,d)) (dll-omalpha ,d)))
     (incf (dll-outpoint ,d))
     (if (= (dll-outpoint ,d) (dll-length ,d)) (setf (dll-outpoint ,d) 0))
     (incf (dll-output ,d) (* (aref (dll-input ,d) (dll-outpoint ,d)) (dll-alpha ,d)))))


;;; now some example instruments

(defcinstrument plucky (beg dur freq amplitude &optional (maxA 1.0))
  (let* ((lowestfreq 100.0)
	 (len (1+ (floor (/ sampling-rate lowestfreq))))
	 (delayLine (make-delayA len))
	 (filter (make-onez))
	 (st (floor (* sampling-rate beg)))
	 (nd (+ st (floor (* sampling-rate dur)))))
    (set-delayA delayLine (- (/ sampling-rate freq) 0.5))
    (loop for i from 0 below len do
      (delayA delayLine (make-short-float (+ (* 0.99 (dla-output delayLine)) (* maxA (- 1.0 (random 2.0)))))))
    (run
     (loop for i from st to nd do
       (outa i (* amplitude (delayA delayLine (onezero filter (dla-output delayLine)))))))))

;;; freq is off in this one (in PRC's original also)
(defcinstrument bow (beg dur frq amplitude &optional (maxA 1.0))
  (let* ((lowestFreq 100.0)
	 (len (1+ (floor (/ sampling-rate lowestFreq))))
	 (neckdelay (make-delayL len))
	 (bridgedelay (make-delayL (floor len 2)))
	 (bowtab (make-bowt :slope 3.0))
	 (filt (make-onep))
	 (rate .001)
	 (bowing t)
	 (bowvelocity rate)
	 (maxvelocity maxA)
	 (attackrate rate)
	 (st (floor (* sampling-rate beg)))
	 (nd (+ st (floor (* sampling-rate dur))))
	 (release (+ st (floor (* .8 nd)))))
    (set-pole filt 0.6)
    (set-gain filt 0.3)
    (let ((ratio 0.8317)
	  (temp (- (/ sampling-rate frq) 4.0)))
      (set-delayL neckdelay (* temp ratio))
      (set-delayL bridgedelay (* temp (- 1.0 ratio))))
    (run
     (loop for i from st to nd do 
       (let* ((bridgerefl 0.0)
	      (nutrefl 0.0) 
	      (veldiff 0.0) 
	      (stringvel 0.0) 
	      (bowtemp 0.0))
	 (if bowing
	     (if (/= maxvelocity bowvelocity)
		 (if (< bowvelocity maxvelocity)
		     (incf bowvelocity attackrate)
		   (decf bowvelocity attackrate)))
	   (if (> bowvelocity 0.0) (decf bowvelocity attackrate)))
	 (setf bowtemp (* 0.3 bowvelocity))
	 (setf bridgerefl (- (onepole filt (dll-output bridgedelay))))
	 (setf nutrefl (- (dll-output neckdelay)))
	 (setf stringvel (+ bridgerefl nutrefl))
	 (setf veldiff (- bowtemp stringvel))
	 (setf veldiff (* veldiff (bowtable bowtab veldiff)))
	 (delayL neckdelay (+ bridgerefl veldiff))
	 (delayL bridgedelay (+ nutrefl veldiff))
	 (outa i (* amplitude 10.0 (onep-output filt)))
	 (when (= i release) 
	   (setf bowing nil)
	   (setf attackrate .0005)))))))

(defcinstrument brass (beg dur freq amplitude &optional (maxA 1.0))
  (let* ((lowestfreq 100.0)
	 (len (1+ (floor (/ sampling-rate lowestfreq))))
	 (delayLine (make-delayA len))
	 (lipFilter (make-biq))
	 (dcBlocker (make-dcb))
	 (blowing t)
	 (rate .001)
	 (breathPressure 0.0) 
	 (maxPressure maxA)
	 (attackRate rate)
	 (st (floor (* sampling-rate beg)))
	 (nd (+ st (floor (* sampling-rate dur))))
	 (release (+ st (floor (* .8 nd)))))
    (set-delayA delayLine (+ 1.0 (/ sampling-rate freq)))
    (lip-set-freq lipFilter freq)
    (run
     (loop for i from st to nd do
       (if blowing
	   (if (/= maxPressure breathPressure)
	       (if (< breathPressure maxPressure)
		   (incf breathPressure attackrate)
		 (decf breathPressure attackrate)))
	 (if (> breathPressure 0.0)
	     (decf breathPressure attackrate)))
       (outa i (* amplitude (delayA delayLine
				    (dcblock dcblocker
					     (lip lipFilter
						  (* 0.3 breathPressure)
						  (* 0.9 (dla-output delayLine)))))))
       (when (= i release) 
	 (setf blowing nil)
	 (setf attackrate .0005))))))

(defcinstrument clarinet (beg dur freq amplitude &optional (maxA 1.0))
  (let* ((lowestfreq 100.0)
	 (len (1+ (floor (/ sampling-rate lowestfreq))))
	 (delayLine (make-delayL len))
	 (reedTable (make-reed :offset 0.7 :slope -0.3))
	 (filter (make-onez))
	 (blowing t)
	 (breathPressure 0.0)
	 (rate .001)
	 (maxPressure maxA)
	 (attackRate rate)
	 (st (floor (* sampling-rate beg)))
	 (nd (+ st (floor (* sampling-rate dur))))
	 (release (+ st (floor (* .8 nd)))))
    (set-delayL delayLine (- (* 0.5 (/ sampling-rate freq)) 1.0))
    (run
     (loop for i from st to nd do
       (let ((pressurediff 0.0))
	 (if blowing
	     (if (/= maxPressure breathPressure)
		 (if (< breathPressure maxPressure)
		     (incf breathPressure attackrate)
		   (decf breathPressure attackrate)))
	   (if (> breathPressure 0.0)
	       (decf breathPressure attackrate)))
	 (setf pressurediff (- (onezero filter (* -0.95 (dll-output delayLine))) breathpressure))
	 (outa i (* amplitude (delayL delayLine (+ breathpressure (* pressurediff (reedtable reedtable pressurediff))))))
	 (when (= i release)
	   (setf blowing nil)
	   (setf attackrate .0005)))))))

(defcinstrument flute (beg dur freq amplitude &optional (maxA 1.0))
  (let* ((lowestfreq 100.0)
	 (len (1+ (floor (/ sampling-rate lowestfreq))))
	 (jetDelay (make-delayL (floor len 2)))
	 (boreDelay (make-delayL len))
	 (filter (make-onep))
	 (dcBlocker (make-dcb))
	 (jetRefl 0.6)
	 (endRefl 0.6)
	 (sinPhase 0.0)
	 (blowing t)
	 (rate .0005)
	 (breathPressure 0.0) 
	 (maxPressure maxA)
	 (attackRate rate)
	 (st (floor (* sampling-rate beg)))
	 (nd (+ st (floor (* sampling-rate dur))))
	 (release (+ st (floor (* .8 nd)))))
    (set-pole filter 0.8)
    (set-gain filter -1.0)
    (let ((ratio 0.8)
	  (temp (- (/ sampling-rate freq) 5.0)))
      (set-delayL boreDelay (* ratio temp))
      (set-delayL jetDelay (* temp (- 1.0 ratio))))
    (run
     (loop for i from st to nd do
       (let ((randPressure (* 0.1 breathPressure (random 1.0)))
	     (temp 0.0) 
	     (pressurediff 0.0))
	 (incf sinPhase 0.0007)		;5 Hz vibrato?
	 (if (> sinPhase 6.28) (decf sinPhase 6.28))
	 (incf randPressure (* 0.05 breathPressure (sin sinPhase)))
	 (if blowing
	     (if (/= maxPressure breathPressure)
		 (if (< breathPressure maxPressure)
		     (incf breathPressure attackrate)
		   (decf breathPressure attackrate)))
	   (if (> breathPressure 0.0) (decf breathPressure attackrate)))
	 (setf temp (dcblock dcBlocker (onepole filter (dll-output boreDelay))))
	 (setf pressurediff (+ (jettable (delayL jetDelay (+ breathpressure (- randpressure (* jetrefl temp))))) (* endrefl temp)))
	 (outa i (* 0.3 amplitude (delayL boreDelay pressurediff)))
	 (when (= i release)
	   (setf blowing nil)
	   (setf attackrate .0005)))))))


#|
      (with-sound () 
	(plucky 0 .3 440 .2)
	(bow .5 .3 220 .2)
	(brass 1 .3 440 .2)
	(clarinet 1.5 .3 440 .2)
	(flute 2 .3 440 .2))
|#