2

Standard non-recursive radix-2 algorithm. I realise what I have here isn't optimal for performance (like repeated trig calls), I just want to get it right before optimising and porting to VHDL. I'm using the complex.h header.

I get correct results for arrays of length 8, but not for 16 or 32. HOWEVER, for real signals, we still get the correct result of conjugal symmetry, just the values are wrong. This leads me to believe that the sorting is correct, but maybe there's something wrong with the trig calls for smaller values? Can anyone help?

I apologise that my code isn't fantastically commented, but honestly it wouldn't help much for the FFT algorithm.

void FFT(double complex x[], int N){

    //bit reverse

int s = log2(N);

for(int i = 0; i < N/2; i++){
    int h = bitrev(i, s);

    printf("%u ", h);
    double complex temp = x[i];
    x[i] = x[h];
    x[h] = temp;
}

unsigned int Np = 2; //NUM POINTS IN EACH BLOCK. INITIALLY 2
unsigned int Bp = (N/2);
for(int i = 0; i < s; i++){
    int NpP = Np>>1; //num butterflies
    int BaseT = 0;
    for (int j = 0; j < Bp; j++){

        int BaseB = BaseT + NpP;

        for(int k = 0; k < NpP; k++){

            double complex top = x[BaseT + k];
            double complex bot = (ccos(2*pi*k/Np) - I*csin(2*pi*k/Np))*x[BaseB + k];
            x[BaseT + k] = top+bot;
            x[BaseB + k] = top-bot;
        }
        BaseT = BaseT + Np;
    }
    Bp =  Bp>>1;
    Np = Np<<1;
}

}

Output is printed with:

 for(int i = 0; i < LENGTH; i++){
    printf("%f + %fj\n", creal(x[i]), cimag(x[i]));
}

Here's a length 8 input, and the correct output:

double complex x[LENGTH] = {1.0, -1.0, 1.0, -1.0, 10.0, -1.0, 5.0, -1.0};

Output:

13.000000 + 0.000000j
-9.000000 + 4.000000j
5.000000 + 0.000000j
-9.000000 + -4.000000j
21.000000 + 0.000000j
-9.000000 + 4.000000j
5.000000 + 0.000000j
-9.000000 + -4.000000j

Here's a length 16 input, the incorrect output that I get, and what it should be:

double complex x[SIZE] = {10.0, -1.0, 5.0, -1.0, 4.0, 4.0, 2.0, 6.0, 9.0, 3.0, 8.0, 4.0, 5.0, 4.0, 3.0, 8.0};

Output (incorrect):

73.000000 + 0.000000j
-4.882497 + 10.451032j
12.535534 + 5.020815j
6.975351 + 7.165394j
12.000000 + 7.000000j
-0.732710 + 0.094326j
5.464466 + 19.020815j
2.639856 + 3.379964j
19.000000 + 0.000000j
2.639856 + -3.379964j
5.464466 + -19.020815j
-0.732710 + -0.094326j
12.000000 + -7.000000j
6.975351 + -7.165394j
12.535534 + -5.020815j
-4.882497 + -10.451032j

Correct output:

73.000000 + 0.000000j
-4.175390 + 10.743925j
13.535534 + 4.020815j
6.268244 + 5.458287j
10.000000 + 7.000000j
-1.439817 + 1.801433j
6.464466 + 20.020815j
3.346963 + 3.087071j
19.000000 + 0.000000j
3.346963 + -3.087071j
6.464466 + -20.020815j
-1.439817 + -1.801433j
10.000000 + -7.000000j
6.268244 + -5.458287j
13.535534 + -4.020815j
-4.175390 + -10.743925j

For this length 16 example, it gives the correct output (!?):

double complex x[SIZE] = {30.0, -1.0, 4.0, -6.0, 4.0, 9.0, 2.0, 6.0, 9.0, 3.0, -8.0, 4.0, -5.0, 4.0, 3.0, 8.0};

66.000000 + 0.000000j
22.604378 + -9.862676j
43.535534 + 28.091883j
24.900438 + 4.851685j
37.000000 + -3.000000j
-1.285214 + 2.408035j
36.464466 + 10.091883j
37.780399 + 23.693673j
12.000000 + 0.000000j
37.780399 + -23.693673j
36.464466 + -10.091883j
-1.285214 + -2.408035j
37.000000 + 3.000000j
24.900438 + -4.851685j
43.535534 + -28.091883j
22.604378 + 9.862676j

EDIT: I've realised my problem: my bit reversing loop was completely wrong. I should have thought about that seemingly simple part and tested it more thoroughly. The reason it worked for the last case was because it had repeated values, which coincidentally gave the same vector as the working reversal. Here's the fixed loop:

int s = log2(N);

    for(int i = 0; i < N; i++){
        int h = bitrev(i, s);
        if(i < h){
        double complex temp = x[i];
        x[i] = x[h];
        x[h] = temp;
        }
    }
6
  • 2
    I'm not familiar with FFT code. Can you provide a sample data set of size 16 or 32, along with the output you get (and the code that prints that answer), and indicate what answer you think you should get. It might help. Sep 29, 2014 at 15:17
  • @JonathanLeffler I've added the requested items Sep 29, 2014 at 15:27
  • wild guess: int overflow? Sep 29, 2014 at 15:28
  • Largest integer used is the size of the array (8, 16, 32). The doubles used are within ~2 degrees of magnitude too. Using long doubles didn't help improve anything. Sep 29, 2014 at 15:33
  • Turn on all compiler warnings and fix all those issues first. That should be your first step. After that, I don't know. Maybe your arithmetic order isn't particularly floating-point friendly.
    – indiv
    Sep 29, 2014 at 15:42

1 Answer 1

1

This part 2*pi*k/Np is letting the compiler decide type casting. I would also check these values are being computed properly. For example, I get different results when using 2*pi*k/Np; and 2.0*pi*(double) k*(double) Np;

(Mods: I rather add this as comment but I don't have enough reputation. Feel free to move my response :) )

1
  • Thanks,I've changed to explicit casts. I've realised what my actual problem was: my bit reversal loop was not correct. I've edited with the correct code. Sep 30, 2014 at 1:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.