快速卷积（C++）
快速卷积
最近在看EqualizerAPO的源码发现其中使用的卷积算法效率很高，使用的是overlap-add卷积法和非均匀分割卷积算法的结合。
简单研究了一下发现快速卷积的方法有很多，比如FFT卷积、overlap-add卷积、overlap-save卷积、非均匀分割卷积算法和均匀分割卷积算法。
本文结合代码分析EqualizerAPO中libHybridConv.c的处理方法。
卷积
卷积是一种数学运算，在信号处理、图像处理、机器学习等领域都有广泛应用。
y [ n ] = ∑ k = − ∞ ∞ x [ k ] h [ n − k ] y[n] = \sum_{k=-\infty}^{\infty} x[k] h[n-k] y[n]=k=−∞∑∞x[k]h[n−k]
混合卷积(Hybrid Convolution)
- 整体功能:
这个库主要用于实现音频信号的卷积运算,通常用于音频效果处理(如混响、均衡器等)。它提供了三种不同级别的卷积实现:
- Single(单级)卷积
- Dual(双级)卷积
- Triple(三级)卷积
- 主要特点:
- 使用FFT(快速傅里叶变换)来加速卷积运算
- 支持实时处理
- 针对不同长度的输入做了优化处理
- 使用SSE指令集优化(在Windows平台)
- 核心数据结构:
- HConvSingle : 单级卷积的结构体,包含FFT缓冲区、滤波器系数等
- HConvDual : 双级卷积结构体,组合了长短两种不同长度的单级卷积
- HConvTripple : 三级卷积结构体,组合了中短两种长度的卷积
单级卷积方法运用了overlap-add方法，双级和三级卷积方法在单级卷积方法的基础上对输入的脉冲响应做了非均匀的分割。
单级卷积算法
第一步：初始化函数 hcInitSingle
将一个长的脉冲响应(滤波器) h[n] 划分成多个适合快速 FFT 处理的小段，全部转换为频域存储，分阶段准备好后续实时音频流的频域卷积操作。
1.把长脉冲响应 h[n] 分成多个小段，每段 flen 长度，方便后续进行 分块 FFT 计算。
2.计算每段的 FFT 结果，并存储到频域缓冲区，避免重复计算，加速卷积。
3.初始化各种缓冲区（时域、频域、混合缓冲区、历史缓冲区），为 Overlap-Add 做准备。
4.设置 steptask，让长脉冲响应 h 进行分步处理，提高实时性，减少计算负载峰值。
void hcInitSingle(HConvSingle *filter, float *h, int hlen, int flen, int steps)
{
int i, j, size, num, pos;
float gain;
// 初始化处理步骤计数器，表示当前处理到第几步
filter->step = 0;
// 设定每帧音频的最大处理步数
filter->maxstep = steps;
// 初始化混合缓冲区索引，表示当前混合缓冲区的位置
filter->mixpos = 0;
// 每帧音频的采样点数（FFT 计算基数）
filter->framelength = flen;
// 分配时域缓冲区（长度为 2 * flen，因为 FFT 需要扩展长度）
size = sizeof(float) * 2 * flen;
filter->dft_time = (float *)fftwf_malloc(size);
// 分配频域缓冲区（复数数组，flen + 1 点存储实数FFT结果）
size = sizeof(fftwf_complex) * (flen + 1);
filter->dft_freq = (fftwf_complex*)fftwf_malloc(size);
// 分配输入频域缓冲区（分别存储实部和虚部）
size = sizeof(float) * (flen + 1);
filter->in_freq_real = (float*)fftwf_malloc(size);
filter->in_freq_imag = (float*)fftwf_malloc(size);
// 计算滤波器的分段数量（每段 flen 长度）
filter->num_filterbuf = (hlen + flen - 1) / flen;
// 分配每个处理步需要处理的滤波器分段任务
size = sizeof(int) * (steps + 1);
filter->steptask = (int *)malloc(size);
num = filter->num_filterbuf / steps;
for (i = 0; i <= steps; i++)
filter->steptask[i] = i * num;
// 处理步数调整，确保任务均匀分布
if (filter->steptask[1] == 0)
pos = 1;
else
pos = 2;
num = filter->num_filterbuf % steps;
for (j = pos; j < pos + num; j++)
{
for (i = j; i <= steps; i++)
filter->steptask[i]++;
}
// 分配滤波器频域缓冲区（存储各段 FFT 结果）
size = sizeof(float*) * filter->num_filterbuf;
filter->filterbuf_freq_real = (float**)fftwf_malloc(size);
filter->filterbuf_freq_imag = (float**)fftwf_malloc(size);
for (i = 0; i < filter->num_filterbuf; i++)
{
size = sizeof(float) * (flen + 1);
filter->filterbuf_freq_real[i] = (float*)fftwf_malloc(size);
filter->filterbuf_freq_imag[i] = (float*)fftwf_malloc(size);
}
// 计算混合缓冲区数量（滤波器分段数 + 1，用于累积结果）
filter->num_mixbuf = filter->num_filterbuf + 1;
// 分配混合缓冲区（用于累积频域结果）
size = sizeof(float*) * filter->num_mixbuf;
filter->mixbuf_freq_real = (float**)fftwf_malloc(size);
filter->mixbuf_freq_imag = (float**)fftwf_malloc(size);
for (i = 0; i < filter->num_mixbuf; i++)
{
size = sizeof(float) * (flen + 1);
filter->mixbuf_freq_real[i] = (float*)fftwf_malloc(size);
filter->mixbuf_freq_imag[i] = (float*)fftwf_malloc(size);
memset(filter->mixbuf_freq_real[i], 0, size);
memset(filter->mixbuf_freq_imag[i], 0, size);
}
// 分配历史缓冲区（用于 Overlap-Add 技术）
size = sizeof(float) * flen;
filter->history_time = (float *)fftwf_malloc(size);
memset(filter->history_time, 0, size);
// 创建 FFT 和 IFFT 计划（FFTW 预计算优化，提高性能）
filter->fft = fftwf_plan_dft_r2c_1d(2 * flen, filter->dft_time, filter->dft_freq, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT);
filter->ifft = fftwf_plan_dft_c2r_1d(2 * flen, filter->dft_freq, filter->dft_time, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT);
// 预计算滤波器的频域数据
gain = 0.5f / flen;
size = sizeof(float) * 2 * flen;
memset(filter->dft_time, 0, size);
for (i = 0; i < filter->num_filterbuf - 1; i++)
{
for (j = 0; j < flen; j++)
filter->dft_time[j] = gain * h[i * flen + j];
fftwf_execute(filter->fft);
for (j = 0; j < flen + 1; j++)
{
filter->filterbuf_freq_real[i][j] = filter->dft_freq[j][0];
filter->filterbuf_freq_imag[i][j] = filter->dft_freq[j][1];
}
}
// 处理最后一段滤波器（可能不足 flen，需补 0）
for (j = 0; j < hlen - i * flen; j++)
filter->dft_time[j] = gain * h[i * flen + j];
size = sizeof(float) * ((i + 1) * flen - hlen);
memset(&(filter->dft_time[hlen - i * flen]), 0, size);
fftwf_execute(filter->fft);
for (j = 0; j < flen + 1; j++)
{
filter->filterbuf_freq_real[i][j] = filter->dft_freq[j][0];
filter->filterbuf_freq_imag[i][j] = filter->dft_freq[j][1];
}
}
第二步：接收输入音频数据 hcPutSingle
核心作用：
1.接收一帧音频数据 x[n]（flen 个采样点）
2.将 x[n] 拷贝到 dft_time，并补零（变成 2 * flen 点）
3.执行 FFT，将数据转换到频域，存入 in_freq_real 和 in_freq_imag
4.存储为频域数据，后续 hcProcessSingle 进行频域卷积计算
void hcPutSingle(HConvSingle *filter, float *x)
{
int j, flen, size;
// 获取每帧音频数据的长度 flen
flen = filter->framelength;
// 复制输入数据 x 到 FFT 时域缓冲区 dft_time
size = sizeof(float) * flen;
memcpy(filter->dft_time, x, size);
// 在后半部分补零（2 * flen 点 FFT 需要补零）
memset(&(filter->dft_time[flen]), 0, size);
// 执行 FFT，把时域信号变换到频域
fftwf_execute(filter->fft);
// 提取 FFT 结果，存入 in_freq_real 和 in_freq_imag
for (j = 0; j < flen + 1; j++)
{
filter->in_freq_real[j] = filter->dft_freq[j][0]; // 取 FFT 结果的实部
filter->in_freq_imag[j] = filter->dft_freq[j][1]; // 取 FFT 结果的虚部
}
}
第三步 处理音频 hcProcessSingle
1.在频域执行复数乘法（输入信号的 FFT 结果 X(f) 和滤波器的 FFT H(f) 相乘）。
2.累加计算结果到混合缓冲区，用于后续 hcGetSingle 做 IFFT 变回时域。
3.利用 SIMD (SSE 指令) 进行并行优化，提高计算效率。
频域：
Y
(
f
)
=
X
(
f
)
⋅
H
(
f
)
Y(f) = X(f) \cdot H(f)
Y(f)=X(f)⋅H(f)
其中 (X(f)) 和 (H(f)) 都是复数，所以需要进行复数乘法：
Y real = X real ⋅ H real − X imag ⋅ H imag Y_{\text{real}} = X_{\text{real}} \cdot H_{\text{real}} - X_{\text{imag}} \cdot H_{\text{imag}} Yreal=Xreal⋅Hreal−Ximag⋅Himag
Y imag = X real ⋅ H imag + X imag ⋅ H real Y_{\text{imag}} = X_{\text{real}} \cdot H_{\text{imag}} + X_{\text{imag}} \cdot H_{\text{real}} Yimag=Xreal⋅Himag+Ximag⋅Hreal
void hcProcessSingle(HConvSingle *filter)
{
#if 0
int s, n, start, stop, flen;
float *x_real;
float *x_imag;
float *h_real;
float *h_imag;
float *y_real;
float *y_imag;
flen = filter->framelength;
x_real = filter->in_freq_real;
x_imag = filter->in_freq_imag;
start = filter->steptask[filter->step];
stop = filter->steptask[filter->step + 1];
for (s = start; s < stop; s++)
{
n = (s + filter->mixpos) % filter->num_mixbuf;
y_real = filter->mixbuf_freq_real[n];
y_imag = filter->mixbuf_freq_imag[n];
h_real = filter->filterbuf_freq_real[s];
h_imag = filter->filterbuf_freq_imag[s];
for (n = 0; n < flen + 1; n++)
{
y_real[n] += x_real[n] * h_real[n] -
x_imag[n] * h_imag[n];
y_imag[n] += x_real[n] * h_imag[n] +
x_imag[n] * h_real[n];
}
}
filter->step = (filter->step + 1) % filter->maxstep;
#endif
int s, n, start, stop, flen, flen4;
__m128 *x4_real;
__m128 *x4_imag;
__m128 *h4_real;
__m128 *h4_imag;
__m128 *y4_real;
__m128 *y4_imag;
float *x_real;
float *x_imag;
float *h_real;
float *h_imag;
float *y_real;
float *y_imag;
flen = filter->framelength;
x_real = filter->in_freq_real;
x_imag = filter->in_freq_imag;
x4_real = (__m128*)x_real;
x4_imag = (__m128*)x_imag;
start = filter->steptask[filter->step];
stop = filter->steptask[filter->step + 1];
for (s = start; s < stop; s++)
{
n = (s + filter->mixpos) % filter->num_mixbuf;
y_real = filter->mixbuf_freq_real[n];
y_imag = filter->mixbuf_freq_imag[n];
y4_real = (__m128*)y_real;
y4_imag = (__m128*)y_imag;
h_real = filter->filterbuf_freq_real[s];
h_imag = filter->filterbuf_freq_imag[s];
h4_real = (__m128*)h_real;
h4_imag = (__m128*)h_imag;
flen4 = flen / 4;
for (n = 0; n < flen4; n++)
{
#ifdef WIN32
__m128 a = _mm_mul_ps(x4_real[n], h4_real[n]);
__m128 b = _mm_mul_ps(x4_imag[n], h4_imag[n]);
__m128 c = _mm_sub_ps(a, b);
y4_real[n] = _mm_add_ps(y4_real[n], c);
a = _mm_mul_ps(x4_real[n], h4_imag[n]);
b = _mm_mul_ps(x4_imag[n], h4_real[n]);
c = _mm_add_ps(a, b);
y4_imag[n] = _mm_add_ps(y4_imag[n], c);
#else
y4_real[n] += x4_real[n] * h4_real[n] -
x4_imag[n] * h4_imag[n];
y4_imag[n] += x4_real[n] * h4_imag[n] +
x4_imag[n] * h4_real[n];
#endif
}
y_real[flen] += x_real[flen] * h_real[flen] -
x_imag[flen] * h_imag[flen];
y_imag[flen] += x_real[flen] * h_imag[flen] +
x_imag[flen] * h_real[flen];
}
filter->step = (filter->step + 1) % filter->maxstep;
}
第四步：返回处理的音频流
1.从混合缓冲区 (mixbuf) 取出频域卷积结果，用于逆变换 (IFFT)。
2.执行 IFFT（逆快速傅里叶变换），把数据从 频域转换回时域。
3.应用 Overlap-Add 技术，合并当前帧的输出，使分帧处理后的音频数据平滑衔接。
4.存储本次计算的后半部分到历史缓冲区 (history_time)，用于下一帧计算。
void hcGetSingle(HConvSingle *filter, float *y)
{
int flen, mpos;
float *out;
float *hist;
int size, n, j;
// 获取帧长（flen 表示每帧音频数据的采样点数）
flen = filter->framelength;
// 获取当前混合缓冲区索引（用于从频域数据中提取时域信号）
mpos = filter->mixpos;
// 指向 IFFT 变换后的时域输出缓冲区
out = filter->dft_time;
// 指向历史缓冲区（Overlap-Add 技术使用）
hist = filter->history_time;
// 从混合缓冲区获取频域数据，并清空混合缓冲区
for (j = 0; j < flen + 1; j++)
{
// 复制频域数据到 IFFT 输入缓冲区
filter->dft_freq[j][0] = filter->mixbuf_freq_real[mpos][j]; // 实部
filter->dft_freq[j][1] = filter->mixbuf_freq_imag[mpos][j]; // 虚部
// 清空当前混合缓冲区，避免数据残留影响下一次计算
filter->mixbuf_freq_real[mpos][j] = 0.0;
filter->mixbuf_freq_imag[mpos][j] = 0.0;
}
// 执行 IFFT（逆快速傅里叶变换），将频域数据变回时域
fftwf_execute(filter->ifft);
// 采用 Overlap-Add 技术进行帧间平滑拼接，保证连续性
for (n = 0; n < flen; n++)
{
// 叠加当前 IFFT 变换结果与历史数据
y[n] = out[n] + hist[n];
}
// 存储当前 IFFT 计算结果的后 flen 个点到历史缓冲区
size = sizeof(float) * flen;
memcpy(hist, &(out[flen]), size);
// 更新混合缓冲区索引，确保循环使用不同的缓冲区存储数据
filter->mixpos = (filter->mixpos + 1) % filter->num_mixbuf;
}
单级算法总结：
hcInitSingle（初始化）：加载 脉冲响应（滤波器），进行 分段 FFT 预处理。
hcPutSingle（输入处理）：接收输入数据，转换到 频域（FFT），存入缓冲区。
hcProcessSingle（频域卷积）：执行 频域复数乘法 X(f) * H(f)，累加到混合缓冲区 mixbuf。
hcGetSingle（获取输出）：执行 IFFT（逆变换），恢复 时域音频数据，并应用 Overlap-Add，得到 最终输出音频
双级卷积算法
可以看到双级处理算法在初始化的时候把脉冲响应分成一个短的部分和长的部分分开处理
void hcInitDual(HConvDual *filter, float *h, int hlen, int sflen, int lflen)
{
int size;
float *h2 = NULL;
int h2len;
// sanity check: minimum impulse response length
h2len = 2 * lflen + 1;
if (hlen < h2len)
{
size = sizeof(float) * h2len;
h2 = (float*)fftwf_malloc(size);
memset(h2, 0, size);
size = sizeof(float) * hlen;
memcpy(h2, h, size);
h = h2;
hlen = h2len;
}
// processing step counter
filter->step = 0;
// number of processing steps per long audio frame
filter->maxstep = lflen / sflen;
// number of samples per long audio frame
filter->flen_long = lflen;
// number of samples per short audio frame
filter->flen_short = sflen;
// input buffer (long frame)
size = sizeof(float) * lflen;
filter->in_long = (float *)fftwf_malloc(size);
memset(filter->in_long, 0, size);
// output buffer (long frame)
size = sizeof(float) * lflen;
filter->out_long = (float *)fftwf_malloc(size);
memset(filter->out_long, 0, size);
// convolution filter (short segments)
size = sizeof(HConvSingle);
filter->f_short = (HConvSingle *)malloc(size);
hcInitSingle(filter->f_short, h, 2 * lflen, sflen, 1);
// convolution filter (long segments)
size = sizeof(HConvSingle);
filter->f_long = (HConvSingle *)malloc(size);
hcInitSingle(filter->f_long, &(h[2 * lflen]), hlen - 2 * lflen, lflen, lflen / sflen);
if (h2 != NULL)
{
fftwf_free(h2);
}
}
void hcProcessDual(HConvDual *filter, float *in, float *out)
{
int lpos, size, i;
// convolution with short segments
hcPutSingle(filter->f_short, in);
hcProcessSingle(filter->f_short);
hcGetSingle(filter->f_short, out);
// add contribution from last long frame
lpos = filter->step * filter->flen_short;
for (i = 0; i < filter->flen_short; i++)
out[i] += filter->out_long[lpos + i];
// convolution with long segments
if (filter->step == 0)
hcPutSingle(filter->f_long, filter->in_long);
hcProcessSingle(filter->f_long);
if (filter->step == filter->maxstep - 1)
hcGetSingle(filter->f_long, filter->out_long);
// add current frame to long input buffer
lpos = filter->step * filter->flen_short;
size = sizeof(float) * filter->flen_short;
memcpy(&(filter->in_long[lpos]), in, size);
// increase step counter
filter->step = (filter->step + 1) % filter->maxstep;
}
三级卷积算法
可以从下面代码看出来三级是一个单级和一个双极结合而成的
void hcInitTripple(HConvTripple *filter, float *h, int hlen, int sflen, int mflen, int lflen)
{
int size;
float *h2 = NULL;
int h2len;
// sanity check: minimum impulse response length
h2len = mflen + 2 * lflen + 1;
if (hlen < h2len)
{
size = sizeof(float) * h2len;
h2 = (float*)fftwf_malloc(size);
memset(h2, 0, size);
size = sizeof(float) * hlen;
memcpy(h2, h, size);
h = h2;
hlen = h2len;
}
// processing step counter
filter->step = 0;
// number of processing steps per medium audio frame
filter->maxstep = mflen / sflen;
// number of samples per medium audio frame
filter->flen_medium = mflen;
// number of samples per short audio frame
filter->flen_short = sflen;
// input buffer (medium frame)
size = sizeof(float) * mflen;
filter->in_medium = (float *)fftwf_malloc(size);
memset(filter->in_medium, 0, size);
// output buffer (medium frame)
size = sizeof(float) * mflen;
filter->out_medium = (float *)fftwf_malloc(size);
memset(filter->out_medium, 0, size);
// convolution filter (short segments)
size = sizeof(HConvSingle);
filter->f_short = (HConvSingle *)malloc(size);
hcInitSingle(filter->f_short, h, mflen, sflen, 1);
// convolution filter (medium segments)
size = sizeof(HConvDual);
filter->f_medium = (HConvDual *)malloc(size);
hcInitDual(filter->f_medium, &(h[mflen]), hlen - mflen, mflen, lflen);
if (h2 != NULL)
{
fftwf_free(h2);
}
}
void hcProcessTripple(HConvTripple *filter, float *in, float *out)
{
int lpos, size, i;
// convolution with short segments
hcPutSingle(filter->f_short, in);
hcProcessSingle(filter->f_short);
hcGetSingle(filter->f_short, out);
// add contribution from last medium frame
lpos = filter->step * filter->flen_short;
for (i = 0; i < filter->flen_short; i++)
out[i] += filter->out_medium[lpos + i];
// add current frame to medium input buffer
lpos = filter->step * filter->flen_short;
size = sizeof(float) * filter->flen_short;
memcpy(&(filter->in_medium[lpos]), in, size);
// convolution with medium segments
if (filter->step == filter->maxstep - 1)
hcProcessDual(filter->f_medium,
filter->in_medium,
filter->out_medium);
// increase step counter
filter->step = (filter->step + 1) % filter->maxstep;
}
总结
可以根据输入的脉冲响应来决定使用哪种卷积算法。算法有些细节虽然还是不太能理解，先尝试调用。
我将找找多种脉冲响应，在下一节尝试卷积看看效果。
作者水平有限，文章如果有问题欢迎指正。
扫一扫
1927
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
