;;; -*- syntax: common-lisp; package: clm; base: 10; mode:lisp -*-
;;;
;;; HRTF localization instrument
;;; Fernando Lopez-Lezcano
;;;
;;; April 13 1998 ... April 21 1998
;;; Feb 2001: translated to clm-2 [nando]
;;; Feb 2004: translated to clm-3 [nando]

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; HRTF data loading (for the MIT kemar hrtf's)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defparameter hrtf-fir-file "/usr/ccrma/snd/nando/hrtf/fir-coeffs.lisp")
(defparameter hrtf-fir-data '(()()))
(defparameter hrtf-dir "/usr/ccrma/snd/nando/hrtf/diffuse/*")
(defparameter hrtf-map "/usr/ccrma/snd/nando/hrtf/diffuse-map.lisp")
(defparameter hrtf-size 128)
(defparameter hrtf-data nil)
(defparameter hrtf-raw nil)

(defparameter hrtf-impulses (make-array (* 2 400 hrtf-size)
					:element-type 'float
					:adjustable t))

(defun load-fir-data (&key
		      (file hrtf-fir-file)
		      (verbose))
  (if (not verbose)
      (format t "; loading fir filter coefficients...~%"))
  (load file))

(defun load-hrtf-data (&key 
		       (dir hrtf-dir)
		       (size hrtf-size)
		       (verbose))
  (if (not verbose)
      (format t "; loading hrtf data set...~%"))
  ;; load the components of the hrtf in list form
  (load hrtf-map)
  (loop with offset = 0
	with current
	with data
	for oneelevation in hrtf-raw do
	(loop for description in oneelevation 
	      for name = (first description)
	      for elevation = (second description)
	      for azimut = (third description)
	      do
	      (if verbose (format t "reading ~a --> [~s:~s]...~%" name elevation azimut))
	      (with-open-file (in name
				  :direction :INPUT
				  :element-type '(unsigned-byte 8))
		;; read in the impulse response for the given elevation and azimut
		(loop for i from 0 below size do
		      (let* ((ia (+ (* 256 (read-byte in))(read-byte in)))
			     (ib (+ (* 256 (read-byte in))(read-byte in))))
			(setf (aref hrtf-impulses (+ offset i))
			      (/ (float (if (< ia 32768) ia (- ia 65536))) 32768.0)
			      (aref hrtf-impulses (+ offset size i))
			      (/ (float (if (< ib 32768) ib (- ib 65536))) 32768.0))))
		;; merge each azimut list into the current elevation
		(setf current (if (not current)
				  (list (list azimut offset (+ offset size)))
				  (merge 'list 
					 current 
					 (list (list azimut offset (+ offset size)))
					 #'< :key #'car)))
		(incf offset (* 2 size)))
	      finally
	      ;; merge each elevation into the dataset
	      (setf data (if (not data)
			     (list (list elevation current))
			     (merge 'list
				    data
				    (list (list elevation current))
				    #'< :key #'car))))
	finally
	;; set the global data, maybe we should just return the list
	(if verbose (format t "~s hrtf sets loaded~%" (/ offset 2 hrtf-size)))
	(setf hrtf-data data)))


(definstrument hrtf (start-time amplitude file
		     &key 
		     (elevation 0)
		     (azimut 0)
		     (reverb-amount 0.1)
		     (headphones nil))
  (let* (e-min e-max e-pos e-pos-gt a-pos a-pos-gt e a impulse-l impulse-r)
    ;; load all impulse responses if necessary
    (if (not hrtf-data)
	(load-hrtf-data))
    ;; load the crosstalk cancellation filter coefficients if necessary
    (if (not (first hrtf-fir-data))
	(load-fir-data))
    ;; sanity checks for elevation data
    (setf e-min (first (first hrtf-data))
	  e-max (first (first (last hrtf-data))))
    (if (< elevation e-min)
	(format t "; Warning: elevation [~s] less than minimum ~s~%" elevation e-min)
      (if (> elevation e-max)
	  (format t "; Warning: elevation [~s] more than maximum ~s~%" elevation e-max)))
    ;; find the closest elevation impulse response set
    (setf e-pos (position-if #'(lambda (x) (<= (first x) elevation)) hrtf-data :from-end t))
    (if (null e-pos)
	(setf e-pos 0))
    (setf e-pos-gt (+ e-pos 1))
    (if (nth e-pos-gt hrtf-data)
	(if (< (abs (- elevation (first (nth e-pos hrtf-data))))
	       (abs (- (first (nth e-pos-gt hrtf-data)) elevation)))
	    (setf e (nth e-pos hrtf-data))
	  (setf e (nth e-pos-gt hrtf-data)))
      (setf e (nth e-pos hrtf-data)))
    ;; find the closest azimut impulse response
    (setf azimut (mod azimut 360))
    (if (> azimut 180)
	(setf azimut (- azimut 360)))
    (setf a-pos (position-if #'(lambda (x) (<= (first x)(abs azimut))) (second e) :from-end t))
    (if (null a-pos)
	(setf a-pos 0))
    (setf a-pos-gt (+ a-pos 1))
    (if (nth a-pos-gt (second e))
	(if (< (abs (- elevation (first (nth a-pos (second e)))))
	       (abs (- (first (nth a-pos-gt (second e))) elevation)))
	    (setf a (nth a-pos (second e)))
	  (setf a (nth a-pos-gt (second e))))
      (setf a (nth a-pos (second e))))
    (setf impulse-l (second a)
	  impulse-r (third a))
    
    ;; create unit generators
    (let* ((beg (floor (* start-time *srate*)))
	   (in-file (open-input file))
	   (in-file-dur (sound-duration in-file))
	   (end (+ beg (floor (* in-file-dur *srate*))))
	   (filter-l (make-double-float-array hrtf-size))
	   (filter-r (make-double-float-array hrtf-size))
	   ;; create the convolvers
	   (conv-l (make-convolve :input in-file 
				  :fft-size hrtf-size
				  :filter (loop
					    for j from 0 below hrtf-size
					    for k from (if (< azimut 0) impulse-l impulse-r) do
					    (setf (aref filter-l j)
						  (double (aref hrtf-impulses k)))
					    finally (return filter-l))))
	   (conv-r (make-convolve :input in-file
				  :fft-size hrtf-size
				  :filter (loop
					    for j from 0 below hrtf-size
					    for k from (if (< azimut 0) impulse-r impulse-l) do
					    (setf (aref filter-r j)
						  (double (aref hrtf-impulses k)))
					    finally (return filter-r))))
	   ;; create the fir filters for crosstalk cancellation
	   (fir-l (make-fir-filter (length (first hrtf-fir-data))
				   (first hrtf-fir-data)))
	   (fir-r (make-fir-filter (length (second hrtf-fir-data))
				   (second hrtf-fir-data))))
      (run
       (loop for i from beg to end do 
	     (let* ((left (* amplitude (convolve conv-l)))
		    (right (* amplitude (convolve conv-r)))
		    (firl (fir-filter fir-l (+ left right)))
		    (firr (fir-filter fir-r (- left right)))
		    (outl (+ firl firr))
		    (outr (- firl firr)))
	       (outa i (if headphones left outl))
	       (outb i (if headphones right outr))
	       (if *reverb* (out-any  i (* (/ (+ (if headphones left outl) 
						 (if headphones right outr)) 
					    2) reverb-amount)
				      0 *reverb*)))))
      (close-input in-file))))



#|

(with-sound(:channels 2 :statistics t)
  (loop 
      for time from 0 by 0.3
      for az from 0 by 5 below 181 do
	(hrtf time 0.25 "/usr/ccrma/web/html/courses/220b/lectures/5/examples/sounds/knife-11.snd" :azimut az :headphones t)))

|#