// Illustrate circular convolution thru fft. #include #include typedef double _Complex c; // In this code ~ (which is not -) is the complex conjugate operator. // http://gcc.gnu.org/onlinedocs/gcc-4.8.2/gcc/Complex.html#Complex void fft(const int n, c a[n]){ {int j=0, i; for(i=0; i void pc(c a){printf("%16.12f + %16.12fi", creal(a), cimag(a));} #define D (1<<5) int main(){c a[D], b[D], c[D]; {int j=D; while(j--) {a[j] = (j==5 || j==29); b[j] = (j==2 || j==11);}} fft(D, a); fft(D, b); {int j=D; while(j--) c[j] = ~(a[j]*b[j])/D;} fft(D, c); {int j; for(j=0; j 1.e-8) {pc(c[j]); printf(" %d\n", j);}} return 0;} // c[k] = sum[i+j=k mod D] a[i]*b[j] ; convolution /* 1.000000000000 + -0.000000000000i 7 1.000000000000 + -0.000000000000i 8 1.000000000000 + -0.000000000000i 16 1.000000000000 + -0.000000000000i 31 */ // Correlation: // Substituting (~a[j])*b[j] for ~(a[j]*b[j]) is for c[k] = sum[i-j=k mod D] a[i]*b[j] /* 1.000000000000 + 0.000000000000i 3 1.000000000000 + 0.000000000000i 18 1.000000000000 + 0.000000000000i 26 1.000000000000 + 0.000000000000i 27 */