;;; -*- Lisp -*-

(defun pqw-normalize-to-power-of-two (coeffs)
  (let ((maxval (abs (aref coeffs 0)))
	(arrsiz (array-total-size coeffs))
	(scaler 0.0))
    (loop for i from 1 below arrsiz do
      (setf maxval (max maxval (abs (aref coeffs i)))))
    (setf maxval (expt 2 (+ 3 (ceiling (/ (log maxval) (log 2))))))
    ;; 3 here to give us plenty of headroom (to accomodate fractional arithmetic on 56000)
    (setf scaler (/ 1.0 maxval))
    (loop for i from 0 below arrsiz do
      (setf (aref coeffs i) (* scaler (aref coeffs i))))
    maxval))

(defvar pqw-degree 0.0)
(defvar pqw-distance 1.0)
(defvar pqw-reverb-amount 0.005)

(definstrument Pqw (start duration spacing-freq carrier-freq amplitude ampfun indexfun partials
		    &key (degree pqw-degree)
			 (distance pqw-distance)
			 (reverb-amount pqw-reverb-amount))
  ;; phase-quadrature waveshaping used to create asymmetric (i.e. single side-band) spectra.
  ;; The basic idea here is a variant of sin x sin y - cos x cos y = cos (x + y)
  (let* ((normalized-partials (normalize-partials partials))
	 (spacing-cos (make-oscil :frequency spacing-freq :initial-phase half-pi))
	 (spacing-sin (make-oscil :frequency spacing-freq))
	 (carrier-cos (make-oscil :frequency carrier-freq :initial-phase half-pi))
	 (carrier-sin (make-oscil :frequency carrier-freq))
	 (sin-coeffs (get-chebychev-coefficients normalized-partials 0))
	 (cos-coeffs (get-chebychev-coefficients normalized-partials 1))
	 (cos-norm (pqw-normalize-to-power-of-two cos-coeffs))
	 (sin-norm (pqw-normalize-to-power-of-two sin-coeffs))
	 (amp-env (make-func :envelope ampfun :scaler amplitude :start-time start :duration duration))
	 (ind-env (make-func :envelope indexfun :scaler 1.0 :start-time start :duration duration))
	 (loc (make-locsig :degree degree :distance distance :revscale reverb-amount))
	 (vib 0.0)
	 (ax 0.0)
	 (fax 0.0)
	 (yfax 0.0)
	 (r (/ carrier-freq spacing-freq))
	 (tr (make-triangle-wave :frequency 5 :amplitude (in-Hz (* .005 spacing-freq))))
	 (rn (make-randi :frequency 12 :amplitude (in-Hz (* .005 spacing-freq))))
	 (beg (floor (* start sampling-rate)))
	 (end (+ beg (floor (* duration sampling-rate)))))
    (Run
     (loop for i from beg to end do
       (setf vib (+ (triangle-wave tr) (randi rn)))
       (setf ax (* (min 1.0 (env ind-env)) (oscil spacing-cos vib)))
       (setf fax (* cos-norm (polynomial cos-coeffs ax)))
       (setf yfax (* (oscil spacing-sin vib) sin-norm (polynomial sin-coeffs ax)))
       (locsig loc i (* (env amp-env)
			(- (* (oscil carrier-sin (* vib r)) yfax) 
			   (* (oscil carrier-cos (* vib r)) fax))))))))

; (pqw 0 .5 200 1000 .2 '(0 0 25 1 100 0) '(0 1 100 0) '(2 .1 3 .3 6 .5))
; to see the asymmetric spectrum most clearly, set the index function above to '(0 1 100 1)
