如何为长信号和长核卷积编写C代码

Question

我想对长度为4000*270的信号做一个线性卷积，内核长度为16000。信号不固定，内核是固定的。为了我的目的，这需要重复多次，所以我想尽快提高速度。我可以在 R 或 C 中实现这个卷积。

起初，我尝试在R中做卷积，但速度不能满足我的需要。我试着通过迭代来做，但速度太慢了。我也试过用FFT做，但是因为signal和kernel都很长，FFT并没有提升多少速度。

然后我决定在C中迭代地做卷积。但是C似乎无法处理这样的计算量并且经常报错。即使它可以工作，它仍然很慢。我也试过在 C 中做 fft 卷积，但程序总是关闭。

我从我的一个朋友那里找到了这段代码，但不确定原始来源。有版权我就删了 issue.This 是我在C中做fft的C代码，但是程序无法处理长度为2097152的长向量（大于等于信号的2的最小次方）矢量长度）。

#define q    3        /* for 2^3 points */
#define N    2097152        /* N-point FFT, iFFT */

typedef float real;
typedef struct{real Re; real Im;} complex;

#ifndef PI
# define PI    3.14159265358979323846264338327950288
#endif


void fft( complex *v, int n, complex *tmp )
                   {
if(n>1) {            /* otherwise, do nothing and return */
  int k,m;    
  complex z, w, *vo, *ve;
  ve = tmp; 
  vo = tmp+n/2;

  for(k=0; k<n/2; k++) {
          ve[k] = v[2*k];
          vo[k] = v[2*k+1];
      }
  fft( ve, n/2, v );        /* FFT on even-indexed elements of v[] */
  fft( vo, n/2, v );        /* FFT on odd-indexed elements of v[] */
  for(m=0; m<n/2; m++) {
      w.Re = cos(2*PI*m/(double)n);
      w.Im = -sin(2*PI*m/(double)n);
      z.Re = w.Re*vo[m].Re - w.Im*vo[m].Im;    /* Re(w*vo[m]) */
      z.Im = w.Re*vo[m].Im + w.Im*vo[m].Re;    /* Im(w*vo[m]) */
      v[  m  ].Re = ve[m].Re + z.Re;
      v[  m  ].Im = ve[m].Im + z.Im;
      v[m+n/2].Re = ve[m].Re - z.Re;
      v[m+n/2].Im = ve[m].Im - z.Im;
      }
      }
  return;
       }

  void ifft( complex *v, int n, complex *tmp )
          {
    if(n>1) {            /* otherwise, do nothing and return */
           int k,m;    
           complex z, w, *vo, *ve;
           ve = tmp; 
           vo = tmp+n/2;
    for(k=0; k<n/2; k++) {
          ve[k] = v[2*k];
          vo[k] = v[2*k+1];
          }
    ifft( ve, n/2, v );        /* FFT on even-indexed elements of v[] */
    ifft( vo, n/2, v );        /* FFT on odd-indexed elements of v[] */
    for(m=0; m<n/2; m++) {
            w.Re = cos(2*PI*m/(double)n);
            w.Im = sin(2*PI*m/(double)n);
            z.Re = w.Re*vo[m].Re - w.Im*vo[m].Im;    /* Re(w*vo[m]) */
            z.Im = w.Re*vo[m].Im + w.Im*vo[m].Re;    /* Im(w*vo[m]) */
            v[  m  ].Re = ve[m].Re + z.Re;
            v[  m  ].Im = ve[m].Im + z.Im;
            v[m+n/2].Re = ve[m].Re - z.Re;
            v[m+n/2].Im = ve[m].Im - z.Im;
            }
            }
        return;
        }

我发现这个页面在谈论长信号卷积 https://ccrma.stanford.edu/~jos/sasp/Convolving_Long_Signals.html 但我不确定如何在其中使用这个想法。如果有任何想法，我们将不胜感激，我准备好提供有关我的问题的更多信息。

Answer 1

根据您引用的 CCRMA 论文，最常见的高效长 FIR 滤波器方法是使用 FFT/IFFT 重叠添加（或重叠保存）快速卷积。只需将您的数据分成更适合您的 FFT 库和处理器数据缓存大小的更短的块，至少按过滤器内核长度、FFT 过滤器进行零填充，然后在每个 IFFT 之后依次重叠添加 remainder/tails。

非常长的 FFT 很可能会破坏处理器的缓存，这可能会超过任何算法 O(NlogN) 加速。

如何为长信号和长核卷积编写C代码

How to write C code for a long signal and long kernel convolution

c

r

fft

convolution

dft