(bind-func static hilbert_c
"n = filter order (should be 2^n + 1)"
(lambda (n)
(let ((i 0)
(h_n:SAMPLE* (alloc n))
(x_n:SAMPLE* (alloc n))
(x_ptr 0))
;; h[i] = (2/(i*pi))*sin^2((i*pi)/2) for i=-n/2,...,n/2
(dotimes (i n (/ (* n -1) 2))
(pset! h_n
(+ i (/ n 2))
(* (/ 2.0 (* (convert i) SPI))
(pow (sin (/ (* (convert i) SPI) 2.0)) 2.0))))
;; h[0] = 0
(pset! h_n (+ 1 (/ n 2)) 0.0)
(lambda (x)
(pset! x_n x_ptr x)
(let ((out 0.0))
(dotimes (i n)
(set! out (+ out (* (pref h_n i)
(pref x_n (% (+ i x_ptr (- n 1)) n))))))
(% (set! x_ptr (+ x_ptr 1)) n)
out)))))