扫一扫
均衡器系列
参数均衡器(PEQ)制作指南（一）
参数化均衡器(PEQ)制作指南（二）
文章目录
前言
经常调音的同学肯定知道autoeq这个功能，autoeq的功能就是通过算法自动生成 EQ 参数，目的是使音频设备的频响接近目标曲线（比如哈曼曲线），避免手动调节复杂且不精确的问题。
下文代码可以在github上获得：C++ AutoEQ 实现可以直接下载调试
一、autoeq是如何实现的
autoeq的整体步骤大概分为以下几点：
- 测量与获得频响数据：得到需要修改的原始频响和目标频响。
- 根据频响偏差在指定频段内寻找最大偏差点，构建滤波器的初始值(fc、gain、Q)
- 使用梯度下降法微调滤波器参数，使整体 EQ 响应更接近目标曲线
- 输出 EQ 配置
以上5点就是我总结的autoeq的大概实现流程，具体实现方法我们下面讲。
二、数据的前期处理
这里不具体讲怎么测量设备的频响曲线了，但是需要注意频响曲线需要经过处理后才能传入算法开始计算。
2.1频率轴对齐
需要针对两条频响曲线的频率轴做对齐，或者说重采样到指定的频率轴上。例如，原始频响的频率轴是20，20.14，20.29，20.43；目标频响的频率轴是：20.64，20.79，20.93，21.08。这时需要把原始频响的幅度数据对应到目标频响的频率轴上。例如下面的resampleTo函数。
/**
* 线性插值函数
* 根据给定的目标频率，通过线性插值计算对应的幅度值
* @param target_freq 目标频率点（Hz）
* @return 插值得到的幅度值（dB）
*/
double interpolate(double target_freq) const {
if (frequencies.empty()) {
return 0.0; // 空数据返回0
}
// 如果目标频率超出数据范围，返回边界值
if (target_freq <= frequencies.front()) {
return magnitudes.front(); // 返回最低频率的幅度
}
if (target_freq >= frequencies.back()) {
return magnitudes.back(); // 返回最高频率的幅度
}
// 在数据范围内查找插值区间
for (size_t i = 0; i < frequencies.size() - 1; ++i) {
if (target_freq >= frequencies[i] && target_freq <= frequencies[i + 1]) {
// 执行线性插值计算
// 公式：y = y1 + (x - x1) * (y2 - y1) / (x2 - x1)
double ratio = (target_freq - frequencies[i]) / (frequencies[i + 1] - frequencies[i]);
return magnitudes[i] + ratio * (magnitudes[i + 1] - magnitudes[i]);
}
}
return 0.0; // 理论上不应该到达这里
}
/**
* 重新采样到指定的频率点
* 将当前频率响应数据重新采样到新的频率网格上
* @param target_frequencies 目标频率点数组
*/
void resampleTo(const std::vector<double>& target_frequencies) {
std::vector<double> new_magnitudes;
new_magnitudes.reserve(target_frequencies.size()); // 预分配内存提高性能
// 对每个目标频率点进行插值计算
for (double freq : target_frequencies) {
new_magnitudes.push_back(interpolate(freq));
}
// 更新当前对象的数据
frequencies = target_frequencies;
magnitudes = new_magnitudes;
}
2.2响度对齐
需要在指定频率点对齐两个频响曲线的幅度基准，例如下图，两个频响曲线在100hz做了响度对齐。
2.3 倍频程平滑
需要对两组数据做平滑，避免极端峰值影响计算结果。我这里用的是1/24倍频程平滑。
三、滤波器的实现和初始化
本方法目前只使用峰值滤波器来实现autoeq的效果，峰值滤波器的具体实现原理及方法可以回顾第一节的内容。这里需要重点讲一下绘制滤波器的频响曲线，和初始化滤波器系数。
3.1绘制滤波器频响曲线
要计算滤波器对频响曲线的影响，首先需要计算出滤波器的频响。根据双二阶传递函数公式：
void updateFrequencyResponse() {
frequency_response.resize(frequencies.size());
for (size_t i = 0; i < frequencies.size(); ++i) {
double f = frequencies[i];
double omega = 2.0 * M_PI * f / fs; // 归一化角频率
// 计算 z = e^(jω)
std::complex<double> z = std::exp(std::complex<double>(0.0, omega));
std::complex<double> z_inv = 1.0 / z; // z^(-1)
std::complex<double> z_inv2 = z_inv * z_inv; // z^(-2)
// 获取双二阶滤波器系数
auto coeffs = getBiquadCoefficients();
// 双二阶传递函数 H(z) = (b0 + b1*z^-1 + b2*z^-2) / (a0 + a1*z^-1 + a2*z^-2)
std::complex<double> numerator = coeffs.b0 + coeffs.b1 * z_inv + coeffs.b2 * z_inv2;
std::complex<double> denominator = coeffs.a0 + coeffs.a1 * z_inv + coeffs.a2 * z_inv2;
std::complex<double> H = numerator / denominator;
// 转换为dB幅度
double magnitude = std::abs(H);
frequency_response[i] = 20.0 * std::log10(std::max(magnitude, 1e-10));
}
}
该方法对所有的双二阶滤波器有效。
3. 2 滤波器参数的初始化
- 用目标曲线减去原始曲线，得到误差曲线。
- 遍历频响曲线误差，找到误差最大值所在的索引
- 把这个索引对应的频率作为初始的 fc (中心频率)
- 把该点的误差作为初始的 gain (增益)
- 根据这个峰值点的 -3dB 宽度来估算 Q 值
- 根据得到的参数生成滤波器频响，并将其从误差曲线中抵消（误差曲线减去滤波器频响）
- 将更新后的误差曲线作为新的输入，进入下一轮迭代。
/**
* 为指定频率范围初始化滤波器
* @param filter 要初始化的滤波器
* @param target_response 目标响应
* @param freq_min 频率范围下限
* @param freq_max 频率范围上限
*/
void initializeFilterForFreqRange(PeakingFilter& filter, const std::vector<double>& target_response,
double freq_min, double freq_max) {
if (target_response.size() != target.frequencies.size()) {
return;
}
// 在指定频率范围内寻找最大偏差
double max_deviation = 0.0;
size_t max_index = 0;
bool found_in_range = false;
for (size_t i = 0; i < target_response.size(); ++i) {
double freq = target.frequencies[i];
if (freq >= freq_min && freq <= freq_max) {
double deviation = std::abs(target_response[i]);
if (deviation > max_deviation) {
max_deviation = deviation;
max_index = i;
found_in_range = true;
}
}
}
if (found_in_range && max_deviation > 0.1) {
double fc = target.frequencies[max_index];
double gain = target_response[max_index];
// 在指定频率范围内估算Q值
double half_gain = gain * 0.7071067811;
size_t left_index = max_index;
size_t right_index = max_index;
// 向左寻找半增益点
while (left_index > 0 && target.frequencies[left_index] >= freq_min &&
std::abs(target_response[left_index]) > std::abs(half_gain)) {
left_index--;
}
// 向右寻找半增益点
while (right_index < target_response.size() - 1 && target.frequencies[right_index] <= freq_max &&
std::abs(target_response[right_index]) > std::abs(half_gain)) {
right_index++;
}
double q = 1.0;
if (right_index > left_index) {
double bandwidth = target.frequencies[right_index] - target.frequencies[left_index];
if (bandwidth > 0) {
q = fc / bandwidth;
q = std::max(filter.getMinQ(), std::min(filter.getMaxQ(), q));
}
}
filter.setParameters(fc, q, gain);
} else {
// 如果在指定范围内没有找到显著偏差，使用范围中心频率
double center_freq = (freq_min + freq_max) / 2.0;
filter.setParameters(center_freq, 1.0, 0.0);
}
}
/**
* 从剩余目标中减去滤波器的影响
* @param filter 滤波器
* @param remaining_target 剩余目标响应
*/
void subtractFilterResponse(const PeakingFilter& filter, std::vector<double>& remaining_target) {
const auto& filter_response = filter.getFrequencyResponse();
for (size_t i = 0; i < remaining_target.size(); ++i) {
remaining_target[i] -= filter_response[i];
}
}
四、优化滤波器参数
优化算法的第一步，也是最关键一步，就是明确定义需要优化的「目标函数」，也称为损失函数。我这里使用的是加权均方根误差（Weighted RMSE, Root Mean Squared Error）。
4.1 确定损失函数
- 计算所有初始化滤波器应用在原始频响曲线后的新的频响曲线，得到当前频响曲线。
- 计算误差：计算目标曲线和当前频响曲线的差。
- 根据听觉感知做加权策略，人耳对中低频(4khz以下)最为敏感。
- 计算感知加权后的 RMSE。
- 最终目标就是不断优化使得Loss的值最小，Loss的值越小说明越接近目标曲线
4.2 计算梯度
梯度可以告诉我们该往哪个方向调参可以让损失下降。
- 梯度为正 → 增大该参数会让损失变大 → 应该减小参数
- 梯度为负 → 增大该参数会让损失减小 → 应该增大参数
由于滤波器参数（fc, Q, gain）对应的损失函数没有显式解析式，我们采用 数值微分（Finite Difference）来近似梯度：
- 对某个参数施加一个很小的扰动 +epsilon，重新计算一次损失 → loss_plus
- 对该参数施加一个相反的扰动 -epsilon，再计算一次损失 → loss_minus
梯度近似为：
θ \theta θ 表示 fc、Q 或 gain 中的某一个
/**
* 计算指定参数的数值梯度
* 使用有限差分方法近似计算损失函数对参数的偏导数
* @param filter 要计算梯度的滤波器
* @param param_type 参数类型（0=中心频率, 1=Q值, 2=增益）
* @param epsilon 微分步长
* @return 该参数的梯度值
*/
double calculateGradient(PeakingFilter& filter, int param_type, double epsilon) {
// 保存原始参数值
double orig_fc = filter.getFc();
double orig_q = filter.getQ();
double orig_gain = filter.getGain();
double loss_plus = 0.0;
double loss_minus = 0.0;
// ----------- 正向扰动 -----------
switch (param_type) {
case 0: // 中心频率
filter.setFc(orig_fc + epsilon * orig_fc);
break;
case 1: // Q值
filter.setQ(orig_q + epsilon);
break;
case 2: // 增益
filter.setGain(orig_gain + epsilon);
break;
}
loss_plus = calculateLoss();
// 恢复原始值
filter.setParameters(orig_fc, orig_q, orig_gain);
// ----------- 反向扰动 -----------
switch (param_type) {
case 0: // 中心频率
filter.setFc(orig_fc - epsilon * orig_fc);
break;
case 1: // Q值
filter.setQ(orig_q - epsilon);
break;
case 2: // 增益
filter.setGain(orig_gain - epsilon);
break;
}
loss_minus = calculateLoss();
// 恢复原始值
filter.setParameters(orig_fc, orig_q, orig_gain);
// ----------- 中心差分计算梯度 -----------
double gradient = 0.0;
switch (param_type) {
case 0: // 中心频率的梯度
gradient = (loss_plus - loss_minus) / (2.0 * epsilon * orig_fc);
break;
case 1: // Q值的梯度
gradient = (loss_plus - loss_minus) / (2.0 * epsilon);
break;
case 2: // 增益的梯度
gradient = (loss_plus - loss_minus) / (2.0 * epsilon);
break;
}
return gradient;
}
我上面代码使用的是中心差分的方式来计算梯度，这样的方式相对于前向差分和后向差分精度更高一些，但是增加了计算。
4.3 参数更新
得到梯度之后，就可以用 梯度下降法（Gradient Descent） 来更新参数：
η
\eta
η表示学习率，决定更新的步长大小
/**
* 更新所有滤波器的参数
* 使用数值微分计算梯度，然后应用梯度下降更新规则
*/
void updateParameters() {
const double epsilon = 1e-5; // 数值微分的小步长
// 对每个滤波器进行参数更新
for (auto& filter : filters) {
// 计算三个参数（中心频率、Q值、增益）的梯度
double grad_fc = calculateGradient(filter, 0, epsilon); // 中心频率梯度
double grad_q = calculateGradient(filter, 1, epsilon); // Q值梯度
double grad_gain = calculateGradient(filter, 2, epsilon); // 增益梯度
// 使用梯度下降更新参数（负梯度方向）
double new_fc = filter.getFc() - learning_rate * grad_fc * filter.getFc() * 0.1;
double new_q = filter.getQ() - learning_rate * grad_q * 0.1;
double new_gain = filter.getGain() - learning_rate * grad_gain;
// 应用新参数到滤波器
filter.setParameters(new_fc, new_q, new_gain);
}
}
五、完成计算
可以看到蓝色的是经过处理后频响曲线，和橘色目标频响曲线十分接近。
总结
梯度下降法只是比较基础的优化算法，可以改进成更复杂的遗传算法、粒子群算法等，提高全局收敛性。
github:C++ AutoEQ 实现
1446
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
