(defpackage cold (:export fftw-convolve fft-kernel))
(in-package cold)
;;;;LENGTH is hardcoded to 4096.

(ffi:clines "
#include <stdio.h>
#include <stdlib.h>
#include <complex.h>

#include <fftw3.h>

#define         LENGTH  4096
")

(ffi:clines "
double a[LENGTH];
double k[LENGTH];
fftw_plan k_pin;
fftw_complex *ck;
")

(defparameter *length* (ffi:c-inline () () :int "@(return)=LENGTH;"))

(ffi:c-inline () () nil "
ck = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * LENGTH);
")

(defun fft-kernel () "
"
 (ffi:c-inline () () nil "
k_pin  = fftw_plan_dft_r2c_1d(LENGTH, k, ck, FFTW_ESTIMATE);
fftw_execute(k_pin);
fftw_destroy_plan(k_pin);
"))

(defun gen-kernel (function) "
"
 (let ((x 0) (f 0.0d0))
  (declare (:int x) (:double f))
  (dotimes (s *length*)
  (setf f (* 1.0d0 (funcall function)))
  (ffi:c-inline (x f) (:int :double ) nil "
k[#0] = #1; ")
   (incf x)))
 (fft-kernel))

;;;;INITIALIZE KERNEL TO SOMETHING DUMB
(gen-kernel (lambda () (* 1.0d0 (1- (random 3)))))

(defun kernth (n)
 (ffi:c-inline (n) (:int) :double "
@(return)=k[#0];
"))

(defun set-kernth (n f)
 (ffi:c-inline (n f) (:int :double) :double "
k[#0]=#1;
@(return)=(k[#0]);
"))

(defun anth (n)
 (ffi:c-inline (n) (:int) :double "
@(return)=a[#0];
"))

(defun set-anth (n f)
 (ffi:c-inline (n f) (:int :double) :double "
a[#0]=#1;
@(return)=(a[#0]);
"))


(defun fftw-convolve () "
(fftw-convolve)
"
 (ffi:c-inline () 
               ()
               nil"
fftw_complex *ca;
fftw_plan a_pin;
fftw_plan a_pout;

int idx;


ca = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * LENGTH);

a_pin  = fftw_plan_dft_r2c_1d(LENGTH, a, ca, FFTW_ESTIMATE);
fftw_execute(a_pin);
fftw_destroy_plan(a_pin);

for (idx=0; idx<LENGTH;++idx) { 
        ca[idx] = ((complex double *)ck)[idx] * ((complex double *)ca)[idx]; 
}

a_pout = fftw_plan_dft_c2r_1d(LENGTH, ca, a, FFTW_ESTIMATE);

fftw_execute(a_pout);

fftw_destroy_plan(a_pout);
fftw_free(ca);
"))