FFTW reverb with half-complex format

3 files changed, 47 insertions, 17 deletions

diff --git a/fftw.bqn b/fftw.bqn
index 55f0cbe..1808db7 100644
--- a/fftw.bqn
+++ b/fftw.bqn
@@ -1,4 +1,4 @@
-# FFTW-based replacement for fft.bqn
+# FFTW-based replacement for fft.bqn (unused; see reverb.bqn)
 
 plan ← "*:i32"
 
diff --git a/mix.bqn b/mix.bqn
index 5aafdbd..4d0c7e6 100644
--- a/mix.bqn
+++ b/mix.bqn
@@ -1,6 +1,5 @@
 "mix.bqn takes a single option namespace, or no arguments" ! 1≥≠•args
 o ← ≠◶⟨•Import∘"options.bqn", ⊑⟩ •args
-fft ← o •Import "fft.bqn"
 
 # List of length |𝕩 that starts at 0 and stops just before 1,
 # reversed if 𝕩<0.
@@ -45,21 +44,8 @@ _crossfade ⇐ {
 # Apply reverb impulse response (IR) 𝕨 to 𝕩.
 # Argument rows are extended, giving result length (𝕩+○≠𝕨)-1 (for lists).
 # Equivalent to (+´¨+⌜○↕○≠⊔×⌜)⎉1 except for numeric precision.
-Reverb ⇐ {
-  M ← -˝∘× ≍ +˝∘×⟜⌽ # Complex multiplication
-  lw‿lx ← (¯1⊑≢)¨ 𝕨‿𝕩
-  ! 0<lw
-  # Use the overlap-add method.
-  o ← lw-1          # Overlap length
-  n ← ⌈⌾(2⋆⁼⊢) 3×o  # Window length, including overlap
-  l ← n-o           # Without overlap
-  k ← lx+o          # Result length
-  k0← ⌈⌾(÷⟜l) k     # Rounded up
-  𝕨 {
-    CW ← ⊏ · (FFT n↑𝕨)⊸M⌾FFT n⊸↑
-    {t←0 ⋄ k↑⥊ {r‿s←(-o)(t⊸+⌾(o⊸↑)∘↓⋈↑)CW𝕩⋄t↩s⋄r}˘ ∘‿l⥊k0↑𝕩}⎉1 𝕩
-  }⎉(1≍1+0⌈-˜○=) 𝕩
-}
+reverb ⇐ •Import "reverb.bqn"
+
 # Apply stereo reverb by keeping early IR sides separate and later
 # blending them together; can avoid some balance issues
 ReverbStereoMix ⇐ {
diff --git a/reverb.bqn b/reverb.bqn
new file mode 100644
index 0000000..867ea26
--- /dev/null
+++ b/reverb.bqn
@@ -0,0 +1,44 @@
+# Reverb, using FFTW if available and BQN-based FFT if not
+
+Reverb ← {
+  lw‿lx ← (¯1⊑≢)¨ 𝕨‿𝕩
+  ! 0<lw
+  # Use the overlap-add method.
+  o ← lw-1          # Overlap length
+  n ← ⌈⌾(2⋆⁼⊢) 3×o  # Window length, including overlap
+  l ← n-o           # Without overlap
+  k ← lx+o          # Result length
+  k0← ⌈⌾(÷⟜l) k     # Rounded up
+  𝕨 {
+    CW ← (n↑𝕨) _rev1
+    {t←0 ⋄ k↑⥊ {r‿s←(-o)(t⊸+⌾(o⊸↑)∘↓⋈↑)CW𝕩⋄t↩s⋄r}˘ ∘‿l⥊k0↑𝕩}⎉1 𝕩
+  }⎉(1≍1+0⌈-˜○=) 𝕩
+}
+
+# Convolve 𝕗 and 𝕗≠⊸↑𝕩, assuming length of 𝕗 is a power of 2
+rev1 ← {𝕊
+  # Use the half-complex form of FFTW
+  # Converts e.g. 8 reals to r0,r1,r2,r3,r4,i3,i2,i1 and back
+  pl ← ">" ∾ plan ← "*:i32"
+  Fn ← "/usr/lib/libfftw3.so.3"⊸•FFI
+  createPlan  ← Fn plan‿"fftw_plan_r2r_1d"‿"i32"‿"*f64"‿"&f64"‿"i32"‿"i32"
+  destroyPlan ← Fn ""‿"fftw_destroy_plan"‿pl
+  executePlan ← Fn ""‿"fftw_execute"‿pl
+  FFTR ← { 𝕊⁼𝕩: 1𝕊𝕩 ;
+    plan‿out ← CreatePlan ⟨≠𝕩,𝕩,0¨𝕩,𝕨⊣0,2⋆6⟩
+    ExecutePlan plan
+    DestroyPlan plan
+    out
+  }
+  {
+    mh←-nh←1-˜2÷˜n←≠𝕗
+    M ← (mh(↓-0∾0∾˜⌽∘↑)×) ∾ nh(⌽∘↑+-⊸↑)×⟜⌽○(1⊸↓)  # Scrambled complex ×
+    (n ÷˜ FFTR 𝕗)⊸M⌾FFTR n⊸↑
+  }
+}⎊{𝕊
+  fft ← •Import "fft.bqn"
+  M ← -˝∘× ≍ +˝∘×⟜⌽  # Complex multiplication
+  {⊏ · (FFT 𝕗)⊸M⌾FFT (≠𝕗)⊸↑}
+}@
+
+Reverb