;;; -*- 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]

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; 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-size 128)
(defparameter hrtf-data nil)

(defvar hrtf-impulses (make-array (* 2 400 hrtf-size)
				     :element-type 'short-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...~%"))
  (loop 
      with data
      with offset = 0
      for elevation-dir in (directory dir) do
	;; parse the elevation value [elevXXX.DAT]
	(multiple-value-bind (elevation index)
	    (parse-integer (file-namestring elevation-dir) :start 4 :junk-allowed t)
	  ;; each elevation is a directory, process the files inside
	  (loop 
	      with current
	      with azimut
	      with name
	      with impulse-file
	      with impulse-in-a
	      with impulse-in-b
	      for azimut-file in (directory (concatenate 'string (namestring elevation-dir) "/")) do
		(setf name (file-namestring azimut-file))
		;; parse the filename's elevation and azimut [hXXXeYYYa.DAT]
		(multiple-value-bind (elevation index)
		    (parse-integer name :start 1 :junk-allowed t)
		  (setf azimut (parse-integer name :start (+ index 1) :junk-allowed t))
		  (if verbose (format t "reading ~a --> [~s:~s]...~%" name elevation azimut))
		  (setf impulse-file (open-input azimut-file)
			impulse-in-a (make-readin impulse-file 0)
			impulse-in-b (make-readin impulse-file 1))
		  ;; read in the impulse response for the given elevation and azimut
		  (loop for i from 0 below size do
		    (setf (aref hrtf-impulses (+ offset i))(readin impulse-in-a)
			  (aref hrtf-impulses (+ offset size i))(readin impulse-in-b)))
		  (close-input impulse-file)
		  ;; merge each azimut list into the current elevation
		  (setf current (if (not current)
				    (list (list azimut offset (+ offset size)))
				  (lisp: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))
			     (lisp: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))
    (format t "at [~s:~s] ~s:~s~%" elevation azimut (first e) (first 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-float (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-float (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 *reverb*) (if headphones right outr)) 
					    2) reverb-amount))))))
      (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)))

|#