191 lines
6.4 KiB
Python
191 lines
6.4 KiB
Python
"""Active noise control: the step size of an adaptive filter.
|
|
|
|
The third domain, and the one where the object we have is most exactly the
|
|
object the field argues about.
|
|
|
|
In NLMS the coefficient update is
|
|
|
|
w <- w + mu * e * x / (||x||^2 + eps)
|
|
|
|
and mu = 1 is not an arbitrary scale: it is the step that drives the a
|
|
posteriori error at the current sample exactly to zero. Every practical
|
|
system uses mu < 1, because zeroing the current sample fits the measurement
|
|
noise as well as the signal. So mu is precisely a SAFETY FACTOR on an exact
|
|
step - the same object as the safety factor in ODE step control, and the same
|
|
object as the projection budget in the root-finding safeguard.
|
|
|
|
The whole variable-step-size LMS literature exists because mu is hard to
|
|
choose, and it answers with more constants: Kwong and Johnston's rule has
|
|
four. What this file asks is whether the quantity mu is guessing can be
|
|
measured instead.
|
|
|
|
predicted a posteriori error (1 - mu) * e if the model were exact
|
|
observed a posteriori error d - w_new . x after the update
|
|
|
|
Their ratio says how much of e was signal the filter could actually remove
|
|
and how much was noise it could not. That is the safety factor, measured.
|
|
|
|
The acoustic paths and the noise regimes are taken verbatim from our earlier
|
|
ANC benchmark (branch feature/sem-anc-benchmark of the SEM repository), so
|
|
the environment is the established one and not built to suit the result.
|
|
"""
|
|
from __future__ import annotations
|
|
|
|
import math
|
|
|
|
import numpy as np
|
|
|
|
FS = 8000
|
|
PRIMARY_PATH_LEN = 64
|
|
SECONDARY_PATH_LEN = 32
|
|
FILTER_LEN = 64
|
|
|
|
|
|
# ---------------------------------------------------------------- paths ---
|
|
def make_primary_path(rng):
|
|
h = np.zeros(PRIMARY_PATH_LEN)
|
|
delay = rng.integers(5, 10)
|
|
h[delay] = 0.8 + 0.2 * rng.standard_normal()
|
|
for k in range(3):
|
|
idx = delay + rng.integers(3, 15)
|
|
if idx < PRIMARY_PATH_LEN:
|
|
h[idx] = (0.3 - 0.1 * k) * rng.standard_normal()
|
|
for i in range(delay + 5, PRIMARY_PATH_LEN):
|
|
h[i] += 0.05 * math.exp(-0.1 * (i - delay)) * rng.standard_normal()
|
|
return h
|
|
|
|
|
|
def make_secondary_path(rng):
|
|
h = np.zeros(SECONDARY_PATH_LEN)
|
|
delay = rng.integers(2, 5)
|
|
h[delay] = 0.9 + 0.1 * rng.standard_normal()
|
|
for k in range(2):
|
|
idx = delay + rng.integers(2, 8)
|
|
if idx < SECONDARY_PATH_LEN:
|
|
h[idx] = (0.2 - 0.05 * k) * rng.standard_normal()
|
|
for i in range(delay + 3, SECONDARY_PATH_LEN):
|
|
h[i] += 0.03 * math.exp(-0.15 * (i - delay)) * rng.standard_normal()
|
|
return h
|
|
|
|
|
|
# ------------------------------------------------------------- regimes ---
|
|
def gen_tonal(n, rng, f0=200.0):
|
|
t = np.arange(n) / FS
|
|
return np.sin(2 * math.pi * f0 * t + rng.uniform(0, 2 * math.pi))
|
|
|
|
|
|
def gen_broadband(n, rng):
|
|
pink = np.zeros(n)
|
|
rows = np.zeros(5)
|
|
for i in range(n):
|
|
for j in range(5):
|
|
if i % (2 ** j) == 0:
|
|
rows[j] = rng.standard_normal()
|
|
pink[i] = np.sum(rows)
|
|
return pink / (np.std(pink) + 1e-12)
|
|
|
|
|
|
def gen_multitonal(n, rng, f0=100.0):
|
|
t = np.arange(n) / FS
|
|
ph = rng.uniform(0, 2 * math.pi, 3)
|
|
sig = (np.sin(2 * math.pi * f0 * t + ph[0])
|
|
+ 0.5 * np.sin(2 * math.pi * 2 * f0 * t + ph[1])
|
|
+ 0.25 * np.sin(2 * math.pi * 3 * f0 * t + ph[2]))
|
|
return sig / (np.max(np.abs(sig)) + 1e-12)
|
|
|
|
|
|
def gen_silence(n, rng):
|
|
return rng.standard_normal(n) * 0.01
|
|
|
|
|
|
def gen_impulse_train(n, rng):
|
|
sig = np.zeros(n)
|
|
for _ in range(rng.integers(5, 12)):
|
|
pos = rng.integers(0, n)
|
|
width = rng.integers(10, 40)
|
|
amp = rng.uniform(1.0, 3.0) * rng.choice([-1, 1])
|
|
end = min(pos + width, n)
|
|
sig[pos:end] = amp * np.exp(-5.0 * np.arange(end - pos) / width)
|
|
return sig
|
|
|
|
|
|
def scen_stable(n, rng):
|
|
seg = n // 4
|
|
return np.concatenate([gen_tonal(seg, rng), gen_broadband(seg, rng),
|
|
gen_multitonal(seg, rng),
|
|
gen_silence(n - 3 * seg, rng)])
|
|
|
|
|
|
def scen_gradual(n, rng):
|
|
a = np.linspace(0, 1, n)
|
|
return (1 - a) * gen_tonal(n, rng) + a * gen_broadband(n, rng)
|
|
|
|
|
|
def scen_abrupt(n, rng):
|
|
out = np.zeros(n)
|
|
gens = [gen_tonal, gen_broadband, gen_multitonal]
|
|
pos = 0
|
|
while pos < n:
|
|
seg = min(int(rng.integers(n // 20, n // 8)), n - pos)
|
|
out[pos:pos + seg] = gens[rng.integers(0, 3)](seg, rng)
|
|
pos += seg
|
|
return out
|
|
|
|
|
|
def scen_cabin(n, rng):
|
|
engine = gen_multitonal(n, rng, f0=80.0)
|
|
road = 0.6 * gen_broadband(n, rng)
|
|
imp = 0.4 * gen_impulse_train(n, rng)
|
|
s = engine + road + imp
|
|
return s / (np.max(np.abs(s)) + 1e-12)
|
|
|
|
|
|
def scen_impulsive(n, rng):
|
|
return gen_impulse_train(n, rng) + 0.2 * gen_broadband(n, rng)
|
|
|
|
|
|
SCENARIOS = [("stable regimes", scen_stable), ("gradual transition", scen_gradual),
|
|
("abrupt switching", scen_abrupt), ("car cabin", scen_cabin),
|
|
("impulsive", scen_impulsive)]
|
|
|
|
|
|
# --------------------------------------------------------- the filter ---
|
|
def run_nlms(noise, prim, sec, rule, snr_db=30.0, seed=0):
|
|
"""One NLMS/FxLMS run. Everything is fixed except `rule`, which returns
|
|
the step size for the current sample."""
|
|
rng = np.random.default_rng(seed)
|
|
n = len(noise)
|
|
d = np.convolve(noise, prim)[:n] # what reaches the ear
|
|
xf = np.convolve(noise, sec)[:n] # filtered reference
|
|
p = float(np.mean(d * d))
|
|
d = d + rng.standard_normal(n) * math.sqrt(p / (10 ** (snr_db / 10.0)))
|
|
|
|
w = np.zeros(FILTER_LEN)
|
|
buf = np.zeros(FILTER_LEN)
|
|
err = np.zeros(n)
|
|
rule.reset()
|
|
for i in range(n):
|
|
buf[1:] = buf[:-1]
|
|
buf[0] = xf[i]
|
|
y = float(np.dot(w, buf))
|
|
e = d[i] - y
|
|
err[i] = e
|
|
nrm = float(np.dot(buf, buf)) + 1e-6
|
|
mu = rule.step(e, buf, nrm)
|
|
if mu != 0.0:
|
|
w += mu * e * buf / nrm
|
|
# the a posteriori error, which is what the rule needs to learn
|
|
# how much of e it was able to remove
|
|
rule.observe(e, d[i] - float(np.dot(w, buf)), mu)
|
|
else:
|
|
rule.observe(e, e, 0.0)
|
|
# noise reduction: power of the residual against the power of what would
|
|
# have reached the ear untreated, measured on the second half so the
|
|
# transient is not counted
|
|
half = n // 2
|
|
p_in = float(np.mean(d[half:] ** 2))
|
|
p_out = float(np.mean(err[half:] ** 2))
|
|
if p_out <= 0.0 or not math.isfinite(p_out):
|
|
return float("nan")
|
|
return 10.0 * math.log10(p_in / p_out)
|