topic FFT CCS format - 1D Fourier series in Intel® oneAPI Math Kernel Library

FFT CCS format - 1D Fourier series

Jozsef — Mon, 09 Aug 2010 19:52:33 GMT

Hello,

I am computing the FFT of a real-valued 1D function in N discrete points. After computing the coefficients, stored in CCS format in array "y", I'd like to evaluate the Fourier series with something like:

double Y;

for ( i = 0; i < 1000; i++ ) {
Y = y[0]/N;
for ( j = 1; j <= N/2; j++ )
Y += y[j*2+0]*2.0/N*cos(2.0*PI*j*i/n) - y[j*2+1]*2.0/N*sin(2.0*PI*j*i/n);
for ( j = N/2+1; j < N; j++ )
Y += y[(N-j)*2+0]*2.0/N*cos(2.0*PI*j*i/n) + y[(N-j)*2+1]*2.0/N*sin(2.0*PI*j*i/n);

// approximated function value now in Y ...
}

As I understand, the CCS format stores the half of the symmetric conjugate-even sequence. Accordingly, I do two loops to get the full sequence for one value of the approximate function. However, I get very oscillatory results, so obviously something is wrong with the scaling and/or with the indexing of the coefficient-array.

Could someone point me to the right direction?

Thanks,
Jozsef

FFT CCS format - 1D Fourier series

barragan_villanueva_ — Tue, 10 Aug 2010 09:18:52 GMT

Hi,

For N even: CCS-result contains

FFT Real	0	1	2	3	...	n-2	n-1	n	n+1
CCS	R₀	0	R₁	I₁	...	R_n/2-1	I_n/2-1	R_n/2	0

For N odd (N=s*2 + 1): CCS-result contains

FFT Real	0	1	2	3	...	n-4	n-3	n-2	n-1	n	n+1
CCS	R₀	0	R₁	I₁	...	I_s-2	R_s-1	I_s-1	R_s	I_s

So, CCS-format has N+2 elements which should be taken into account.

As to your function, I haveseveral questions:
1) why Y = y[0]/N; is inside of i-loop? So in every iteratin Y value is reset...
2) why N and n ? But, N is number of points...

FFT CCS format - 1D Fourier series

Jozsef — Tue, 10 Aug 2010 14:30:41 GMT

Victor,

1) why Y = y[0]/N; is inside of i-loop? So in every iteratin Y value is reset...

The i-loop goes over the domain on which I am evaluating the series: n=1000 points, so it's nice and smooth to plot. The i-loop should look like (my mistake)

for ( i = 0; i < n; i++ )

2) why N and n ? But, N is number of points...

N would be the number of discrete points of the finite series approximation.

-------
Let me try to give a simple example of what I'm trying to do:

Let the input function be sin(x) of which I'm evaluating at N=4 discrete points over the interval x=[0...1), to get the discrete sequence:

double y[N+2];

for ( i = 0; i < N; i++ )
y = sin(2.0*PI*i/N);

which fills the first N elements of the array as y[] = {0.0, 1.0, 0.0, -1.0}.

Then, I do an in-place forward transform and I get the CCS array:

y[] = {0.0, 0.0, -0.0, -2.0, 0.0, 0.0},

which is to be interpreted as the complex sequence of (N+2)/2 = 3 complex numbers as

Z[] ={ 0.00+0.00i, -0.00-2.00i, 0.00+0.00i}.

Now, if I understand correctly, the above array Z contains:
- the first complex coefficient: R0 + 0.0i = 0.0 + 0.0i,
- the first half of the full sequence, which for N=4 is: R1 + I1 = -0.0 - 2.0i and R2 + I2 = 0.0 + 0.0i,
- plus the (would-be-redundant) symmetric conjugate half "sequence": R3 + I3 = -0.0 + 2.0i.

I think so far my interpretation is correct and this is the correct result.

Now, my question is how to interpret (and use) this sequence to evaluate the first N members of the sum in the discrete Fourier series:

f(x) ~ sin(x) = sum_{j=0}^{N-1} Z_j * exp(2 * PI * i * j * x).

(I'm evaluating the x in n=1000 discrete points with the i-loop just to plot f(x).)

Maybe j should go from -N/2 to N/2-1 ?

Thanks for your help,
Jozsef

FFT CCS format - 1D Fourier series

Jozsef — Tue, 10 Aug 2010 20:46:38 GMT

Oh ok. I figured it out in the end. I didn't exactly understand how I'm supposed to traverse the half-sequence and from which direction. Here is a piece of code that does what I needed:

// suppose result of FFT in y[] in CCS format

// now evaluate the discrete Fourier series at n locations

int j, k;
double Y;

for ( j = 0; j < n; j++ )
{
Y =0.0;
for ( k = -N/2; k < N/2; k++ )
{
Y += y[abs(k)*2+0]*cos(2.0*PI*j/n*k) - sgn(k)*y[abs(k)*2+1]*sin(2.0*PI*j/n*k);
// now function value in Y
}
}

Really pretty simple, once one understands it...

Cheers,
Jozsef