;; STFT
;;
;; r, n and k should be powers of 2
;;
;; r is the rate at which you want forward fft's to be calculated
;; n is the size of the fft window 'unpadded'
;; k is the size of the fft window 'padded'
;;
;; r <= n
;; k >= n
(bind-func stft_c 10000000
(lambda (r:i64 n:i64 k:i64)
(if (< k n) (set! k n))
(if (> r n) (set! r n))
(let ((buf:SAMPLE* (alloc n))
(sample_size_in_bytes 4) ;; assuming float SAMPLE
(spectrum:SComplex* (alloc k))
(pad:SAMPLE* (alloc k))
(segments (/ n r))
(pad_offset (dtoi64 (- (* .5 (i64tod k)) (* .5 (i64tod n)))))
(padbuf (cast (pref-ptr pad pad_offset) i8*))
(f_fftr:FFTR_Config* (fftr_config k))
(n2 (* n 2)) (nhalf (/ n 2)) (i 0) (idx (- segments 1)) (t 0) (t2 0))
(lambda (x)
(if (= t2 0)
(let ((start (* idx r))
(samps (- n start)))
;; copy 'n' samples to padded buffer
(memcpy padbuf (convert (pref-ptr buf start)) (* sample_size_in_bytes samps))
(memcpy padbuf (convert (pref-ptr buf 0)) (* sample_size_in_bytes start))
;; apply windowing function
(window_hamming (cast padbuf SAMPLE*) n)
;; forward fft (pad is padbuf)
(fft pad (convert spectrum) f_fftr)
;; inc index
(set! idx (% (+ idx 1) segments))
))
(pset! buf t x)
(set! t (+ t 1))
(if (= t n) (set! t 0))
(set! t2 (+ t2 1))
(if (= t2 r) (set! t2 0))
x))))