#!/usr/bin/env python3 """soothe2 поведенческий симулятор v5 — STFT-детектор (архитектура фазы 3-6). Модель (верифицировано рендерами фазы 4-6): red(f) = base_red(level_f) + 6.02 * min(1, |sens_i|/12) * H(f; fc_i, Qeff_i) * sat(level) H(f; fc, Qeff) = 1/sqrt(1 + (Qeff*A)^2), A = f/fc - fc/f Qeff = 1.54 * q^1.33 (измерено: q=1->1.54, 2->3.32, 3->5.64, 6->16.8) sat(level) = 1/(1+exp(-0.117*(level_db - (-32.99)))) (boost насыщается к 6.02 на высоком уровне) base_red(level): 2.22/3.36/4.85/5.71/6.65/7.65/9.80 @ -30/-24/-18/-15/-12/-9/-3 dBFS (neut, depth0.864) Синтез нотча: per detected peak - спектральная маска на кадре FFT: X[k] *= 1 - g_peak * bell(k; f_peak, Q_notch), g_peak = 1 - 10^(-red(f_peak)/20) Q_notch = 1.15 * sharpness (измерено psh_{1,3,5,10}: Q ~11/6/3/1.2) Детектор: окно 2048/hop 512 Hann; |X| на кадре -> лок. луднесс -> red(f) -> пики (лок. макс. выше порога; selectivity = min-peak-spacing + кол-во). Env attack/release на g_peak. """ import numpy as np import wave SR = 44100 # ---- измеренные LUT ---- _NLEVEL_DB = np.array([-30.02, -24.02, -18.02, -15.02, -12.02, -9.02, -3.05]) _NLEVEL_RED = np.array([2.22, 3.36, 4.85, 5.71, 6.65, 7.65, 9.80]) _DEPTH_MULT = {0.5: 0.56, 0.864: 1.0, 1.0: 1.03, 2.0: 1.24, 3.0: 1.36, 5.0: 1.51, 10.0: 1.58, 20.0: 1.60} # base_red * mult(depth) (грубо, калибровать) _DEPTH_XS = np.array([0.5, 0.864, 1.0, 2.0, 3.0, 5.0, 10.0, 20.0]) _DEPTH_YS = np.array([0.56, 1.0, 1.03, 1.24, 1.36, 1.51, 1.58, 1.60]) def base_red(level_db, depth=0.864): b = np.interp(level_db, _NLEVEL_DB, _NLEVEL_RED) return b * float(np.interp(depth, _DEPTH_XS, _DEPTH_YS)) def sat_level(level_db): return 1.0 / (1 + np.exp(-0.117 * (level_db - (-32.99)))) def qeff(q): return 1.54 * q ** 1.33 def bell_h(f, fc, q): a = f / fc - fc / f qe = qeff(q) return 1.0 / np.sqrt(1 + (qe * a) ** 2) def eq_weight_db(freqs, bands, sens_global=0.0): """Сумма весов по включённым полосам (on=1). bands: list of dict(freq,q,sens,on).""" w = np.zeros_like(freqs, dtype=float) for b in bands: if not b.get('on', 0): continue s = b.get('sens', 0.0) if abs(s) < 1e-9: continue coeff = 6.02 * min(1.0, abs(s) / 12.0) h = bell_h(freqs, b.get('freq', 500.0), b.get('q', 1.0)) w = w + coeff * h return w _QSHARP = (np.array([1.0, 3.0, 5.0, 10.0]), np.array([2.0, 6.0, 10.0, 12.0])) def q_notch_from_sharp(sharp): """Базовый Q нотча от sharpness (умеренный драйв).""" sh = max(float(sharp), 0.05) return float(np.interp(sh, _QSHARP[0], _QSHARP[1], left=2.0, right=12.0)) def notch_bell(freqs, f_peak, qn): """Колокол нотча (gain-форма), peak 1 на f_peak, асимметричный BP-форм.""" a = freqs / f_peak - f_peak / freqs return 1.0 / np.sqrt(1 + (qn * a) ** 2) def stft_frames(x, nfft=2048, hop=512): w = np.hanning(nfft) n = len(x) nf = max(0, (n - nfft) // hop + 1) fr = np.fft.rfftfreq(nfft, 1 / SR) return fr, w, nf, hop def detect_peaks(red, freqs, spacing_hz, min_red_db=3.0, max_peaks=16, floor_db=1.0): """Локальные максимумы red(f) (нормир. к 0..max) выше min_red_db + min spacing.""" if red.size == 0: return [] r = red.copy() # локальные максимумы d = np.diff(np.sign(np.diff(r))) idx = np.where(d < 0)[0] + 1 out = [] for i in idx: f = freqs[i] if r[i] < min_red_db: continue # spacing if any(abs(f - o) < spacing_hz for o in out): continue out.append(f) if len(out) >= max_peaks: break return out def lp_smooth(x, alpha): y = np.empty_like(x) s = 0.0 for i in range(x.size): s += alpha * (x[i] - s) y[i] = s return y def simulate(x, sr=SR, bands=None, depth=0.864, sharp=10.0, sel=10.0, attack=0.0, release=0.0, mix=100.0, nfft=2048, hop=512, thr_db=6.0, level_ref=1.0): """bands: [{'on':1,'freq':500,'q':1,'sens':12}, ...] дефолт = factory trk.""" if bands is None: bands = [dict(on=1, freq=500.0, q=1.0, sens=12.0)] x = np.asarray(x, dtype=np.float64) fr, w, nf, H = stft_frames(x, nfft, hop) W = eq_weight_db(fr, bands) # loudness per bin: |X| на кадре -> dBFS относительно level_ref out = np.zeros(len(x) + nfft) win_cnt = np.zeros(len(x) + nfft) # параметры from sim_v5 import q_notch_from_sharp as qn_map qn = qn_map(sharp) # Q растёт с sharpness, сужается при глубоком нотче (g->1): Q_eff = qn*(1-g)+1 # (умеренный драйв Q=12 при sharp10; глубокий нотч сильного тона Q~4.5) spacing = max(80.0, nfft / sr * 2.0) # ~ bin-разрешение*2 по умолчанию # selectivity -> spacing & порог (sel больше = меньше ложных пиков) spacing = 6.02 * sel * (nfft / sr) if False else max(80.0, 90.0 * sel) minred = thr_db * (1.0) + (10 - sel) * 0.0 # level_db на кадре = 20log10( |X| / (level_ref * nfft/2) ) - 3dB (peak-bin -> rms) norm = level_ref * nfft / 4.0 # per-frame smoothing alpha (кадр = hop сэмплов) ta = max(0.020 * np.exp(attack / 1.955), 1e-9) tr = max(float(np.interp(release, [0, 1, 2, 5, 10, 100], [0.027, 0.08, 0.11, 0.14, 15, 15])), 1e-9) a_alpha = 1 - np.exp(-hop / (ta * sr)) r_alpha = 1 - np.exp(-hop / (tr * sr)) # state per peak: current gain; мы будем сглаживать маску целиком (nfft/2+1) mask_prev = np.ones(fr.size) mask_cur = np.ones(fr.size) for fi in range(nf): s0 = fi * hop seg = x[s0:s0 + nfft] seg = seg * w X = np.fft.rfft(seg) mag = np.abs(X) # per-bin уровень: |X|/(level_ref*nfft/4) -> dBFS (peak-bin ~ амплитуда тона) lvl_bin = 20 * np.log10(mag / norm + 1e-12) # кадровый уровень для sat: rms кадра lvl = np.sqrt(np.mean(x[s0:s0 + nfft] ** 2) + 1e-12) level_db = 20 * np.log10(lvl / level_ref + 1e-12) # red(f): база от per-bin уровня + EQ-вес с sat(кадровый уровень) red = base_red(lvl_bin, depth) + W * sat_level(level_db) red = np.clip(red, 0, 60) # пики детектим по lvl_bin (тон выделяется над шумовым дном на любом уровне) thr = float(np.max(lvl_bin)) - 40.0 peaks = detect_peaks(lvl_bin, fr, spacing, thr) if not peaks: im = int(np.argmax(lvl_bin)) if lvl_bin[im] > thr and red[im] > 1.0: peaks = [float(fr[im])] factor = np.ones(fr.size) for f_peak in peaks: rr = float(np.interp(f_peak, fr, red)) g = 1 - 10 ** (-rr / 20.0) qeff = max(qn * (1 - g) + 1.0, 0.5) bell = notch_bell(fr, f_peak, qeff) factor = factor * (1 - g * bell) factor = np.clip(factor, 1e-6, 1.0) # env smoothing (целевая маска -> сглаженная) alpha = a_alpha if (factor < mask_prev).mean() > 0.5 else r_alpha mask = mask_prev + alpha * (factor - mask_prev) mask_prev = mask Y = X * mask y = np.fft.irfft(Y, nfft) out[s0:s0 + nfft] += y * w win_cnt[s0:s0 + nfft] += w * w out = out[:len(x)] / np.maximum(win_cnt[:len(x)], 1e-9) if mix < 100: out = (100 - mix) / 100.0 * x + mix / 100.0 * out return out # ---- wav IO (как в sim.py) ---- def read_wav(path): w = wave.open(path, 'rb') sw, nc, n = w.getsampwidth(), w.getnchannels(), w.getnframes() d = np.frombuffer(w.readframes(n), dtype=np.uint8).reshape(n, nc, sw) w.close() ch = d[:, 0, :] if sw == 3: v = (ch[:, 0].astype(np.int64) | (ch[:, 1].astype(np.int64) << 8) | (ch[:, 2].astype(np.int64) << 16)) v = (v ^ (1 << 23)) - (1 << 23) return v.astype(np.float64) / (1 << 23) elif sw == 2: v = (ch[:, 0].astype(np.int64) | (ch[:, 1].astype(np.int64) << 8)) v = (v ^ (1 << 15)) - (1 << 15) return v.astype(np.float64) / (1 << 15) v = ch[:, 0].astype(np.float64) return (v - 128) / 128.0 if __name__ == '__main__': import sys import os # test: tone1k через factory trk (band1@500 sens12) base = '/home/m/soothe-bt' x = read_wav(os.path.join(base, 'tone1k.wav')) y = simulate(x, bands=[dict(on=1, freq=500.0, q=1.0, sens=12.0)]) # ред на 1k sr = SR i0, i1 = int(1.5 * sr), int(2.0 * sr) red = -20 * np.log10(np.sqrt(np.mean(y[i0:i1] ** 2)) / np.sqrt(np.mean(x[i0:i1] ** 2))) print("sim_v5 tone1k factory(b1@500,s12): red = %.2f dB (real trk_b1_500 = 11.91)" % red)