so basically I have a function that calculates the complex discrete Fourier series of a given vector of size N. (So I have been provided with the vector y and I have to find x)
The code uses Daniel-Lanczo's algorithm for FFT and works perfectly when N=2^m. However, when I test it with something like N=3.2^m or N=5.2^m, I start to encounter problems. The simplest one I've considered so far is N=6.
if((N%3)==0 && (N%2)==0 && (N%5)!=0)
is the condition that makes N=6 fall into that specific if statement. Inside this, I first separate the elements of y into even entries of y (y_e) and then odd entries of y (y_o) and load these into another array w of size 2*N. Once x is found I also separate out the x_e and x_o in a similar way and load them into different parts of w also.
The algorithm works by splitting N into smaller and smaller components then doing the fourier transform of those. for eg if N=4 then it splits it into 2 vectors of size 2 (namely y_e and y_o), calls the function twice within itself now for size N=2 and then returns the result. The result is then combined/built up to get the final fourier series.
The problem that is turning out for N=6 is that it splits it up into 2, so I have 2 vectors y_e and y_o of size 3, so when I call the function from within itself now with size 3 instead of 6, it falls into the the first if statement below (if N==3) and does a direct matrix multiply to find the series. Notice that I have called this vector 'x' inside "If N==3" It is now supposed to return to the position where I called the function
FastDFS(x, y_e, w, Wp, (N/2), 1);
(where FastDFS is just the name of the function) and print out the result 'x' again. But for some reason it fails to do so. It tells me the entries of x are [0+0i,0+0i,0+0i], i.e. it is empty. I don't understand why this could be happening. I tried to create a new vector, assign it the entries of 'x' after it had found x for N=3, and then printed the new vector inside the N=6 case and it works.. but its impractical doesn't work for when I try higher N like 12 or 24.
Does anyone know why it might be setting the entries of x to zero?
I understand this is confusing but if anyone can help I would really appreciate it!
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]));
}
}
المحلول
The line
x=complex_matrix_vector_multiply(Cn, y, 3, 3, 3, 1);
will overwrite the local value of x with the result of the computation. I assume that x is a pointer, and that initially that pointer indicated the location of where you want the result to go. But after this line, it points to some other location, and the x you originally passed in will become unaccessible from within the function.
You might want to either use some form of the multiplication which will write its output to a pre-allocated array. Or you want to store the result using a temporary pointer, and then copy from that result to x one element at a time. Remember to free the temporary result afterwards.
Adding the function header to the code snippet you pastet might have made such things clearer, and will still make the question clearer to those who might want to read it in the future.
I didn't look too closely at the rest of the implementation, but the line Cn2=MakeCn(N/2); looks superfluous. You're not using that value, are you?
https://stackoverflow.com/questions/16268933
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
Russian