c - FFTW 相关中的交替误差

Question

// Initialize variables, ZALLOC allocates and fills with zeroes
ZALLOC(z0q, numqbins, double) ;  
ZALLOC(z1q, numqbins, double) ;  
ZALLOC(z2q, numqbins, double) ;  

z_in = (double*)fftw_malloc(sizeof(double) * numqbins * 2) ;
z_out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * numqbins * 2) ;

r2c_plan = fftw_plan_dft_r2c_1d(2 * numqbins, z_in, z_out, FFTW_ESTIMATE) ;
c2r_plan = fftw_plan_dft_c2r_1d(2 * numqbins, z_out, z_in, FFTW_ESTIMATE) ;


... (ddbins is initialized to all 0s, z0q, z1q, and z2q are filled with values) ...


for (k1=0; k1<numx; ++k1) { 
    if (snpmarkers[xindex[k1]] -> ignore) {
        continue ;
    } else {
        s = snpmarkers[xindex[k1]] -> tagnumber ;
    }
    y = res[k1] * wt[k1] ; 

    z0q[s] += 1 ; 
    z1q[s] += y ; 
    z2q[s] += y*y ; 
  }
// Perform convolution.
fft_ddadd(ddbins, z0q, z1q, z2q, numqbins, diffmax, z_in, z_out, r2c_plan, c2r_plan) ;

void fft_ddadd(double **ddbins, double *z0q, double *z1q, double *z2q, int numqbins, 
            int diffmax, double *z_in, fftw_complex *z_out, fftw_plan r2c_plan, fftw_plan c2r_plan) 
{
    double *dd00, *dd01, *dd10, *dd11, *dd02, *dd20 ; 
    int i;
    fftw_complex *ft_z0q, *ft_z1q, *ft_z2q ;    

    dd00 = ddbins[0] ;
    dd01 = ddbins[1] ; 
    dd10 = ddbins[3] ; 
    dd11 = ddbins[4] ; 
    dd02 = ddbins[2] ; 
    dd20 = ddbins[6] ;     

    ft_z0q = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * numqbins * 2) ;
    ft_z1q = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * numqbins * 2) ;
    ft_z2q = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * numqbins * 2) ;    


    // Copy arrays into z_in and transform, copy complex results out into ft_z's
    for (i = 0; i < numqbins; ++i) {
        z_in[i] = z0q[i] ;
        z_in[i + numqbins] = 0 ;
    }    

    fftw_execute(r2c_plan) ;    

    for (i = 0; i < numqbins; ++i) {
        z_in[i] = z1q[i] ;
        ft_z0q[i][0] = z_out[i][0] ;
        ft_z0q[i][1] = z_out[i][1] ;
    }    

    fftw_execute(r2c_plan) ;    

    for (i = 0; i < numqbins; ++i) {
        z_in[i] = z2q[i] ;
        ft_z1q[i][0] = z_out[i][0] ;
        ft_z1q[i][1] = z_out[i][1] ;
    }    

    fftw_execute(r2c_plan) ;    

    for (i = 0; i < numqbins; ++i) {
        ft_z2q[i][0] = z_out[i][0] ;
        ft_z2q[i][1] = z_out[i][1] ;
    }    

    // Do correlations
    fft_singledd(dd00, ft_z0q, ft_z0q, numqbins, diffmax, z_in, z_out, c2r_plan);
    fft_singledd(dd01, ft_z0q, ft_z1q, numqbins, diffmax, z_in, z_out, c2r_plan);
    fft_singledd(dd10, ft_z1q, ft_z0q, numqbins, diffmax, z_in, z_out, c2r_plan);
    fft_singledd(dd11, ft_z1q, ft_z1q, numqbins, diffmax, z_in, z_out, c2r_plan);
    fft_singledd(dd02, ft_z0q, ft_z2q, numqbins, diffmax, z_in, z_out, c2r_plan);
    fft_singledd(dd20, ft_z2q, ft_z0q, numqbins, diffmax, z_in, z_out, c2r_plan);    


}    

void fft_singledd(double* ddbin, fftw_complex* ftz_0, fftw_complex* ftz_1, int numqbins, 
            int diffmax, double *z_in, fftw_complex *z_out, fftw_plan c2r_plan) {

    int i ;
    // Correlate ftz_0 and ftz_1 together and copy into z_out
    for (i = 0; i < numqbins * 2; ++i) {
        z_out[i][0] =   ftz_0[i][0] * ftz_1[i][0] + ftz_0[i][1] * ftz_1[i][1] ;
        z_out[i][1] = - ftz_0[i][0] * ftz_1[i][1] + ftz_0[i][1] * ftz_1[i][0] ;
    }    
    // Transform back
    fftw_execute(c2r_plan) ;    
    // Copy the desired results out
    for (i = 1; i < diffmax; ++i) {
        ddbin[i] += (z_in[2 * numqbins - i] / (2 * numqbins)) ;
    }    
}

void ddadd(double **ddbins, double *z0q, double *z1q, double *z2q, int numqbins, int diffmax) {
    double *dd00, *dd01, *dd10, *dd11, *dd02, *dd20 ; 
    double a0, a1, a2, b0, b1, b2 ; 
    int i, s, t, tmax, d ;    

    dd00 = ddbins[0] ;
    dd01 = ddbins[1] ; 
    dd10 = ddbins[3] ; 
    dd11 = ddbins[4] ; 
    dd02 = ddbins[2] ; 
    dd20 = ddbins[6] ;     

    for (s=0; s<numqbins; ++s) { 
        a0 = z0q[s] ; 
        a1 = z1q[s] ;
        a2 = z2q[s] ;
        tmax = MIN(numqbins-1, s + diffmax) ; 
        for (t=s+1; t<=tmax; ++t) { 
            b0 = z0q[t] ;
            b1 = z1q[t] ;
            b2 = z2q[t] ;
            d = t - s ; 
            dd00[d] += a0*b0 ;
            dd01[d] += a0*b1 ;
            dd10[d] += a1*b0 ;
            dd11[d] += a1*b1 ;
            dd02[d] += a0*b2 ;
            dd20[d] += a2*b0 ;
        }
    }
}

运行上述代码时，我在结果中遇到了问题。当将输出ddbins（ddadd非傅立叶变换，由 O(n^2) 相关性计算）与 的输出进行比较时fft_ddadd，误差几乎是恒定的并且是交替的：

ddbins[0] - fft_ddbins[0] = [0.000000000000, -0.016764512518, 0.016764511936, -0.016764512286, 0.016764512053, -0.016764512286, 0.016764512053, -0.016764512169, 0.016764511995, -0.016764512227, 0.016764512053, -0.016764512344, 0.016764512053, -0.016764512227, 0.016764511995, -0.016764512227, 0.016764512053, -0.016764512344, 0.016764512053, -0.016764512227, 0.016764512053, -0.016764512344, 0.016764511995, -0.016764512227, 0.016764512053, -0.016764512227, 0.016764512053, -0.016764512344 ...]

如您所见，每个地方的误差大致相同，符号交替，加上或减去一些归因于机器精度的误差。

已经提出的一个可能的问题是，存在一个傅里叶变换分量，它在实际空间中对应于数组[+1, -1, +1, -1 ...]，如果排除，它将解释交替误差。但是，我看不出如何排除这个元素。数组 z*q 被完全复制 (length numqbins) 到计划中并执行，而输出数组 (length 2*numqbins) 也被直接复制，并且在输出复制到所需位置之前引入了错误。

这是我第一次在这个网站上提问，所以请告诉我是否需要更多信息（或更少！），或者如果这不是正确的提问地方。提前致谢。

标签： cfftfftw