Archive DLL-free personalized HRTF binaural renderer

2026-09-05 02:43:54 +08:00
parent 329445ed25
commit 22a16ab60d
25 changed files with 4285 additions and 109 deletions
@@ -0,0 +1,536 @@
+# JustOneCacophony — 双耳渲染数学
+
+[English](binaural.en.md) · [返回 README](../README.md)
+
+本文说明 `pcm16 + ID11/OAMD → stereo` 路径中的计算、状态和时间对齐。双耳渲染直接接在对象重建之后，不经过 ADM BWF 或 AXML。
+
+## 1. 总体路径与记号
+
+```text
+pcm16：LFE + 15 路对象 PCM
+  → 64-band QMF analysis
+  → 77-band hybrid analysis
+  → 逐对象方向、距离、双耳传递函数和 room send
+  → 对象累加 + room network
+  → hybrid synthesis
+  → QMF synthesis
+  → 961-sample 延迟补偿
+  → stereo WAV
+```
+
+主要记号：
+
+| 符号 | 含义 |
+|---|---|
+| $s=0\ldots15$ | 输入源；$s=0$ 为 LFE，$s=1\ldots15$ 为对象 |
+| $e\in\{L,R\}$ | 左右输出耳 |
+| $p=0\ldots63$ | QMF phase / 时域 hop 内采样 |
+| $k=0\ldots63$ | QMF 子带 |
+| $h=0\ldots76$ | hybrid 子带 |
+| $j=0\ldots35$ | 方向 basis 项 |
+| $m$ | 64-sample QMF 时槽 |
+| $n$ | 时域采样位置 |
+
+每个 QMF hop 为 64 samples，每个双耳控制块为
+
+$$
+N_b=512=8\times64,
+$$
+
+每个 E-AC-3/JOC 音频帧为
+
+$$
+N_f=1536=3N_b.
+$$
+
+## 2. 输入与控制块
+
+输入矩阵为
+
+$$
+x_s[n],\qquad s=0\ldots15.
+$$
+
+`pcm16` 的通道约定为：
+
+```text
+ch0      special LFE
+ch1..15  JOC 对象 1..15
+```
+
+渲染器按连续采样流推进。QMF、hybrid、room 和 synthesis 状态不会在 1536-sample 帧边界清零。
+
+## 3. 64-band QMF analysis
+
+令 $a_{p,\ell}$ 为固定的 64×10 polyphase 系数，$r_{s,\ell,p}[m]$ 为当前和前 9 个 hop 的 phase 历史。奇偶 lag 分别累加：
+
+$$
+E_{s,p}[m]
+=\sum_{\substack{\ell=0\\\ell\text{ even}}}^{9}
+a_{p,\ell}r_{s,\ell,p}[m],
+$$
+
+$$
+O_{s,p}[m]
+=\sum_{\substack{\ell=0\\\ell\text{ odd}}}^{9}
+a_{p,\ell}r_{s,\ell,p}[m].
+$$
+
+对任一 64-vector $v[p]$，定义调制变换
+
+$$
+\mathcal Q(v)_k
+=
+\operatorname{FFT}_{128}
+\left(
+\left[v[p]e^{-j\pi p/128}\right]_{p=0}^{63},
+0_{64}
+\right)_k
+
+e^{-j3\pi(k+1/2)/128}.
+$$
+
+analysis 输出为
+
+$$
+X_{s,k}[m]
+=
+\mathcal Q(O_s)_k
+j(-1)^k\mathcal Q(E_s)_k.
+$$
+
+所有历史、乘加和 FFT 结果使用 `float64/complex128`。
+
+## 4. 77-band hybrid analysis
+
+低 3 个 QMF 子带使用 13-slot FIR 拆分为 16 个 hybrid bands。把复数的实部和虚部分量记为 $i,o\in\{0,1\}$，固定核为 $K_{p,i,\ell,h,o}$：
+
+$$
+H_{s,h,o}[m]
+=
+\sum_{p=0}^{2}
+\sum_{i=0}^{1}
+\sum_{\ell=0}^{12}
+X_{s,p,i}[m-\ell]K_{p,i,\ell,h,o},
+\qquad h=0\ldots15.
+$$
+
+其余 61 个 hybrid bands 是 QMF 3..63 的 6-slot 延迟：
+
+$$
+H_{s,16+q}[m]=X_{s,3+q}[m-6],
+\qquad q=0\ldots60.
+$$
+
+因此 hybrid vector 的顺序为：
+
+```text
+0..15   低 3 个 QMF bands 的细分
+16..76  延迟后的 QMF bands 3..63
+```
+
+## 5. OAMD 坐标与时间轴
+
+### 5.1 Q15 坐标到 Cartesian
+
+对象状态中的 $q_1,q_2,q_3$ 先恢复到离散位置网格：
+
+$$
+u_1=\min\left(1,\frac{\operatorname{round}(62q_1/32767)}{62}\right),
+$$
+
+$$
+u_2=\min\left(1,\frac{\operatorname{round}(62q_2/32767)}{62}\right),
+$$
+
+$$
+u_3=\operatorname{clip}\left(
+\frac{\operatorname{round}(15q_3/32767)}{15},-1,1\right).
+$$
+
+ADM Cartesian 坐标为
+
+$$
+(X,Y,Z)=(2u_1-1,\ 1-2u_2,\ u_3).
+$$
+
+### 5.2 更新时间
+
+一条位置更新的编码时刻为
+
+$$
+n_{\mathrm{coded}}
+=n_{\mathrm{frame}}
+n_{\mathrm{outer}}
+n_{\mathrm{block}}.
+$$
+
+首个有效状态作为 sample 0 的初始位置。后续更新加入对象 PCM 延迟 $D_o$，默认
+
+$$
+D_o=1473.
+$$
+
+若 ramp duration 为 $R>64$，连续运动为
+
+$$
+n_{\mathrm{start}}=n_{\mathrm{coded}}+D_o+64,
+$$
+
+$$
+R_{\mathrm{eff}}=R-64,
+$$
+
+$$
+\mathbf p[n]
+=(1-\alpha)\mathbf p_0+\alpha\mathbf p_1,
+\qquad
+\alpha=\frac{n-n_{\mathrm{start}}}{R_{\mathrm{eff}}}.
+$$
+
+当 $R\le64$ 时，目标位置在 $n_{\mathrm{coded}}+D_o$ 直接生效。
+
+Rosella 参数在每个 512-sample block 起点求值，并用于该块的 8 个 hybrid slots。
+
+## 6. 距离 profile 与方向
+
+普通对象只使用 Near、Mid、Far 三个 profile。每个 profile 包含：
+
+- 三轴负/正边界 $b_{x-},b_{x+},b_{y-},b_{y+},b_{z-},b_{z+}$；
+- 距离尺度 $D$ 和倒数尺度 $D^{-1}$；
+- 三轴内部尺度 $a_x,a_y,a_z$；
+- 最小归一化半径 $\rho_{\min}$。
+
+Cartesian 坐标经过 Q15 metadata grid 后换成内部前、侧、上轴，乘以 profile 尺度：
+
+$$
+\mathbf s=(a_zq_f,\ a_xq_l,\ a_yq_v).
+$$
+
+若射线超出 profile 边界，则用单一比例 $\lambda\le1$ 缩放：
+
+$$
+\mathbf s' = \lambda\mathbf s.
+$$
+
+随后
+
+$$
+\rho=\|\mathbf s'\|_2,
+\qquad
+\rho_c=\max(\rho,\rho_{\min}),
+\qquad
+\alpha=\frac{\rho}{\rho_c},
+$$
+
+$$
+\mathbf d=
+\begin{cases}
+\mathbf s'/\rho,&\rho>0,\\
+(1,0,0),&\rho=0,
+\end{cases}
+$$
+
+物理半径为
+
+$$
+R=D\rho.
+$$
+
+## 7. 36 项方向 basis
+
+方向 $\mathbf d=(x,y,z)$ 被展开为 36 项实值多项式：
+
+$$
+\mathbf b(\mathbf d)=
+[1,x,y,z,x^2-\tfrac13,xy,xz,y^2-\tfrac13,yz,\ldots]^T.
+$$
+
+完整顺序由 `rosella_model.direction_basis()` 固定。最高次数为 5；所有 field 系数必须按该顺序点积，不能交换 basis 项。
+
+逐耳 basis 会根据耳偏移重新归一化。令耳偏移标量为 $e$：
+
+$$
+\epsilon=\frac{eD^{-1}}{\rho_c},
+$$
+
+$$
+\mathbf d_-=
+\frac{(x,y-\epsilon,z)}{\|(x,y-\epsilon,z)\|_2},
+\qquad
+\mathbf d_+=
+\frac{(x,y+\epsilon,z)}{\|(x,y+\epsilon,z)\|_2}.
+$$
+
+## 8. 逐耳路径与 ITD
+
+归一化路径长度为
+
+$$
+\ell_-=\rho_c\sqrt{x^2+(y-\epsilon)^2+z^2},
+$$
+
+$$
+\ell_+=\rho_c\sqrt{x^2+(y+\epsilon)^2+z^2}.
+$$
+
+模型允许通过 36-vector 对路径加入非负方向修正：
+
+$$
+\ell'_e
+=
+\ell_e
+
+\max(\mathbf v_e^T\mathbf b_e,0)\,2cD^{-1}.
+$$
+
+耳间延迟为
+
+$$
+\tau
+=|\ell'_+-\ell'_-|\,D\frac{f_s}{343.3}\alpha,
+\qquad f_s=48000.
+$$
+
+路径较长的一耳应用 hybrid-band 相位：
+
+$$
+P_h=e^{j\omega_h\tau},
+$$
+
+其中 $\omega_h$ 由模型的 20 个 hybrid group 参数递推到 77 个 bands。
+
+## 9. 方向 field 与直达增益
+
+左右耳各有一个 77×36 complex field：
+
+$$
+C_{e,h}(\mathbf d_e)
+=
+\sum_{j=0}^{35}F_{e,h,j}b_j(\mathbf d_e).
+$$
+
+另一次耳路径计算给出左右权重：
+
+$$
+w_L=\frac{\ell_+}{\sqrt{\ell_-^2+\ell_+^2}},
+\qquad
+w_R=\frac{\ell_-}{\sqrt{\ell_-^2+\ell_+^2}}.
+$$
+
+有效距离为
+
+$$
+R_e=\rho\,s_dD,
+$$
+
+其中 $s_d$ 为模型距离标量。Mid/Far 的公共衰减和 room send 为
+
+$$
+g_c=\frac{1}{\sqrt{1+s_rR_e^2}},
+$$
+
+$$
+g_{\mathrm{room}}=R_eg_c.
+$$
+
+Near 使用
+
+$$
+g_c=1,
+\qquad
+g_{\mathrm{room}}=0.
+$$
+
+令 $C_{e,h}^{(0)}$ 为 field 的第 0 个 basis 系数，中心保护项为 $1-\alpha$。普通直达传递函数可写成
+
+$$
+G_{L,h}
+=g_c\left[C_{L,h}w_L\alpha+C_{L,h}^{(0)}c_L(1-\alpha)\right],
+$$
+
+$$
+G_{R,h}
+=g_c\left[C_{R,h}w_R\alpha+C_{R,h}^{(0)}c_R(1-\alpha)\right].
+$$
+
+$c_L,c_R$ 由模型的耳权重配置选择；路径较长的一耳再乘 $P_h$。
+
+## 10. LFE 传递函数
+
+LFE 不进入普通对象方向计算。其传递函数为
+
+$$
+G_{L,h}^{\mathrm{LFE}}=G_{R,h}^{\mathrm{LFE}}=
+\begin{cases}
+g_h,&0\le h<16,\\
+0,&16\le h<77.
+\end{cases}
+$$
+
+前 16 个固定系数为
+
+```text
+ 2.60290003,  1.80741799,  0.659342408, -0.0275855921,
+-0.105803289, -0.0699509233, 0.0749056414, -0.00919809937,
+ 0.00349014648,-0.0158600751,-0.000723021978,0.00188189559,
+-0.000421735429,0.0000329252762,0.0000317397971,0.000000580376991
+```
+
+LFE 的 room send 恒为 0。
+
+## 11. 对象累加与 room input
+
+每个 hybrid slot 的直达输出为
+
+$$
+Y^{\mathrm{direct}}_{e,h}
+=
+\sum_{s=0}^{15}H_{s,h}G_{s,e,h}.
+$$
+
+room 输入为
+
+$$
+U_h
+=
+\sum_{s=1}^{15}H_{s,h}g_{\mathrm{room},s}.
+$$
+
+LFE 不进入该和式。
+
+## 12. Room network
+
+room 只处理前 64 个 hybrid bands。输入先乘
+
+$$
+g_0=0.70710677.
+$$
+
+对每级 all-pass，设延迟样本为 $d[n]$、系数为 $a$：
+
+$$
+r[n]=x[n]-ad[n],
+$$
+
+$$
+y[n]=ar[n]+d[n].
+$$
+
+all-pass 输出复制到 4 个 FDN branches。设延迟输出为 $\mathbf d_h[m]$、4×4 混合矩阵为 $M$：
+
+$$
+\mathbf b_h[m]
+=U_h[m]\mathbf 1+M\mathbf d_h[m].
+$$
+
+每个 branch 使用复反馈系数 $f_{h,i}$：
+
+$$
+m_{h,i}[m]=f_{h,i}b_{h,i}[m].
+$$
+
+主 tap、可选额外 tap 和左右输出矩阵合成为
+
+$$
+Y^{\mathrm{room}}_{e,h}[m]
+=
+\sum_{i=0}^{3}O_{e,h,i}z_{h,i}[m].
+$$
+
+最终 hybrid 输出为
+
+$$
+Y_{e,h}=Y^{\mathrm{direct}}_{e,h}+Y^{\mathrm{room}}_{e,h}.
+$$
+
+Python 后端把该递归网络展开为有限 complex FIR 并使用 overlap-add；C++ 后端直接保持递归状态。两者均跨帧连续。
+
+## 13. Hybrid synthesis
+
+hybrid synthesis 是 154 项稀疏即时映射。令映射项为 $(h,i,k,o,w)$，其中 $i,o$ 表示实部或虚部，则
+
+$$
+Q_{e,k,o}[m]
+\mathrel{+}=
+Y_{e,h,i}[m]w.
+$$
+
+输出是每耳 64 个 complex QMF bands。
+
+## 14. QMF synthesis
+
+每耳 QMF vector 先展开为 128 项实向量
+
+$$
+\mathbf q_e=[\Re Q_{e,0},\Im Q_{e,0},\ldots,\Re Q_{e,63},\Im Q_{e,63}]^T.
+$$
+
+对每个 phase $p$ 和 rank $r=0\ldots3$：
+
+$$
+f_{e,p,r}[m]
+=\mathbf b_{p,r}^T\mathbf q_e[m].
+$$
+
+使用 10-slot taps 合成时域样本：
+
+$$
+y_e[64m+p]
+=
+\sum_{\ell=0}^{9}
+\sum_{r=0}^{3}
+t_{p,\ell,r}f_{e,p,r}[m-\ell].
+$$
+
+## 15. 延迟、尾声和输出
+
+完整 filterbank 的固定延迟为
+
+$$
+L=961\ \text{samples}.
+$$
+
+只在连续流起点丢弃一次前 $L$ 个输出 samples。输入结束后继续送零，以释放 QMF、hybrid 和 room 状态。尾声裁切只作用于文件末端：
+
+$$
+\max(|y_L[n]|,|y_R[n]|) > 10^{-8}
+$$
+
+的最后一个 sample 被保留，同时输出长度不得短于源 PCM 长度。
+
+所有内部状态、参数计算和对象累加使用 `float64/complex128`。最终 writer 才转换为 float32 或 PCM24。
+
+## 16. Python 与 C++ 后端
+
+两套后端共享：
+
+- 同一份 QMF/hybrid 固定表；
+- 同一份模型解析结果；
+- 同一套 512-sample 参数更新时间轴；
+- 同一组逐对象 complex gains 和 room sends；
+- 同一 961-sample 延迟补偿与尾声策略。
+
+C++ 后端以 512-sample block 为处理单位，内部持有 QMF、hybrid、room 和 synthesis 状态。Python 只负责模型解析、OAMD 时间轴和每块参数更新。
+
+## 17. 模型文件
+
+默认路径为
+
+```text
+HRTF/binaural.personalized_headphone
+```
+
+也可通过
+
+```text
+--personalized-headphone PATH
+```
+
+指定其它文件。
+
+`.personalized_headphone` 中的 int32/Q15 参数在解析后提升为 float64。任意 SOFA FIR 不能只通过数组重排变成该参数模型；若要转换，需要拟合方向 fields、ITD、距离 profile、耳几何和 room 参数。
+
+## 18. 适用范围
+
+当前路径处理 15 个普通点对象和 1 路 special LFE。对象 extent、spread、diffuse、divergence、channel lock，以及未实现的 OAMD element 变体不在本公式范围内。