0

所以基本上我有一个函数可以计算给定大小为 N 的向量的复杂离散傅立叶级数。(所以我得到了向量 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 的向量(即 y_e 和 y_o),现在在自身内部调用该函数两次,大小为 N=2,然后返回结果。然后将结果组合/建立以获得最终的傅立叶级数。

N = 6的问题是它把它分成2个,所以我有2个大小为3的向量y_e和y_o,所以当我现在用大小3而不是6从自身内部调用函数时,它会下降进入下面的第一个 if 语句(如果 N==3)并进行直接矩阵乘法以找到系列。请注意,我在 "If N==3" 中调用了这个向量 'x' 它现在应该返回到我调用函数的位置

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

(其中 FastDFS 只是函数的名称)并再次打印出结果“x”。但由于某种原因,它没有这样做。它告诉我 x 的条目是 [0+0i,0+0i,0+0i],即它是空的。我不明白为什么会发生这种情况。我尝试创建一个新向量,在找到 N=3 的 x 后为其分配“x”的条目,然后在 N=6 的情况下打印新向量并且它可以工作.. 但它不切实际不起作用因为当我尝试更高的 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]));

    }


}
4

1 回答 1

0

线

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

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

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

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

我没有仔细研究其余的实现,但这条线Cn2=MakeCn(N/2);看起来是多余的。你没有使用那个值,是吗?

于 2013-04-29T05:17:59.380 回答