(require 'asdf)
(mapc 'require '(#:alexandria #:bordeaux-fft))

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

(defvar *array* 
 (make-array 4096 :initial-contents
            (loop for x below 4096 collecting 0.0)))

(loop for x below 4096 do
 (setf (aref *array* x)
      (+ (wave-3 x)
        (cond ((< 100 x 300) (wave-1 x))
             ((< 3700 x 4000) (wave-2 x))
             (t 0)))))

(defvar *kernel* 
 (make-array 4096 :initial-contents
            (loop for x below 4096 collecting ; just change to wave-2
             (if (< x 200) (wave-1 x) 0.0)))) ; for 2nd kernel

(defvar *convolution-1*
 (make-array 4096 :initial-contents
            (loop for x below 4096 collecting #C(0.0 0.0))))

(let ((farray (bordeaux-fft:sfft *array*))
      (fkernel (bordeaux-fft:sfft *kernel*))
      (fconv   (make-array 4096))
      (auto    (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 o across (bordeaux-fft:sifft auto)
               for p across (bordeaux-fft:sifft fconv) do
               (setf (aref *convolution-1* n) p))))

(with-open-file (out "to-graph.txt" :direction :output :if-exists :supersede)
 (loop for c across *convolution-1*
  do (print (realpart c) out)
  finally (terpri)))

(ext:quit)