c - FFT无法重新组合结果？

Question

这个问题不太可能对任何 future 的访客有帮助；它只与一个较小的地理区域、一个特定的时间点或一个非常狭窄的情况相关，通常不适用于全世界的互联网受众。如需帮助使此问题更广泛适用， visit the help center 。

已关闭 10 年前。

所以基本上我有一个函数可以计算给定大小为 N 的 vector 的复离散傅里叶级数。(所以我已经获得了 vector y 并且我必须找到 x)该代码使用 Daniel-Lanczo 的 FFT 算法，并且在 N=2^m 时完美运行。然而，当我用 N=3.2^m 或 N=5.2^m 之类的东西测试它时，我开始遇到问题。到目前为止我考虑过的最简单的是 N=6。

if((N%3)==0 && (N%2)==0 && (N%5)!=0)

是使 N=6 落入该特定 if 语句的条件。在此内部，我首先将 y 的元素分成 y (y_e) 的偶数条目，然后将 y (y_o) 的奇数条目分成并将它们加载到另一个大小为 2*N 的数组 w 中。一旦找到 x，我也会以类似的方式分离出 x_e 和 x_o，并将它们加载到 w 的不同部分。

该算法的工作原理是将 N 分成越来越小的分量，然后对它们进行傅立叶变换。例如，如果 N=4，则它将其分成 2 个大小为 2 的 vector (即 y_e 和 y_o)，现在在其自身内部调用该函数两次，以获得大小 N=2，然后返回结果。然后将结果组合/累积以获得最终的傅立叶级数。

N=6 出现的问题是它把它分成了 2 个，所以我有 2 个大小为 3 的 vector y_e 和 y_o，所以当我现在从自身内部调用该函数时，大小为 3 而不是 6 ，它落入下面的第一个 if 语句(如果 N==3)，并执行直接矩阵乘法来查找级数。请注意，我已在“If N==3”内将此 vector 称为“x”，它现在应该返回到我调用该函数的位置

FastDFS(x, y_e, w, Wp, (N/2), 1);

(其中 FastDFS 只是函数的名称)并再次打印结果“x”。但由于某种原因它未能这样做。它告诉我 x 的条目是 [0+0i,0+0i,0+0i]，即它是空的。我不明白为什么会发生这种情况。我尝试创建一个新 vector ，在找到 N=3 的 x 后为其分配“x”条目，然后在 N=6 的情况下打印新 vector ，它可以工作..但它不切实际不起作用当我尝试更高的 N(例如 12 或 24)时。

有谁知道为什么它可能将 x 的条目设置为零？

我知道这很令人困惑，但如果有人可以提供帮助，我将非常感激!

if(N==3)
{
    Cn=MakeCn(3);

    x=complex_matrix_vector_multiply(Cn, y, 3, 3, 3, 1);
    print_complex_vector(x, 3);

    /*for(i=2*N;i<=3*N-1;i++)
    {
        for(j=0;j<=N-1;j++)
        {
            w[i]=x[j];
            i++;
        }
    }*/
    //printf("\ni am here thuis is w\n");
    //print_complex_vector(w, 12);
}


if((N%3)==0 && (N%2)==0 && (N%5)!=0)
{   
    double complex *x_e,*x_o;
    x_e=make_complex_vector(x_e, N/2);
    x_o=make_complex_vector(x_o, N/2);
    printf("whatup BBBB");
    int j=0,k=0;
    double complex *y_e,*y_o;
    printf("\nthis is N: %d\n",N);
    print_complex_vector(y, N);
    y_e=make_complex_vector(y_e,N/2);
    y_o=make_complex_vector(y_o,N/2);

    /*********************************************************************************************************
     Here were are going to separate the even and odd elements of y into the y_e and y_o. 
     *********************************************************************************************************/

    for(i=0;i<=N-1;i++)
    {
        if(i%2==0)
        {
            y_e[j]=y[i];
            j++;
        }
        else 
        {
            y_o[k]=y[i];
            k++;
        }

    }
    printf("\n These are vectors y_e and y_o: \n");
    print_complex_vector(y_e, N/2);
    print_complex_vector(y_o, N/2);

    /*********************************************************************************************************
     Here were are going to load the even and odd elements of y into the w. w[0...N/2-1]=y_e and w[N/2...N-1]=y_o 
     *********************************************************************************************************/

    for(k=0;k<=(N/2)-1;k++)
    {
        for(i=0;i<=(N-1);i++)
        {

            if(i%2==0)
            {   
                w[k]=y[i];
                k++;
            }
        }
    }

    for(j=N/2;j<=N-1;j++)
    {
        for(i=0;i<=(N-1);i++)
        {

            if(i%2!=0)
            {
                w[j]=y[i];
                j++;
            }   
        }
    }
    printf("\n This is the vector w: \n");
    print_complex_vector(w, N);

    /*********************************************************************************************************
     Here were are going to call FastDFS twice within itself for N/2 with x, y_e, y_o and w. 
     The values of x that are found are the respective x_e and x_o that we load into w. 
     w[N...3N/2-1]=x_e and w[3N/2...2N-1]=x_o
     *********************************************************************************************************/

    Cn2=MakeCn(N/2);

    FastDFS(x, y_e, w, Wp, (N/2), 1);
    // printf("\n w in we ljkdgj\n");
    // print_complex_vector(w, 2*N);
    /* for(i=N;i<=(3*N/2)-1;i++)
     {
     for(j=0;j<=N-1;j++)
     {
     w[i]=x[j];
     i++;
     }
     }*/
    printf("\n here:\n");
    print_complex_vector(x, N);


    for(j=0;j<=(N/2)-1;j++)
    {
        for(i=N;i<=(3*N/2)-1;i++)
        {
            x_e[j]=w[i];
            j++;
        }   
    }
    printf("\n this is x_e:\n");
    print_complex_vector(x_e, N/2);

    FastDFS(x, y_o, w, Wp, (N/2), 1);

    // print_complex_vector(w, 2*N);
    for(j=0;j<=(N/2)-1;j++)
    {
        for(i=N;i<=(3*N/2)-1;i++)
        {
            x_o[j]=w[i];
            j++;
        }   
    }
    printf("\n this is x_o:\n");
    print_complex_vector(x_o, N/2);

    for(k=N;k<=(3*N/2)-1;k++)
    {
        for(j=0;j<=(N/2)-1;j++)
        {

            w[k]=x_e[j];
            k++;
        }
    }

    for(k=(3*N/2);k<=2*N-1;k++)
    {
        for(j=0;j<=(N/2)-1;j++)
        {

            w[k]=x_o[j];
            k++;
        }
    }
    print_complex_vector(w, 2*N);

    /*********************************************************************************************************
     Here were are going to find the final x_j and x_j+N/2 by x_j=[x_e]+W_n^j[x_o] and x_j+N/2+Wn^j+N/2[x_o]
     This is the final answer for the Discrete Fourier Series of N even.
     We do not use x_e and x_o explicitly but use different parts of w. 
     *********************************************************************************************************/

    F=cos(2*pi/N)+I*sin(2*pi/N);


    for(i=0;i<=(N/2)-1;i++)
    {
        for(j=N;j<=2*N-1;j++)
        {   
            x[i]=w[j]+((cpow(F, i))*w[(j+N/2)]);
            i++;
        }
    }

    for(i=0;i<=(N/2)-1;i++)
    {
        for(k=N;k<=(2*N)-1;k++)
        {
            x[(i+N/2)]=w[k]+((cpow(F, (i+N/2)))*w[k+N/2]);
            i++;
        }
    }

    /*for(i=0;i<=(N/2)-1;i++)
     {

     x[i]=x_e[i]+((cpow(F, i))*x_o[i]);
     x[(i+N/2)]=x_e[i]+((cpow(F, (i+N/2)))*x_o[i]);


     }
     */

    printf("\n\nThe Final Discrete Fourier Series, for N = %d, is:\n\n",N);

    for (i=0; i<=N-1; i++) 

    {
        printf ("%9.4f+%.4fi  \n", creal(x[i]),cimag(x[i]));

    }


}

Answer 1

线路

x=complex_matrix_vector_multiply(Cn, y, 3, 3, 3, 1);

将用计算结果覆盖x的本地值。我假设 x 是一个指针，并且最初该指针指示您希望结果到达的位置。但在这一行之后，它指向其他位置，并且您最初传入的 x 将无法从函数内访问。

您可能想要使用某种形式的乘法，将其输出写入预先分配的数组。或者您想使用临时指针存储结果，然后一次从该结果复制到 x 一个元素。请记住之后释放临时结果。

将函数头添加到您粘贴的代码片段中可能会使这些事情变得更加清晰，并且对于将来可能想要阅读它的人来说仍然会使问题更加清晰。

我没有仔细查看其余的实现，但是 Cn2=MakeCn(N/2); 行看起来是多余的。您没有使用该值，是吗？