ANSI C代码中的1d线性卷积？
我想知道是否有人可以将我推荐给ANSI C中的一维线性卷积代码片段而不是重新发明轮子？ 我在谷歌搜索和堆栈溢出，但在CI中找不到任何东西可以使用。
例如，对于数组A，B和C，所有双精度，其中A和B是输入，C是输出，长度分别为len_A ， len_B和len_C = len_A + len_B - 1 。
我的arrays尺寸很小，因此不需要通过FFT实现快速卷积的任何速度增加。 寻找简单的计算。
这是如何做：
#include #include void convolve(const double Signal[/* SignalLen */], size_t SignalLen, const double Kernel[/* KernelLen */], size_t KernelLen, double Result[/* SignalLen + KernelLen - 1 */]) { size_t n; for (n = 0; n < SignalLen + KernelLen - 1; n++) { size_t kmin, kmax, k; Result[n] = 0; kmin = (n >= KernelLen - 1) ? n - (KernelLen - 1) : 0; kmax = (n < SignalLen - 1) ? n : SignalLen - 1; for (k = kmin; k <= kmax; k++) { Result[n] += Signal[k] * Kernel[n - k]; } } } void printSignal(const char* Name, double Signal[/* SignalLen */], size_t SignalLen) { size_t i; for (i = 0; i < SignalLen; i++) { printf("%s[%zu] = %f\n", Name, i, Signal[i]); } printf("\n"); } #define ELEMENT_COUNT(X) (sizeof(X) / sizeof((X)[0])) int main(void) { double signal[] = { 1, 1, 1, 1, 1 }; double kernel[] = { 1, 1, 1, 1, 1 }; double result[ELEMENT_COUNT(signal) + ELEMENT_COUNT(kernel) - 1]; convolve(signal, ELEMENT_COUNT(signal), kernel, ELEMENT_COUNT(kernel), result); printSignal("signal", signal, ELEMENT_COUNT(signal)); printSignal("kernel", kernel, ELEMENT_COUNT(kernel)); printSignal("result", result, ELEMENT_COUNT(result)); return 0; }
输出：
signal[0] = 1.000000 signal[1] = 1.000000 signal[2] = 1.000000 signal[3] = 1.000000 signal[4] = 1.000000 kernel[0] = 1.000000 kernel[1] = 1.000000 kernel[2] = 1.000000 kernel[3] = 1.000000 kernel[4] = 1.000000 result[0] = 1.000000 result[1] = 2.000000 result[2] = 3.000000 result[3] = 4.000000 result[4] = 5.000000 result[5] = 4.000000 result[6] = 3.000000 result[7] = 2.000000 result[8] = 1.000000
没有经过测试，但似乎它会起作用……
void conv(const double v1[], size_t n1, const double v2[], size_t n2, double r[]) { for (size_t n = 0; n < n1 + n2 - 1; n++) for (size_t k = 0; k < max(n1, n2); k++) r[n] += (k < n1 ? v1[k] : 0) * (n - k < n2 ? v2[n - k] : 0); }
提示： 如果重新发明轮子所需的时间少于找到轮子，请考虑前者。
因为，我们正在对2个有限长度序列进行卷积，因此如果执行循环卷积而不是线性卷积，则实现期望的频率响应。 循环卷积的一个非常简单的实现将获得与Alex给出的算法相同的结果。
#define MOD(n, N) ((n<0)? N+n : n) ...... ...... for(n=0; n < signal_Length + Kernel_Length - 1; n++) { out[n] = 0; for(m=0; m < Kernel_Length; m++) { out[n] = h[m] * x[MOD(nm, N)]; } }
我使用了@ Mehrdad的方法，并创建了以下anwer：
void conv(const double v1[], size_t n1, const double v2[], size_t n2, double r[]) { for (size_t n = 0; n < n1 + n2 - 1; n++) for (size_t k = 0; k < max(n1, n2) && n >= k; k++) r[n] += (k < n1 ? v1[k] : 0) * (n - k < n2 ? v2[n - k] : 0); }
当第二个循环k变得大于n时，索引超出下限的问题，因此，猜测应该有额外的条件来防止这种情况。