(in-package sdffft)

(defun wave-1 (x)
 (* (sin (* x 2 pi (/ 200)))
   (sin (* x 2 pi (/ 17)))))

(defun wave-2 (x)
 (* (sin (* x 2 pi (/ 300)))
   (sin (* x 2 pi (/ 23)))))

(defun wave-3 (x)
 (* (sin (* x 2 pi (/ 250)))
   (sin (* x 2 pi (/ 37)))))

(defparameter *array-length* 4096)
(defvar *array* (make-array *array-length* :element-type 'double-float))


(defvar *kernel* (make-array *array-length* :element-type 'double-float))
(setf (documentation '*kernel* 'variable) "
Just change wave-1 to wave-2 to try another kernel
")


(defvar *convolution* (make-array *array-length* 
                                 :element-type '(complex double-float)))

(defun eg-0 () "
"
      ;; Initialize *array*
      (loop for x below *array-length* do
           (setf (aref *kernel* x) (+ (wave-3 x)
                                     (cond ((< 100 x 300) (wave-1 x))
                                          ((< 3700 x 4000) (wave-2 x))
                                          (t 0.0D0)))))

      ;; Set *kernel* to 200 entries of #'wav-1
      (loop for x below *array-length* do
           (setf (aref *array* x) (if (< x 200) (wave-1 x) 0.0D0)))

      ;; setf *convolution* the convolution of *array* and *kernel*
      ;; auto (autocorrelation) unused for now
      (let ((farray     (bordeaux-fft:sfft *array*))
            (fkernel    (bordeaux-fft:sfft *kernel*))
            ;; (auto       (make-array 4096))
            (fconv      (make-array 4096)))
          (loop for a across farray
               for k across fkernel
               for x from 0
               do (setf (aref fconv x) (* a k))
               ;; do (setf (aref auto x) (* a a))
               finally (loop for n from 0
                            for p across (bordeaux-fft:sifft fconv)
                            do (setf (aref *convolution* n) (abs p)))))
      
      ;;Write to a file, "to-graph.txt"
      (with-open-file (out #p"to-graph.txt" :direction :output
                                           :if-exists :supersede)
       (loop for c across *convolution*
            do (prin1 (abs c) out)
            do (terpri out))))