feat: implement bell filter unit for equalizing

CHANGELOG.md

@@ -5,6 +5,10 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
 
 ## [Unreleased]
 ### Added
+- belleq unit: a bell-shaped second-order filter for equalization. Belleq unit
+  takes the center frequency, bandwidth (inverse of Q-factor) and gain (+-40
+  dB). Useful for boosting or reducing specific frequency ranges. Kudos to Reaby
+  for the initial implementation!
 - Multithreaded synths: the user can split the patch up to four threads.
   Selecting the thread can be done on the instrument properties pane.
   Multithreading works only on the multithreaded synths, selectable from the CPU

patch.go

@@ -206,6 +206,11 @@ var UnitTypes = map[string]([]UnitParameter){
 		{Name: "stereo", MinValue: 0, MaxValue: 1, CanSet: true, CanModulate: false},
 		{Name: "channel", MinValue: 0, MaxValue: 6, CanSet: true, CanModulate: false, DisplayFunc: arrDispFunc(channelNames[:])}},
 	"sync": []UnitParameter{},
+	"belleq": []UnitParameter{
+		{Name: "stereo", MinValue: 0, MaxValue: 1, CanSet: true, CanModulate: false},
+		{Name: "frequency", MinValue: 0, MaxValue: 128, CanSet: true, CanModulate: true, DisplayFunc: func(v int) (string, string) { return belleqFrequencyDisplay(v) }},
+		{Name: "bandwidth", MinValue: 0, MaxValue: 128, CanSet: true, CanModulate: true, DisplayFunc: func(v int) (string, string) { return belleqBandwidthDisplay(v) }},
+		{Name: "gain", MinValue: 0, Neutral: 64, MaxValue: 128, CanSet: true, CanModulate: true, DisplayFunc: func(v int) (string, string) { return belleqGainDisplay(v) }}},
 }
 
 // compile errors if interface is not implemented.
@@ -266,6 +271,23 @@ func filterFrequencyDispFunc(v int) (string, string) {
 	return strconv.FormatFloat(f, 'f', 0, 64), "Hz"
 }
 
+func belleqFrequencyDisplay(v int) (string, string) {
+	freq := float64(v) / 128
+	p := 2 * freq * freq
+	f := 44100 * p / math.Pi / 2
+	return strconv.FormatFloat(f, 'f', 0, 64), "Hz"
+}
+
+func belleqBandwidthDisplay(v int) (string, string) {
+	p := float64(v) / 128
+	Q := 1 / (4 * p)
+	return strconv.FormatFloat(Q, 'f', 2, 64), "Q"
+}
+
+func belleqGainDisplay(v int) (string, string) {
+	return strconv.FormatFloat(40*(float64(v)/64-1), 'f', 2, 64), "dB"
+}
+
 func compressorTimeDispFunc(v int) (string, string) {
 	alpha := math.Pow(2, -24*float64(v)/128) // alpha is the "smoothing factor" of first order low pass iir
 	sec := -1 / (44100 * math.Log(1-alpha))  // from smoothing factor to time constant, https://en.wikipedia.org/wiki/Exponential_smoothing
@@ -435,6 +457,7 @@ var stackUseMonoStereo = map[string][2]StackUse{
 		{Inputs: [][]int{{0}}, Modifies: []bool{false}, NumOutputs: 1},
 		{},
 	},
+	"belleq": stackUseEffect,
 }
 var stackUseSendNoPop = [2]StackUse{
 	{Inputs: [][]int{{0}}, Modifies: []bool{true}, NumOutputs: 1},

tests/CMakeLists.txt

@@ -155,6 +155,9 @@ regression_test(test_filter_stereo "VCO_SINE;ENVELOPE;FOP_MULP")
 regression_test(test_filter_freqmod "VCO_SINE;ENVELOPE;FOP_MULP;SEND")
 regression_test(test_filter_resmod "VCO_SINE;ENVELOPE;FOP_MULP;SEND")
 
+regression_test(test_belleq "VCO_SINE;ENVELOPE;FOP_MULP")
+regression_test(test_belleq_stereo "VCO_SINE;ENVELOPE;FOP_MULP")
+
 regression_test(test_delay "ENVELOPE;FOP_MULP;PANNING;VCO_SINE")
 regression_test(test_delay_stereo "ENVELOPE;FOP_MULP;PANNING;VCO_SINE")
 regression_test(test_delay_notetracking "ENVELOPE;FOP_MULP;PANNING;NOISE")

tests/test_belleq.yml

@@ -0,0 +1,24 @@
+bpm: 100
+rowsperbeat: 4
+score:
+    rowsperpattern: 16
+    length: 1
+    tracks:
+        - numvoices: 1
+          order: [0]
+          patterns: [[64, 0, 68, 0, 32, 0, 0, 0, 75, 0, 78, 0, 0, 0, 0, 0]]
+patch:
+    - numvoices: 1
+      units:
+        - type: envelope
+          parameters: {attack: 64, decay: 64, gain: 128, release: 72, stereo: 0, sustain: 64}
+        - type: oscillator
+          parameters: {color: 128, detune: 64, gain: 128, lfo: 0, phase: 0, shape: 64, stereo: 0, transpose: 64, type: 1, unison: 0}
+        - type: mulp
+          parameters: {stereo: 0}
+        - type: belleq
+          parameters: {frequency: 64, bandwidth: 64, gain: 96, stereo: 0}
+        - type: pan
+          parameters: {panning: 64, stereo: 0}
+        - type: out
+          parameters: {gain: 128, stereo: 1}

tests/test_belleq_stereo.yml

@@ -0,0 +1,22 @@
+bpm: 100
+rowsperbeat: 4
+score:
+    rowsperpattern: 16
+    length: 1
+    tracks:
+        - numvoices: 1
+          order: [0]
+          patterns: [[64, 0, 68, 0, 32, 0, 0, 0, 75, 0, 78, 0, 0, 0, 0, 0]]
+patch:
+    - numvoices: 1
+      units:
+        - type: envelope
+          parameters: {attack: 64, decay: 64, gain: 128, release: 72, stereo: 1, sustain: 64}
+        - type: oscillator
+          parameters: {color: 128, detune: 64, gain: 128, lfo: 0, phase: 0, shape: 64, stereo: 1, transpose: 64, type: 1, unison: 0}
+        - type: mulp
+          parameters: {stereo: 1}
+        - type: belleq
+          parameters: {frequency: 64, bandwidth: 64, gain: 96, stereo: 1}        
+        - type: out
+          parameters: {gain: 64, stereo: 1}

tracker/presets.go

@@ -464,6 +464,7 @@ var defaultUnits = map[string]sointu.Unit{
 	"compressor": {Type: "compressor", Parameters: map[string]int{"stereo": 0, "attack": 64, "release": 64, "invgain": 64, "threshold": 64, "ratio": 64}},
 	"send":       {Type: "send", Parameters: map[string]int{"stereo": 0, "amount": 64, "voice": 0, "unit": 0, "port": 0, "sendpop": 1}},
 	"sync":       {Type: "sync", Parameters: map[string]int{}},
+	"belleq":     {Type: "belleq", Parameters: map[string]int{"stereo": 0, "frequency": 64, "bandwidth": 64, "gain": 64}},
 }
 
 var defaultInstrument = sointu.Instrument{

vm/compiler/templates/amd64-386/effects.asm

@@ -198,6 +198,68 @@ su_op_filter_skipneghighpass:
 {{end}}
 
 
+{{- if .HasOp "belleq"}}
+;-------------------------------------------------------------------------------
+;   BELLEQ opcode: perform second order bell eq filtering on the signal
+;-------------------------------------------------------------------------------
+;   Mono:   x   ->  eq(x)
+;   Stereo: l r ->  eq(l) eq(r)
+;-------------------------------------------------------------------------------
+{{.Func "su_op_belleq" "Opcode"}}
+{{- if .Stereo "belleq"}}
+    {{.Call "su_effects_stereohelper"}}
+{{- end}}
+    ; Note: we calculate the gain first because su_power needs temp stack and everything here was crafted to stay altogether below max 4 temp stack
+    ; The cost of staying at max 4 stack was a few extra instructions because of stack juggling.
+    ; The bell filter biquad coefficients (see go_synth.go):
+    ;   b0, b1, b2 = 1+u, -2*cos(w), 1-u
+	;   a0, a1, a2 = 1+v, b1, 1-v
+    ; where w=freq*freq, u=alpha*A, v=alpha/A, alpha=sin(w)*2*bandwidth, A=gain. The filter is implemented as:
+    ;  y = (b0*x+s1)/a0 = ((1+u)*x + s1) / (1+v) = (x+u*x+s1)/(1+v)
+    ;  s1' = b1*x - a1*y + s2 = b1*(x-y)+s2 = 2*cos(w)*(y-x)+s2
+	;  s2' = b2*x - a2*y = (1-u)*x-(1-v)*y = x-y-u*x+v*y
+    fld     dword [{{.Input "belleq" "gain"}}]        ; g x
+	{{- .Float 0.5 | .Prepare | indent 4}}
+    fsub    dword [{{.Float 0.5 | .Use}}]           ; g-0.5 x
+    {{- .Float 6.643856189774724 | .Prepare | indent 4}}
+    fmul    dword [{{.Use (.Float 6.643856189774724)}}]   ; (g-0.5)*6.643856189774724 x
+    {{.Call "su_power"}}                            ; A=2^((g-0.5)*6.643856189774724) x
+    fld     dword [{{.Input "belleq" "frequency"}}] ; f A x
+    fmul    st0, st0                                ; f*f A x
+    fadd    st0, st0                                ; w=2*f*f
+    fsincos                                         ; cos(w) sin(w) A x
+    fadd    st0, st0                                ; r=2*cos(w) sin(w) A x
+    fld     dword [{{.Input "belleq" "bandwidth"}}]   ; b r sin(w) A x
+    fadd    st0, st0                                ; 2*b r sin(w) A x
+    fmulp   st2, st0                                ; r alpha=sin(w)*2*b A x
+	fxch    st0, st1                                ; alpha r A x
+    fdivr   st2, st0                                ; alpha r v=alpha/A x
+	fmul    st0, st0                                ; alpha*alpha r v x
+    fdiv    st0, st2                                ; u=alpha*A r v x
+    fld1                                            ; 1 u r v x
+    faddp   st3, st0                                ; u r v+1 x
+    fmul    st0, st3                                ; u*x r v+1 x
+    fld     dword [{{.WRK}}]                        ; s1 u*x r v+1 x
+    fadd    st0, st1                                ; s1+u*x u*x r v+1 x
+    fadd    st0, st4                                ; s1+u*x+x u*x r v+1 x
+    fdiv    st0, st3                                ; y=(s1+u*x+x)/(v+1) u*x r v+1 x
+    {{- .Float 0.5 | .Prepare | indent 4}}
+    fadd    dword [{{.Float 0.5 | .Use}}]           ; add and sub small offset to prevent denormalization
+    fsub    dword [{{.Float 0.5 | .Use}}]           ; See for example: https://stackoverflow.com/questions/36781881/why-denormalized-floats-are-so-much-slower-than-other-floats-from-hardware-arch
+    fmul    st3, st0                                ; y u*x r v*y+y x
+    fsub    st3, st0                                ; y u*x r v*y x
+    fxch    st4, st0                                ; x u*x r v*y y
+    fsubr   st0, st4                                ; y-x u*x r v*y y
+    fmul    st2, st0                                ; y-x u*x r*(y-x) v*y y
+    fsubp   st3, st0                                ; u*x r*(y-x) x-y+v*y y
+    fsubp   st2, st0                                ; r*(y-x) x-y+v*y-u*x y
+    fadd    dword [{{.WRK}}+4]                      ; s2+r*(y-x) x-y+v*y-u*x y
+    fstp    dword [{{.WRK}}]                        ; x-y+v*y-u*x y
+    fstp    dword [{{.WRK}}+4]                      ; y
+    ret
+{{end}}
+
+
 {{- if .HasOp "clip"}}
 ;-------------------------------------------------------------------------------
 ;   CLIP opcode: clips the signal into [-1,1] range

vm/compiler/templates/wasm/effects.wat

@@ -203,6 +203,62 @@
 {{end}}
 
 
+{{- if .HasOp "belleq"}}
+;;-------------------------------------------------------------------------------
+;;   BELLEQ opcode: perform second order bell eq filtering on the signal
+;;-------------------------------------------------------------------------------
+;;   Mono:   x   ->  eq(x)
+;;   Stereo: l r ->  eq(l) eq(r)
+;;-------------------------------------------------------------------------------
+(func $su_op_belleq (param $stereo i32) (local $sinw f32) (local $A f32) (local $u f32) (local $v f32) (local $x f32) (local $y f32) (local $d f32) (local $alpha f32)
+{{- if .Stereo "belleq"}}
+    (call $stereoHelper (local.get $stereo) (i32.const {{div (.GetOp "belleq") 2}}))
+{{- end}}
+    (global.get $WRK)
+    (local.tee $x (call $pop))                              ;; x WRK
+    (f32.mul
+        (call $input (i32.const {{.InputNumber "belleq" "frequency"}}))
+        (call $input (i32.const {{.InputNumber "belleq" "frequency"}}))
+    )
+    (f32.mul (f32.const 2))
+    (local.tee $sinw (call $sin))                           ;; sinw x WRK
+    (call $input (i32.const {{.InputNumber "belleq" "bandwidth"}})) ;; b sinw x WRK
+    (f32.mul (f32.const 2))                                 ;; 2*b sinw x WRK
+    (local.tee $alpha (f32.mul))                            ;; alpha=sinw*2*b x WRK
+    (f32.sub (call $input (i32.const {{.InputNumber "belleq" "gain"}})) (f32.const 0.5)) ;; g-0.5 alpha x WRK
+    (f32.mul (f32.const 6.643856189774724))
+    (local.tee $A (call $pow2))                             ;; A=2^((g-0.5)*6.643856189774724) alpha x WRK
+    (local.tee $u (f32.mul))                                ;; u=A*alpha x WRK
+    ;; Computing (y=x+u*x+s1)/(1+v)
+    (f32.mul (local.get $x))                                ;; u*x x WRK
+    (f32.add)                                               ;; ux+x WRK
+    (f32.load (global.get $WRK))                            ;; s1 ux+x WRK
+    (f32.add)                                               ;; ux+x+s1 WRK
+    ;; Compute v=alpha/A
+    (local.tee $v (f32.div (local.get $alpha) (local.get $A))) ;; v ux+x+s1 WRK
+    (f32.add (f32.const 1))                                 ;; 1+v ux+x+s1 WRK
+    (local.tee $y (f32.div))                                ;; y WRK
+    ;; s1' = 2*cos(w)*(y-x)+s2
+    (f32.sub (local.get $x))                                ;; y-x WRK
+    ;; need to compute cos(w) as sqrt(1-sin(w)^2)
+    (f32.sqrt (f32.sub (f32.const 1) (f32.mul (local.get $sinw) (local.get $sinw)))) ;; cos(w) y-x WRK
+    (f32.mul)                                               ;; cos(w)*(y-x) WRK
+    (f32.mul (f32.const 2))                                 ;; 2*cos(w)*(y-x) WRK
+    (f32.add (f32.load offset=4 (global.get $WRK)))         ;; s2+2*cos(w)*(y-x) WRK
+    (f32.store)                                             ;; s1'=s2+2*cos(w)*(y-x)
+    ;; s2' = x-y+v*y-u*x
+    (global.get $WRK)
+    (f32.sub (local.get $x) (local.get $y))                 ;; x-y WRK
+    (f32.mul (local.get $v) (local.get $y))
+    (f32.mul (local.get $u) (local.get $x))
+    (f32.sub)                                               ;; v*y-u*x x-y WRK
+    (f32.add)                                               ;; v*y-u*x+x-y WRK
+    (f32.store offset=4)
+    (call $push (local.get $y))
+)
+{{end}}
+
+
 {{- if .HasOp "clip"}}
 ;;-------------------------------------------------------------------------------
 ;;   CLIP opcode: clips the signal into [-1,1] range

vm/go_synth.go

@@ -594,6 +594,25 @@ func (s *GoSynth) Render(buffer sointu.AudioBuffer, maxtime int) (samples int, r
 				if stereo {
 					stack = append(stack, gain)
 				}
+			case opBelleq:
+				// Bell-shaped peaking filter equations based on https://shepazu.github.io/Audio-EQ-Cookbook/audio-eq-cookbook.html:
+				//   alpha = sin(omega0)/(2*Q) where omega0 determines the angular frequency of the peak and Q is the Q-factor
+				//   A = sqrt(10^(dBgain/20)) = 10^(dBgain/40) where dbGain determines the gain at the peak
+				//   b0 = 1 + alpha*A, b1 = -2*cos(omega0), b2 = 1 - alpha*A,
+				//   a0 = 1 + alpha/A, a1 = -2*cos(omega0), a2 = 1 - alpha/A are the biquad filter coefficients
+				omega0 := 2 * params[0] * params[0]                                // square the omega to have a bit more values mapping to bass frequencies
+				alpha := float32(math.Sin(float64(omega0))) * 2 * params[1]        // Q=1/(4*(p/128)) gives a range of Q = 0.25 ... 32
+				A := float32(math.Pow(2, float64(params[2]-.5)*6.643856189774724)) // +-40 dB, reusing same constant as dbgain unit
+				u, v := alpha*A, alpha/A
+				b0, b1, b2 := 1+u, -2*float32(math.Cos(float64(omega0))), 1-u
+				a0, a1, a2 := 1+v, b1, 1-v
+				for i := range channels { // biquad filter in transposed direct from II (https://en.wikipedia.org/wiki/Digital_biquad_filter)
+					x := stack[l-1-i]
+					y := (b0*x + unit.state[i]) / a0 // the biquad was not in normalized form, so we need to divide by a0
+					unit.state[i] = b1*x - a1*y + unit.state[2+i]
+					unit.state[2+i] = b2*x - a2*y
+					stack[l-1-i] = y
+				}
 			case opSync:
 				break
 			default:

vm/go_synth_test.go

@@ -109,6 +109,7 @@ var defaultUnits = map[string]sointu.Unit{
 	"compressor": {Type: "compressor", Parameters: map[string]int{"stereo": 0, "attack": 64, "release": 64, "invgain": 64, "threshold": 64, "ratio": 64}},
 	"send":       {Type: "send", Parameters: map[string]int{"stereo": 0, "amount": 128, "voice": 0, "unit": 0, "port": 0, "sendpop": 1}},
 	"sync":       {Type: "sync", Parameters: map[string]int{}},
+	"belleq":     {Type: "belleq", Parameters: map[string]int{"stereo": 0, "freq": 64, "bandwidth": 64, "gain": 96}},
 }
 
 var defaultInstrument = sointu.Instrument{

vm/opcodes.go

@@ -5,34 +5,35 @@ const (
 	opAdd        = 1
 	opAddp       = 2
 	opAux        = 3
-	opClip       = 4
+	opBelleq     = 4
-	opCompressor = 5
+	opClip       = 5
-	opCrush      = 6
+	opCompressor = 6
-	opDbgain     = 7
+	opCrush      = 7
-	opDelay      = 8
+	opDbgain     = 8
-	opDistort    = 9
+	opDelay      = 9
-	opEnvelope   = 10
+	opDistort    = 10
-	opFilter     = 11
+	opEnvelope   = 11
-	opGain       = 12
+	opFilter     = 12
-	opHold       = 13
+	opGain       = 13
-	opIn         = 14
+	opHold       = 14
-	opInvgain    = 15
+	opIn         = 15
-	opLoadnote   = 16
+	opInvgain    = 16
-	opLoadval    = 17
+	opLoadnote   = 17
-	opMul        = 18
+	opLoadval    = 18
-	opMulp       = 19
+	opMul        = 19
-	opNoise      = 20
+	opMulp       = 20
-	opOscillator = 21
+	opNoise      = 21
-	opOut        = 22
+	opOscillator = 22
-	opOutaux     = 23
+	opOut        = 23
-	opPan        = 24
+	opOutaux     = 24
-	opPop        = 25
+	opPan        = 25
-	opPush       = 26
+	opPop        = 26
-	opReceive    = 27
+	opPush       = 27
-	opSend       = 28
+	opReceive    = 28
-	opSpeed      = 29
+	opSend       = 29
-	opSync       = 30
+	opSpeed      = 30
-	opXch        = 31
+	opSync       = 31
+	opXch        = 32
 )
 
-var transformCounts = [...]int{0, 0, 1, 0, 5, 1, 1, 4, 1, 5, 2, 1, 1, 0, 1, 0, 1, 0, 0, 2, 6, 1, 2, 1, 0, 0, 0, 1, 0, 0, 0}
+var transformCounts = [...]int{0, 0, 1, 3, 0, 5, 1, 1, 4, 1, 5, 2, 1, 1, 0, 1, 0, 1, 0, 0, 2, 6, 1, 2, 1, 0, 0, 0, 1, 0, 0, 0}