mirror of
https://github.com/vsariola/sointu.git
synced 2026-02-03 14:20:17 -05:00
255 lines
7.0 KiB
Go
255 lines
7.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
|
||
}
|
||
)
|
||
|
||
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 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
|
||
// in state-space, the filter has the form:
|
||
// 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 band]
|
||
//
|
||
// 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]
|
||
// Invert it:
|
||
// (z*I-A)^-1 = 1/det * [z+f*r+f*f-1 f; -f z-1], where det = (z-1)*(z+f*r+f*f-1)+f^2 = z^2 + z * (f*r+f*f−2) + 1-f*r
|
||
// (z*I-A)^-1*B = 1/det * [-f*f; f*z-f]
|
||
var a0 float32 = 1
|
||
var a1 float32 = r*f + f*f - 2
|
||
var a2 float32 = 1 - r*f
|
||
var b0, b1, b2 float32
|
||
if p["lowpass"] == 1 {
|
||
b2 = -f * f
|
||
}
|
||
b2 -= f * float32(p["bandpass"])
|
||
b1 += f * float32(p["bandpass"])
|
||
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)
|
||
}
|