有效的比特流卷積

Question

我有兩個浮點時間序列A，B各有長度N。我必須計算循環卷積並找到最大值。這樣做的經典和最快的方式是

C = IFFT（FFT（A）* FFT（B））

現在，讓我們假設A和B是一系列只包含1和0，所以原則上我們可以將它們表示爲比特流。

問題：如果我能以某種方式利用上述事實，是否有更快的卷積方法（並找到其最大值）？

（我已經思考沃爾什很多 - 阿達瑪變換和SSE指令，popcounts，但發現M> 2 ** 20沒有更快的方法這是我的情況。）

感謝， GD

Answer 1

一維卷積兩個陣列a和大小n的b的c是一個數組，使得：

之和的非空的術語不限於變化的bnb的數目：如果b是一個簡單的圖案，這個總和可這個公式可以以迭代的方式被重寫限於幾個條款。的算法現在可以設計用以計算c：

1：計算c[0]（約n操作）

2：使用式（約nb*n操作）

如果nb小0<i<n計算c[i]，這個方法可能比fft更快。請注意，它將爲比特流信號提供精確的結果，而fft需要過採樣和浮點精度來提供準確的結果。

下面是一段代碼，用輸入類型unsigned char實現這個技巧。

#include <stdlib.h> 
#include <math.h> 
#include <string.h> 
#include <time.h> 

#include <fftw3.h> 

typedef struct{ 
    unsigned int nbchange; 
    unsigned int index[1000]; 
    int change[1000]; 
}pattern; 

void topattern(unsigned int n, unsigned char* b,pattern* bp){ 
    //initialisation 
    bp->nbchange=0; 
    unsigned int i; 
    unsigned char former=b[n-1]; 
    for(i=0;i<n;i++){ 
     if(b[i]!=former){ 
      bp->index[bp->nbchange]=i; 
      bp->change[bp->nbchange]=((int)b[i])-former; 
      bp->nbchange++; 
     } 
     former=b[i]; 
    } 
} 

void printpattern(pattern* bp){ 
    int i; 
    printf("pattern :\n"); 
    for(i=0;i<bp->nbchange;i++){ 
     printf("index %d change %d\n",bp->index[i],bp->change[i]); 
    } 
} 

//https://stackoverflow.com/questions/109023/how-to-count-the-number-of-set-bits-in-a-32-bit-integer 

unsigned int NumberOfSetBits(unsigned int i) 
{ 
    i = i - ((i >> 1) & 0x55555555); 
    i = (i & 0x33333333) + ((i >> 2) & 0x33333333); 
    return (((i + (i >> 4)) & 0x0F0F0F0F) * 0x01010101) >> 24; 
} 

//https://stackoverflow.com/questions/2525310/how-to-define-and-work-with-an-array-of-bits-in-c 

unsigned int convol_longint(unsigned int a, unsigned int b){ 
    return NumberOfSetBits(a&b); 
} 

int main(int argc, char* argv[]) { 

    unsigned int n=10000000; 

    //the array a 
    unsigned char* a=malloc(n*sizeof(unsigned char)); 
    if(a==NULL){printf("malloc failed\n");exit(1);} 
    unsigned int i,j; 
    for(i=0;i<n;i++){ 
     a[i]=rand(); 
    } 
    memset(&a[2],5,2); 
    memset(&a[10002],255,20); 

    for(i=0;i<n;i++){ 
     //printf("a %d %d \n",i,a[i]); 
    } 

    //pattern b 
    unsigned char* b=malloc(n*sizeof(unsigned char)); 
    if(b==NULL){printf("malloc failed\n");exit(1);} 
    memset(b,0,n*sizeof(unsigned char)); 
    memset(&b[2],1,20); 


    //memset(&b[120],1,10); 
    //memset(&b[200],1,10); 

    int* c=malloc(n*sizeof(int)); //nb bit in the array 
    memset(c,0,n*sizeof(int)); 

    clock_t begin, end; 
    double time_spent; 

    begin = clock(); 
    /* here, do your time-consuming job */ 


    //computing c[0] 
    for(i=0;i<n;i++){ 
     //c[0]+= convol_longint(a[i],b[i]); 
     c[0]+= ((int)a[i])*((int)b[i]); 
     //printf("c[0] %d %d\n",c[0],i); 
    } 
    printf("c[0] %d\n",c[0]); 

    //need to store b as a pattern. 
    pattern bpat; 
    topattern(n,b,&bpat); 
    printpattern(&bpat); 

    //computing c[i] according to formula 
    for(i=1;i<n;i++){ 
     c[i]=c[i-1]; 
     for(j=0;j<bpat.nbchange;j++){ 
      c[i]+=bpat.change[j]*((int)a[(bpat.index[j]-i+n)%n]); 
     } 
    } 

    //finding max 
    int currmax=c[0]; 
    unsigned int currindex=0; 
    for(i=1;i<n;i++){ 
     if(c[i]>currmax){ 
      currmax=c[i]; 
      currindex=i; 
     } 
     //printf("c[i] %d %d\n",i,c[i]); 
    } 

    printf("c[max] is %d at index %d\n",currmax,currindex); 

    end = clock(); 
    time_spent = (double)(end - begin)/CLOCKS_PER_SEC; 

    printf("computation took %lf seconds\n",time_spent); 


    double* dp = malloc(sizeof (double) * n); 
    fftw_complex * cp = fftw_malloc(sizeof (fftw_complex) * (n/2+1)); 

    begin = clock(); 
    fftw_plan plan = fftw_plan_dft_r2c_1d(n, dp, cp, FFTW_ESTIMATE); 

    end = clock(); 
    time_spent = (double)(end - begin)/CLOCKS_PER_SEC; 

    fftw_execute (plan); 
    printf("fftw took %lf seconds\n",time_spent); 

    free(dp); 
    free(cp); 

    free(a); 
    free(b); 
    free(c); 
    return 0; 
} 

編譯：gcc main.c -o main -lfftw3 -lm

對於n=10 000 000和nb=2（b只是一個 「矩形1D窗口」）在我的電腦上0.65秒該算法運行。使用fftw的雙精度fft大約花費同一時間。與大多數比較相比，這種比較可能是不公平的，因爲：

• nb=2是此答案中提供的算法的最佳情況。
• 基於fft的算法會需要過採樣。
• 基於fft的算法可能不需要雙精度
• 此處公開的實現未優化。這只是基本的代碼。

該實現可以處理n=100 000 000。此時，建議使用long int for c來避免任何溢出風險。

如果信號是比特流，則可以用各種方式優化該程序。對於按位運算，請看this question和this one