快速卷积算法—混响(二)

本文涉及的完整代码已经上传到github Fast-Convolution-Reverb

卷积混响

针对混响的采样通常需要在房间内架设一对立体声传声器，然后录下该房间对于一个很短的冲击信号的房间响应(也就是混响)，比如说一声枪响。由于很难产生完美的冲击信号，因此可以使用另一种方法来替代，这就是用扬声器来播放一个正弦波的扫频信号。在被记录下来的声音中，原始的信号可能会被去除，只留下混响。之后，用一个卷积混响器(Convolution reverb)加载录制得到的冲击响应信号(IR),并分析生成一个巨大的数学矩阵,以便随后进行卷积运算使用。由于冲击响应中的所有声音都通过了卷积混响器，因此最终得到的混响效果与原始声场所形成的混响非常相似。

摘自《混音指南》 Roey Izhaki

一、房间脉冲响应（Room Impulse Response, RIR）

房间脉冲响应可以简单理解为在一个房间里播放一个时间极短，声音极大的脉冲后录制到的回响，可以记录房间的声学特性，包括直达声、反射、混响等。常见获取房间脉冲响应的办法可以分为实测法和仿真模拟。

实测法：

1.脉冲信号法：

短时间内播放脉冲信号，如气球爆炸声、发令枪声，用麦克风记录响应。这种方法很方便快捷，但是这个方法的缺点在于很难产生完美不失真的脉冲信号，容易夹杂多余的高频或者低频信息。

2. 正弦扫频法

指数扫频信号是一种随时间指数升高的正弦波，时域数学表达式为：
$$s(t)=\sin \left(2\pi \cdot f_0\cdot \frac{e^{t\cdot \ln (f_1/f_0)/T}-1}{\ln (f_1/f_0)}\right)$$
参数说明：

• ( f0 )：起始频率（Hz），通常为 20 Hz。
• ( f1 )：终止频率（Hz），通常为 20 kHz。
• ( T )：扫频总时长（秒）。
• ( t )：时间变量，范围 ( t \in [0, T] )。

操作流程：
1.将麦克风置于待测位置，指向声源（或使用全向麦克风）。
2.播放扫频信号，同时同步录制房间响应信号
3.多次测量取平均，降低随机噪声影响。
4.反卷积提取脉冲响应：
录制信号 ( y(t) ) 与原始扫频信号 ( s(t) ) 的反卷积运算：
h ( t ) = F − 1 { F { y ( t ) } F { s ( t ) } } h(t) = \mathcal{F}^{-1}\left\{ \frac{\mathcal{F}\{y(t)\}}{\mathcal{F}\{s(t)\}} \right\} h(t)=F−1{F{s(t)}F{y(t)}​}
其中 F 表示傅里叶变换，( h(t) ) 为提取的脉冲响应。

• REW软件中的measure模式可以测量，用的线性扫频信号和对数扫频信号

仿真模拟法：

1 . 镜像源法

计算声源在墙壁上的镜像反射路径，叠加直达声和反射声，计算高效适合高频模拟，但是缺失衍射和波动效应。详细可参考matlab的Room Impulse Response Simulation with the Image-Source Method and HRTF Interpolation 算法。

2.射线追踪法

模拟声波以涉嫌形式传播，统计反射路径。不过计算量较大。

二、脉冲响应格式

常见的脉冲响应格式有单通道、立体声、多声道环绕声、Ambisonics B-Format、高阶 Ambisonics（HOA）。

脉冲响应可以在网络上找到很多开源的脉冲响应，比如：OpenAir、 room-impulse-responses

我这里下载到的脉冲响应是4通道的B-Format的脉冲响应，B-Format 是 Ambisonics 技术的核心音频格式，通过数学正交分解将声场表示为球谐函数的系数（W, X, Y, Z），能够完整记录三维空间中的声压和粒子速度信息。

通道组成:
W 通道（全向分量）：记录声压的全向信息，类似单声道信号。
X 通道（前-后轴分量）：记录声场在前后方向的空间梯度（偶极子指向）。
Y 通道（左-右轴分量）：记录声场在左右方向的空间梯度。
Z 通道（上-下轴分量，可选）：记录垂直方向梯度，用于三维声场（部分B-Format省略Z，仅支持水平环绕）。

Ambisonics B-Format 到立体声的解码：

由于大部分音乐是双声道，目前的使用场景也是耳机听歌。所以我们需要把B-format格式的脉冲响应转成左右声道。
，我们假设音源在正前方正负30°的位置，所以左右耳的声音是：
$$\begin{array}{rl}L& =\frac{\sqrt{2}}{2}\left(W+X\cos {\theta }_L+Y\sin {\theta }_L\right)\\ R& =\frac{\sqrt{2}}{2}\left(W+X\cos {\theta }_R+Y\sin {\theta }_R\right)\end{array}$$
其中 ${\theta }_L$ 和 ${\theta }_R$ 分别为左右声道的虚拟扬声器方位角（通常 ${\theta }_L=30^{\circ }$, ${\theta }_R=-30^{\circ }$）。

三、代码

完整代码已经上传到github Fast-Convolution-Reverb

1.预处理脉冲响应

 // 拆分 IR 数据，使用所有B格式通道进行更复杂的空间处理
    std::vector<float> irW(irFrameCount);
    std::vector<float> irX(irFrameCount);
    std::vector<float> irY(irFrameCount);
    std::vector<float> irZ(irFrameCount);
    
    for (unsigned i = 0; i < irFrameCount; i++) {
        // B格式通道: W(全向), X(前/后), Y(左/右), Z(上/下)
        irW[i] = irInterleaved[i * irInfo.channels + 0];
        irX[i] = irInterleaved[i * irInfo.channels + 1];
        irY[i] = irInterleaved[i * irInfo.channels + 2];
        irZ[i] = irInterleaved[i * irInfo.channels + 3];
    }
    
    // 创建左右耳的IR，使用B格式的线性组合来模拟双耳效果
    // 使用基于扬声器方向角的解码矩阵
    std::vector<float> irLeft(irFrameCount);
    std::vector<float> irRight(irFrameCount);
    
    // 定义扬声器角度（以弧度为单位）
    const float theta_L = 30.0f * 3.1415926 / 180.0f;  // 左声道 +30度
    const float theta_R = -30.0f * 3.1415926 / 180.0f; // 右声道 -30度
    
    // 计算解码系数
    const float cos_theta_L = std::cos(theta_L);
    const float sin_theta_L = std::sin(theta_L);
    const float cos_theta_R = std::cos(theta_R);
    const float sin_theta_R = std::sin(theta_R);
    
    // 可选的增益调整系数
    const float wGain = 0.7071f; // W通道增益，通常为1/sqrt(2)
    const float xyGain = 0.6f;   // X和Y通道增益
    const float zGain = 0.0f;    // Z通道增益（通常在水平面布局中不使用）
    
    for (unsigned i = 0; i < irFrameCount; i++) {
        // 使用解码矩阵计算左右声道
        // L = W + X*cos(θL) + Y*sin(θL)
        irLeft[i] = wGain * irW[i] + 
                    xyGain * irX[i] * cos_theta_L + 
                    xyGain * irY[i] * sin_theta_L + 
                    zGain * irZ[i];
        
        // R = W + X*cos(θR) + Y*sin(θR)
        irRight[i] = wGain * irW[i] + 
                     xyGain * irX[i] * cos_theta_R + 
                     xyGain * irY[i] * sin_theta_R + 
                     zGain * irZ[i];
    }

2.卷积初始化

    // 初始化卷积滤波器（每个通道一个 HConvSingle 实例）
    HConvSingle filters[2];
    hcInitSingle(&filters[0], irLeft.data(), irFrameCount, FRAME_SIZE, 1);
    hcInitSingle(&filters[1], irRight.data(), irFrameCount, FRAME_SIZE, 1);

3.卷积处理音频

    // 按帧处理数据，每个帧包含 FRAME_SIZE 个样本
    std::cout << "开始处理卷积..." << std::endl;
    unsigned totalFrames = inputFrameCount;
    for (unsigned frameIdx = 0; frameIdx < totalFrames; frameIdx += FRAME_SIZE) {
        // 计算当前帧的实际处理样本数（最后一帧可能不足 FRAME_SIZE）
        unsigned curFrameSize = FRAME_SIZE;
        if (frameIdx + FRAME_SIZE > totalFrames) {
            curFrameSize = totalFrames - frameIdx;
        }

        // 从交错格式中提取左右声道数据
        for (unsigned i = 0; i < curFrameSize; i++) {
            inputBufferLeft[i] = inputData[(frameIdx + i) * 2];
            inputBufferRight[i] = inputData[(frameIdx + i) * 2 + 1];
        }
        
        // 如果是最后一帧且不足FRAME_SIZE，将剩余部分填充为0
        if (curFrameSize < FRAME_SIZE) {
            for (unsigned i = curFrameSize; i < FRAME_SIZE; i++) {
                inputBufferLeft[i] = 0.0f;
                inputBufferRight[i] = 0.0f;
            }
        }

        // 处理左声道
        hcPutSingle(&filters[0], inputBufferLeft.data());
        hcProcessSingle(&filters[0]);
        hcGetSingle(&filters[0], outputBufferLeft.data());

        // 处理右声道
        hcPutSingle(&filters[1], inputBufferRight.data());
        hcProcessSingle(&filters[1]);
        hcGetSingle(&filters[1], outputBufferRight.data());

4.调整干湿比

        // 将处理后的左右声道数据合并回交错格式，并应用干湿比
        for (unsigned i = 0; i < curFrameSize; i++) {
            // 应用干湿比混合
            float dryLeftSample = inputBufferLeft[i];
            float dryRightSample = inputBufferRight[i];
            float wetLeftSample = outputBufferLeft[i];
            float wetRightSample = outputBufferRight[i];
            
            // 混合干湿信号
            outputData[(frameIdx + i) * 2] = (1.0f - wetDryRatio) * dryLeftSample + wetDryRatio * wetLeftSample;
            outputData[(frameIdx + i) * 2 + 1] = (1.0f - wetDryRatio) * dryRightSample + wetDryRatio * wetRightSample;
        }

总结

本文大致讲述了我在做卷积混响时的处理方法。文中如果有错误欢迎指正。
文中提到了Ambisonics B-Format、HRTF，这些方法也是目前空间音频的主要技术方案，后续的文章我也会围绕空间音频展开技术讨论。