(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))))