CODEWITHSUNDEEP
c++fft
ectuser
ectuser

Reputation: 171

What could parametrs of FFT function mean

I'm trying to understand the FFT algorithm. Here's a code 

void fft(double *a, double *b, double *w, int m, int l)
{
    int i, i0, i1, i2, i3, j;
    double u, v, wi, wr;
    for (j = 0; j < l; j++) {
        wr = w[j << 1];
        wi = w[j << 1 + 1];
        for (i = 0; i < m; i++) {
            i0 = (i << 1) + (j * m << 1);
            i1 = i0 + (m * l << 1);
            i2 = (i << 1) + (j * m << 2);
            i3 = i2 + (m << 1);
            u = a[i0] - a[i1];
            v = a[i0 + 1] - a[i1 + 1];
            b[i2] = a[i0] + a[i1];
            b[i2 + 1] = a[i0 + 1] + a[i1 + 1];
            b[i3] = wr * u - wi * v;
            b[i3 + 1] = wr * v + wi * u;
        }
    }
}

If I get it right, array W is input, where every odd number is real and even is imag. A and B are imag and real parts of complex result Also I found that l = 2**m

But when i'm trying to do this:

double a[4] = { 0, 0, 0, 0 };
double b[4] = { 0, 0, 0, 0 };
double w[8] = { 1, 0, 0, 0, 0, 0, 0, 0 };
int m = 3;
int l = 8;

fft(a, b, w, m, l);

There's error.

Upvotes: 1

Views: 57

Answers (1)

Eric Postpischil
Eric Postpischil

Reputation: 222933

This code is only part of an FFT. a is input. b is output. w contains precomputed weights. l is a number of subdivisions at the current point in the FFT. m is the number of elements per division. The data in a, b, and w is interleaved complex data—each pair of double elements from the array consists of the real part and the imaginary part of one complex number.

The code performs one radix-two butterfly pass over the data. To use it to compute an FFT, it must be called multiple times with specific values for l, m, and the weights in w. Since, for each call, the input is in a and the output is in b, the caller must use at least two buffers and alternate between them for successive calls to the routine.

From the indexing performed in i0 and i2, it appears the data is being rearranged slightly. This may be intended to produce the final results of the FFT in “natural” order instead of the bit-reversed order that occurs in a simple implementation.

But when i'm trying to do this:

double a[4] = { 0, 0, 0, 0 };
double b[4] = { 0, 0, 0, 0 };
double w[8] = { 1, 0, 0, 0, 0, 0, 0, 0 };
int m = 3;
int l = 8;
 
fft(a, b, w, m, l);

There's error.

From for (j = 0; j < l; j++), we see the maximum value of j in the loop is l-1. From for (i = 0; i < m; i++), we see the maximum value of i is m-1. Then in i0 = (i << 1) + (j * m << 1), we have i0 = ((m-1) << 1) + ((l-1) * m << 1) = (m-1)*2 + (l-1) * m * 2 = 2*m - 2 + l*m*2 - m*2 = 2*m*l - 2. And in i1 = i0 + (m * l << 1), we have i1 = 2*m*l - 2 + (m * l * 2) = 4*m*l - 2. When the code uses a[i1 + 1], the index is i1 + 1 = 4*m*l - 2 + 1 = 4*m*l - 1.

Therefore a must have an element with index 4*m*l - 1, so it must have at least 4*m*l elements. The required size for b can be computed similarly and is the same.

When you call fft with m set to 3 and l set to 8, a must have 4•3•8 = 96 elements. Your sample code shows four elements. Thus, the array is overrun, and the code fails.

I do not believe it is correct that l should equal 2m. More likely, 4*m*l should not vary between calls to fft in the same complete FFT computation, and, since a and b contain two double elements for every complex number, 4*m*l should be twice the number of complex elements in the signal being transformed.

Upvotes: 4

Related Questions