;; uniform partitioned convolution in frq domain
;; expects mono aif/wav audio file at samplerate
(bind-func creverb_c
(lambda (filename:i8* len:i64)
(let ((ir_len_unpadded (sf_samples_from_file filename)) ;; get size of ir file (assumes mono!!)
(ir_len (+ ir_len_unpadded (- len (% ir_len_unpadded len)))) ;; pad to whole multiple of len
(parts (/ ir_len len)) ;; how many partitions
(po 0) ;; partition offset
(irs:SComplex* (zalloc (* 2 ir_len))) ;; storage for all IR FFT 'partitions'
(ffts:SComplex* (zalloc (* 2 ir_len))) ;; storage for all total forward ffts
(padlgth (* len 2)) ;; padlgth is for convolution buffers
(aorb:i1 #t) ;; #t for a #f for b
(irtaila:SAMPLE* (zalloc len)) ;; buffer ir tail
(irtailb:SAMPLE* (zalloc len))
(f_fft:FFT_Config* (fft_config padlgth)) ;; forward complex
(i_fft:IFFT_Config* (ifft_config padlgth)) ;; inverse complex
(f_fftr:FFTR_Config* (fftr_config padlgth)) ;; forward real
(i_fftr:IFFTR_Config* (ifftr_config padlgth)) ;; inverse real
(drya:SAMPLE* (zalloc padlgth)) ;; buffer input signal
(dryb:SAMPLE* (zalloc padlgth))
(t1:i64 0) (i:i64 0) (out:SAMPLE 0.0) (j:i64 0)
(A:SComplex* (zalloc padlgth))
(B:SComplex* (zalloc padlgth))
(a:SAMPLE* (zalloc padlgth))
(t2:double 0.0)
(output:SAMPLE* (zalloc len))) ;; buffer output
(printf "IR:(%.1f seconds) Partitions(size:%lld num:%lld)\n" (ftod (/ (convert ir_len) (convert SAMPLE_RATE))) len parts)
;; partition and store IR spectra (circlar convolution so padded)
(dotimes (i parts)
(memset (cast A i8*) 0 (convert (* 8 padlgth))) ;; zero pad, assumes Complexf (8 bytes per struct)
(sf_read_file_into_buffer filename output (* i len) len #f)
(Complex_bufferize output A len)
(fft A (pref-ptr irs (* i padlgth)) f_fft))
(lambda (in:SAMPLE dry:SAMPLE wet:SAMPLE)
(set! t2 (clock_clock))
;; store 'input' data which is 'ahead of out' by len samples
(pset! (if aorb drya dryb) t1 in)
;; set output which is len samples behind 'in'
(set! out ;; out is always delayed by len
(+ (* dry (pref (if aorb dryb drya) t1)) ;; playback 'delayed' input
(* 4.0 wet ;; wet output signal
(+ (pref output t1) ;; overlap-add current-output with ...
(pref (if aorb irtaila irtailb) t1))))) ;; delayed-irtail
;; increment time
(set! t1 (+ t1 1))
;; if we have buffered len samples then run convolution
(if (= t1 len)
(let ((_fft:SComplex* (pref-ptr ffts (* po padlgth)))
(_ir:SComplex* null))
;; forward FFT of incoming signal
(fft (if aorb drya dryb) B f_fftr)
;; store FFT to use 'now' and also for 'delayed' use
(memcpy (cast _fft i8*) (cast B i8*) (convert (* padlgth 8)))
;; run convolution over all partions
(dotimes (i parts)
(set! j (% (+ parts (- po i)) parts))
(set! _fft (pref-ptr ffts (* j padlgth)))
(set! _ir (pref-ptr irs (* i padlgth)))
;; '*' and '+' are overloaded for buffers of Complex{f,d}
(Complex_multiplication_bybuf _fft _ir A padlgth)
(Complex_addition_bybuf A B B padlgth))
;; after convolution is complete back to time domain!
(ifft B a i_fftr)
(let ((scale:SAMPLE (/ 1.0 (convert (* len parts))))
(tail (if aorb irtaila irtailb)))
(dotimes (i len)
(pset! output i (* (pref a i) scale)))
(dotimes (i len)
(pset! tail i (* (pref a (+ len i)) scale))))
;; reset everything for next cycle
(set! aorb (if aorb #f #t)) ;; swap buffers
(set! po (% (+ po 1) parts))
(set! t1 0)))
;; (if (= t1 0)
;; (println (- (clock_clock) t2)))
;; simply return out
out))))