This commit is contained in:
5684185+vsariola@users.noreply.github.com
2025-12-30 22:43:35 +02:00
parent f765d75fde
commit 2303e89bbd
10 changed files with 330 additions and 110 deletions

View File

@ -4,34 +4,35 @@ import (
"math"
"math/cmplx"
"github.com/viterin/vek/vek32"
"github.com/vsariola/sointu"
)
type (
SpecAnalyzer struct {
settings SpecSettings
settings SpecAnSettings
broker *Broker
chunker chunker
temp specTemp
}
SpecSettings struct {
Channels SpecChannels
Smooth SpecSmooth
SpecAnSettings struct {
ChnMode SpecChnMode
Smooth SpecSmoothing
Resolution int
}
SpecChannels int
SpecSmooth int
Spectrum []Decibel
SpecResult []Spectrum
SpecChnMode int
SpecSmoothing int
Spectrum [2][]float32
specTemp struct {
spectra []Spectrum
chunk sointu.AudioBuffer
weight []float32 // window weighting function
power [2][]float32
window []float32 // window weighting function
normFactor float32 // normalization factor, to account for the windowing
perm []int // bit-reversal permutation table
tmp []complex128 // temporary buffer for FFT
bitPerm []int // bit-reversal permutation table
tmpC []complex128 // temporary buffer for FFT
tmp1, tmp2 []float32 // temporary buffers for processing
}
)
@ -41,23 +42,35 @@ const (
)
const (
SpecChannelsOff SpecChannels = iota // no spectrum analysis is done to save CPU resources
SpecChannelsCombined // calculate a single combined spectrum for both channels
SpecChannelsSeparated // calculate separate spectrums for left and right channels
SpecChannelsLeft // calculate spectrum only for the left channel
SpecChannelsRight // calculate spectrum only for the right channel
SpecChnModeOff SpecChnMode = iota // no spectrum analysis is done to save CPU resources
SpecChnModeCombine // calculate a single combined spectrum for both channels
SpecChnModeSeparate // calculate separate spectrums for left and right channels
SpecChnModeLeft // calculate spectrum only for the left channel
SpecChnModeRight // calculate spectrum only for the right channel
NumSpecChnModes
)
const (
SpecSmoothSlow SpecSmooth = iota
SpecSmoothMedium
SpecSmoothFast
SpecSmoothingSlow SpecSmoothing = iota
SpecSmoothingMedium
SpecSmoothingFast
NumSpecSmoothing
)
var spectrumSmoothingMap map[SpecSmooth]float32 = map[SpecSmooth]float32{
SpecSmoothSlow: 0.05,
SpecSmoothMedium: 0.2,
SpecSmoothFast: 1.0,
var spectrumSmoothingMap map[SpecSmoothing]float32 = map[SpecSmoothing]float32{
SpecSmoothingSlow: 0.05,
SpecSmoothingMedium: 0.2,
SpecSmoothingFast: 1.0,
}
func NewSpecAnalyzer(broker *Broker) *SpecAnalyzer {
ret := &SpecAnalyzer{broker: broker}
ret.init(SpecAnSettings{
ChnMode: SpecChnModeCombine,
Smooth: SpecSmoothingMedium,
Resolution: 10,
})
return ret
}
func (s *SpecAnalyzer) Run() {
@ -72,37 +85,45 @@ func (s *SpecAnalyzer) Run() {
}
}
func (s *SpecAnalyzer) handleMsg(msg any) {
switch m := msg.(type) {
case SpecSettings:
if s.settings != m {
s.init(m)
func (s *SpecAnalyzer) handleMsg(msg MsgToSpecAn) {
if msg.HasSettings {
s.init(msg.SpecSettings)
}
switch m := msg.Data.(type) {
case *sointu.AudioBuffer:
if s.settings.ChnMode != SpecChnModeOff {
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)})
})
}
case sointu.AudioBuffer:
s.update(m)
s.broker.PutAudioBuffer(m)
default:
// unknown message type; ignore
}
}
func (a *SpecAnalyzer) init(s SpecSettings) {
func (a *SpecAnalyzer) init(s SpecAnSettings) {
s.Resolution = min(max(s.Resolution, SpecResolutionMin), SpecResolutionMax)
a.settings = s
n := 1 << s.Resolution
a.temp = specTemp{
spectra: make([]Spectrum, 0),
chunk: make(sointu.AudioBuffer, n),
weight: make([]float32, n),
perm: make([]int, n),
tmp: make([]complex128, n),
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.weight[i] = w
a.temp.window[i] = w
a.temp.normFactor += w
// initialize the bit-reversal permutation table
a.temp.perm[i] = i
a.temp.bitPerm[i] = i
}
// compute the bit-reversal permutation
for i, j := 1, 0; i < n; i++ {
@ -111,36 +132,86 @@ func (a *SpecAnalyzer) init(s SpecSettings) {
j ^= bit
}
j ^= bit
if i < j {
a.temp.perm[i], a.temp.perm[j] = a.temp.perm[j], a.temp.perm[i]
a.temp.bitPerm[i], a.temp.bitPerm[j] = a.temp.bitPerm[j], a.temp.bitPerm[i]
}
}
}
func (sd *spectrumDetector) update(input []float32) {
c := sd.tmp
for i := range sd.tmp {
p := sd.perm[i]
c[i] = complex(float64(input[p]*sd.window[p]), 0)
func (s *SpecAnalyzer) update(buf sointu.AudioBuffer) *Spectrum {
ret := s.broker.GetSpectrum()
switch s.settings.ChnMode {
case SpecChnModeLeft:
s.process(buf, 0)
ret[0] = append(ret[0], s.temp.power[0]...)
case SpecChnModeRight:
s.process(buf, 1)
ret[1] = append(ret[1], s.temp.power[1]...)
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 SpecChnModeCombine:
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])
vek32.MulNumber_Inplace(ret[0], 0.5)
}
// convert to decibels
for c := range 2 {
vek32.MaximumNumber_Inplace(ret[c], 1e-8)
vek32.MinimumNumber_Inplace(ret[c], 1e8)
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] = 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 l := 2; l <= n; l <<= 1 {
ang := 2 * math.Pi / float64(l)
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 += l {
for i := 0; i < n; i += len {
w := complex(1, 0)
for j := 0; j < l/2; j++ {
for j := 0; j < len/2; j++ {
u := c[i+j]
v := c[i+j+l/2] * w
v := c[i+j+len/2] * w
c[i+j] = u + v
c[i+j+l/2] = u - v
c[i+j+len/2] = u - v
w *= wlen
}
}
}
for i := range input {
a := cmplx.Abs(c[i])
power := float32(a*a) / sd.normFactor
sd.spectrum[i] += sd.alpha * (power - sd.spectrum[i])
// 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 := spectrumSmoothingMap[sd.settings.Smooth]
vek32.MulNumber_Inplace(t2, alpha)
vek32.Add_Inplace(sd.temp.power[channel], t2)
}