所以基本上我有一个函数可以计算给定大小为 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]));
}
}