mirror of
https://github.com/vsariola/sointu.git
synced 2026-02-01 05:10:19 -05:00
268 lines
8.0 KiB
Go
268 lines
8.0 KiB
Go
package tracker
|
|
|
|
import (
|
|
"math"
|
|
"math/cmplx"
|
|
|
|
"github.com/viterin/vek/vek32"
|
|
"github.com/vsariola/sointu"
|
|
)
|
|
|
|
type (
|
|
SpecAnalyzer struct {
|
|
settings SpecAnSettings
|
|
broker *Broker
|
|
chunker chunker
|
|
temp specTemp
|
|
}
|
|
|
|
SpecAnSettings struct {
|
|
ChnMode SpecChnMode
|
|
Smooth int
|
|
Resolution int
|
|
}
|
|
|
|
SpecChnMode int
|
|
Spectrum [2][]float32
|
|
|
|
specTemp struct {
|
|
power [2][]float32
|
|
window []float32 // window weighting function
|
|
normFactor float32 // normalization factor, to account for the windowing
|
|
bitPerm []int // bit-reversal permutation table
|
|
tmpC []complex128 // temporary buffer for FFT
|
|
tmp1, tmp2 []float32 // temporary buffers for processing
|
|
}
|
|
|
|
BiquadCoeffs struct {
|
|
b0, b1, b2 float32
|
|
a0, a1, a2 float32
|
|
}
|
|
|
|
SpecAnEnabled Model
|
|
)
|
|
|
|
const (
|
|
SpecResolutionMin = -3
|
|
SpecResolutionMax = 3
|
|
)
|
|
|
|
const (
|
|
SpecSpeedMin = -3
|
|
SpecSpeedMax = 3
|
|
)
|
|
|
|
const (
|
|
SpecChnModeSum SpecChnMode = iota // calculate a single combined spectrum for both channels
|
|
SpecChnModeSeparate // calculate separate spectrums for left and right channels
|
|
NumSpecChnModes
|
|
)
|
|
|
|
func (m *Model) SpecAnEnabled() Bool { return MakeBoolFromPtr(&m.specAnEnabled) }
|
|
|
|
func NewSpecAnalyzer(broker *Broker) *SpecAnalyzer {
|
|
ret := &SpecAnalyzer{broker: broker}
|
|
ret.init(SpecAnSettings{})
|
|
return ret
|
|
}
|
|
|
|
func (m *Model) BiquadCoeffs() (coeffs BiquadCoeffs, ok bool) {
|
|
i := m.d.InstrIndex
|
|
u := m.d.UnitIndex
|
|
if i < 0 || i >= len(m.d.Song.Patch) || u < 0 || u >= len(m.d.Song.Patch[i].Units) {
|
|
return BiquadCoeffs{}, false
|
|
}
|
|
switch m.d.Song.Patch[i].Units[u].Type {
|
|
case "filter":
|
|
p := m.d.Song.Patch[i].Units[u].Parameters
|
|
f := float32(p["frequency"]) / 128
|
|
f *= f
|
|
r := float32(p["resonance"]) / 128
|
|
// The equations for the filter are:
|
|
// s1[n+1] = s1[n] + f*s2[n]
|
|
// h = u - s1[n+1] - r*s2[n]
|
|
// s2[n+1] = s2[n] + f*h = s2[n] + f*(u-s1[n]-f*s2[n]-r*s2[n]) = - f*s1[n]+(1-f*r-f*f)*s2[n] + f*u
|
|
// y_low[n] = s1[n+1], y_band[n] = s2[n+1], y_high[n] = -s1[n+1]-r*s2[n]+u
|
|
// This gives state space representation
|
|
// s(n+1) = A*s(n)+B*u, where A = [1 f;-f 1-f*r-f*f] and B = [0;f]
|
|
// y(n) = C*s(n)+D*u, where
|
|
// C_low = [z 0], C_band = [0 z], C_high = [-z -r], D_high = [1] (note we use those z:s in C to account for those 1 sample time shifts)
|
|
// The transfer function is then H(z) = C*(zI-A)^-1*B + D
|
|
// z*I-A = [z-1 -f; f z+f*r+f*f-1]
|
|
// Calculate (zI-A)^-1*B:
|
|
// (z*I-A)^-1*B = 1/det * [z+f*r+f*f-1 f; -f z-1] * [0;f] = 1/det * f * [f; z-1], where
|
|
// det = (z+f*r+f*f-1)*(z-1)+f^2 = z*z+z*f*r+z*f*f-z-z-f*r-f*f+1+f^2 = z*z + (r*f+f*f-2)*z + 1-f*r = a0*z^2 + a1*z + a2
|
|
// Low: [z 0]*f*[f;z-1] / det = f*f*z / det = b1 * z / det
|
|
// Band: [0 z]*f*[f;z-1] / det = (f*z^2-f*z) / det = (b0*z^2 + b1*z) / det
|
|
// High: [-z -r]*f*[f;z-1] / det + 1 = ((-f*f-r*f)*z+r*f)/det + 1 = ((-f*f-r*f)*z+r*f+det)/det = (z^2-2*z+1)/det = (b0*z^2 + b1*z + b2)/det
|
|
// Negative versions have only b coefficients negated
|
|
var a0 float32 = 1
|
|
var a1 float32 = r*f + f*f - 2
|
|
var a2 float32 = 1 - f*r
|
|
var b0, b1, b2 float32
|
|
b1 += f * f * float32(p["lowpass"])
|
|
b0 += f * float32(p["bandpass"])
|
|
b1 -= f * float32(p["bandpass"])
|
|
b0 += float32(p["highpass"])
|
|
b1 += -2 * float32(p["highpass"])
|
|
b2 += float32(p["highpass"])
|
|
return BiquadCoeffs{a0: a0, a1: a1, a2: a2, b0: b0, b1: b1, b2: b2}, true
|
|
case "belleq":
|
|
f := float32(m.d.Song.Patch[i].Units[u].Parameters["frequency"]) / 128
|
|
band := float32(m.d.Song.Patch[i].Units[u].Parameters["bandwidth"]) / 128
|
|
gain := float32(m.d.Song.Patch[i].Units[u].Parameters["gain"]) / 128
|
|
omega0 := 2 * f * f
|
|
alpha := float32(math.Sin(float64(omega0))) * 2 * band
|
|
A := float32(math.Pow(2, float64(gain-.5)*6.643856189774724))
|
|
u, v := alpha*A, alpha/A
|
|
return BiquadCoeffs{
|
|
b0: 1 + u,
|
|
b1: -2 * float32(math.Cos(float64(omega0))),
|
|
b2: 1 - u,
|
|
a0: 1 + v,
|
|
a1: -2 * float32(math.Cos(float64(omega0))),
|
|
a2: 1 - v,
|
|
}, true
|
|
default:
|
|
return BiquadCoeffs{}, false
|
|
}
|
|
}
|
|
|
|
func (c *BiquadCoeffs) Gain(omega float32) float32 {
|
|
e := cmplx.Rect(1, -float64(omega))
|
|
return float32(cmplx.Abs((complex(float64(c.b0), 0) + complex(float64(c.b1), 0)*e + complex(float64(c.b2), 0)*(e*e)) /
|
|
(complex(float64(c.a0), 0) + complex(float64(c.a1), 0)*e + complex(float64(c.a2), 0)*e*e)))
|
|
}
|
|
|
|
func (s *SpecAnalyzer) Run() {
|
|
for {
|
|
select {
|
|
case <-s.broker.CloseSpecAn:
|
|
close(s.broker.FinishedSpecAn)
|
|
return
|
|
case msg := <-s.broker.ToSpecAn:
|
|
s.handleMsg(msg)
|
|
}
|
|
}
|
|
}
|
|
|
|
func (s *SpecAnalyzer) handleMsg(msg MsgToSpecAn) {
|
|
if msg.HasSettings {
|
|
s.init(msg.SpecSettings)
|
|
}
|
|
switch m := msg.Data.(type) {
|
|
case *sointu.AudioBuffer:
|
|
buf := *m
|
|
l := len(s.temp.window)
|
|
// 50% overlap with the windows
|
|
s.chunker.Process(buf, l, l>>1, func(chunk sointu.AudioBuffer) {
|
|
TrySend(s.broker.ToModel, MsgToModel{Data: s.update(chunk)})
|
|
})
|
|
s.broker.PutAudioBuffer(m)
|
|
default:
|
|
// unknown message type; ignore
|
|
}
|
|
}
|
|
|
|
func (a *SpecAnalyzer) init(s SpecAnSettings) {
|
|
s.Resolution = min(max(s.Resolution, SpecResolutionMin), SpecResolutionMax) + 10
|
|
a.settings = s
|
|
n := 1 << s.Resolution
|
|
a.temp = specTemp{
|
|
power: [2][]float32{make([]float32, n/2), make([]float32, n/2)},
|
|
window: make([]float32, n),
|
|
bitPerm: make([]int, n),
|
|
tmpC: make([]complex128, n),
|
|
tmp1: make([]float32, n),
|
|
tmp2: make([]float32, n),
|
|
}
|
|
for i := range n {
|
|
// Hanning window
|
|
w := float32(0.5 * (1 - math.Cos(2*math.Pi*float64(i)/float64(n-1))))
|
|
a.temp.window[i] = w
|
|
a.temp.normFactor += w
|
|
// initialize the bit-reversal permutation table
|
|
a.temp.bitPerm[i] = i
|
|
}
|
|
// compute the bit-reversal permutation
|
|
for i, j := 1, 0; i < n; i++ {
|
|
bit := n >> 1
|
|
for ; j&bit != 0; bit >>= 1 {
|
|
j ^= bit
|
|
}
|
|
j ^= bit
|
|
|
|
if i < j {
|
|
a.temp.bitPerm[i], a.temp.bitPerm[j] = a.temp.bitPerm[j], a.temp.bitPerm[i]
|
|
}
|
|
}
|
|
}
|
|
|
|
func (s *SpecAnalyzer) update(buf sointu.AudioBuffer) *Spectrum {
|
|
ret := s.broker.GetSpectrum()
|
|
switch s.settings.ChnMode {
|
|
case SpecChnModeSeparate:
|
|
s.process(buf, 0)
|
|
s.process(buf, 1)
|
|
ret[0] = append(ret[0], s.temp.power[0]...)
|
|
ret[1] = append(ret[1], s.temp.power[1]...)
|
|
case SpecChnModeSum:
|
|
s.process(buf, 0)
|
|
s.process(buf, 1)
|
|
ret[0] = append(ret[0], s.temp.power[0]...)
|
|
vek32.Add_Inplace(ret[0], s.temp.power[1])
|
|
}
|
|
// convert to decibels
|
|
for c := range 2 {
|
|
vek32.Log10_Inplace(ret[c])
|
|
vek32.MulNumber_Inplace(ret[c], 10)
|
|
}
|
|
return ret
|
|
}
|
|
|
|
func (sd *SpecAnalyzer) process(buf sointu.AudioBuffer, channel int) {
|
|
for i := range buf { // de-interleave
|
|
sd.temp.tmp1[i] = removeNaNsAndClamp(buf[i][channel])
|
|
}
|
|
vek32.Mul_Inplace(sd.temp.tmp1, sd.temp.window) // apply windowing
|
|
vek32.Gather_Into(sd.temp.tmp2, sd.temp.tmp1, sd.temp.bitPerm) // bit-reversal permutation
|
|
// convert into complex numbers
|
|
c := sd.temp.tmpC
|
|
for i := range c {
|
|
c[i] = complex(float64(sd.temp.tmp2[i]), 0)
|
|
}
|
|
// FFT
|
|
n := len(c)
|
|
for len := 2; len <= n; len <<= 1 {
|
|
ang := 2 * math.Pi / float64(len)
|
|
wlen := complex(math.Cos(ang), math.Sin(ang))
|
|
for i := 0; i < n; i += len {
|
|
w := complex(1, 0)
|
|
for j := 0; j < len/2; j++ {
|
|
u := c[i+j]
|
|
v := c[i+j+len/2] * w
|
|
c[i+j] = u + v
|
|
c[i+j+len/2] = u - v
|
|
w *= wlen
|
|
}
|
|
}
|
|
}
|
|
// take absolute values of the first half, including nyquist frequency but excluding DC
|
|
m := n / 2
|
|
t1 := sd.temp.tmp1[:m]
|
|
t2 := sd.temp.tmp2[:m]
|
|
for i := 0; i < m; i++ {
|
|
t1[i] = float32(cmplx.Abs(c[1+i])) // do not include DC
|
|
}
|
|
// square the amplitudes to get power
|
|
vek32.Mul_Into(t2, t1, t1)
|
|
vek32.DivNumber_Inplace(t2, sd.temp.normFactor*sd.temp.normFactor) // normalize for windowing
|
|
// Since we are using a real-valued FFT, we need to double the values except for Nyquist (and DC, but we don't have that here)
|
|
vek32.MulNumber_Inplace(t2[:m-1], 2)
|
|
// calculate difference to current spectrum and add back, multiplied by smoothing factor
|
|
vek32.Sub_Inplace(t2, sd.temp.power[channel])
|
|
alpha := float32(math.Pow(2, float64(sd.settings.Smooth-SpecSpeedMax)))
|
|
vek32.MulNumber_Inplace(t2, alpha)
|
|
vek32.Add_Inplace(sd.temp.power[channel], t2)
|
|
}
|