Algorithme FFT récursif en Java renvoyant null ?

Question

Je suis actuellement en train de mettre en œuvre l'algorithme FFT en Java et j'ai un peu de mal avec ça! J'ai testé toutes les autres parties de l'algorithme et elles semblent fonctionner correctement.

Le problème que je rencontre est que dans le cas de base, un tableau de nombres complexes est retourné, dans le cas de base A[0] est renseigné. Après l'exécution des cas de base, la boucle for est exécutée où y0[0] et y1[0] sont trouvés pour être null, malgré leur attribution dans les cas de base, assez déconcerté par cela. Cela est visible dans la ligne System.out.println

Quelqu'un peut-il me dire où se trouvent mes erreurs?

    //Cette méthode implémente l'algorithme FFT récursif, elle suppose que la longueur d'entrée
//N est une puissance de deux
private static Complex[] FFT(Complex[] A, int N) { 
    double real, imag;
    Complex A0[] = new Complex[((int) Math.ceil(N/2))]; 
    Complex A1[] = new Complex[((int) Math.ceil(N/2))];
    Complex[] omega = new Complex[N];
    Complex[] y = new Complex[N];
    Complex[] y0 = new Complex[((int) Math.ceil(N/2))];
    Complex[] y1 = new Complex[((int) Math.ceil(N/2))]; 

    //cas de base
    if (N == 1) {
        return A;
    }
    else {
        real = Math.cos((2*Math.PI)/N); if (real < 1E-10 && real > 0) real = 0;
        imag = Math.sin((2*Math.PI)/N); if (imag < 1E-10 && imag > 0) imag = 0;
        omega[N-1] = new Complex(real, imag);
        omega[0] = new Complex(1, 0);
        A0 = splitInput(A, 1);
        A1 = splitInput(A, 0);
        //appels récursifs
        y0 = FFT(A0, N/2);
        y1 = FFT(A1, N/2);
        for (int k = 0; k < ((N/2)-1); k++) {
            System.out.print("k: " + k + ", y0: " + y0[k]); System.out.println(", y1: " + y1[k]);
            y[k] = y0[k].plus(omega[k].times(y1[k]));
            y[k+(N/2)] = y0[k].minus(omega[k].times(y1[k]));
            omega[0] = omega[0].times(omega[N]);
        }
        return y;
    }
}

Voici le code de ma méthode splitInput comme demandé

//Cette méthode prend un tableau de nombres complexes en argument et retourne chaque élément pair ou impair
//selon le deuxième argument int étant 1 ou 0
private static Complex[] splitInput(Complex[] input, int even) {
    Complex[] newArray = new Complex[(input.length/2)];

    //Renvoie tous les éléments pairs du tableau double, y compris 0
    if (even == 1) {
        for (int i = 0; i < (input.length/2); i++) {
            newArray[i] = new Complex(input[i*2].re, 0.0);
        }
        return newArray;
    }
    //Renvoie tous les éléments impairs du tableau double
    else {
        for (int i = 0; i < (input.length/2); i++) {
            newArray[i] = new Complex (input[(i*2) + 1].re, 0.0);
        }
    return newArray;
    }
}

MODIFICATION: J'ai mis à jour mon code selon vos suggestions, je reçois toujours une exception de pointeur nul à la ligne y[k] = y0[k].plus(omega[k].times(y1[k])); alors que y0 et y1 sont toujours null après le cas de base :( des idées supplémentaires? Voici l'algorithme mis à jour

//Cette méthode implémente l'algorithme FFT récursif, elle suppose que la longueur d'entrée
//N est une puissance de deux
private static Complex[] FFT(Complex[] A, int N) { 
    double real, imag;
    Complex[] omega = new Complex[N];
    Complex[] y = new Complex[N];
    Complex[] A0;
    Complex[] A1;
    Complex[] y0;
    Complex[] y1;

    //cas de base
    if (N == 1) {
        return A;
    }
    else {
        real = Math.cos((2*Math.PI)/N); if (real < 1E-10 && real > 0) real = 0;
        imag = Math.sin((2*Math.PI)/N); if (imag < 1E-10 && imag > 0) imag = 0;;
        omega[N-1] = new Complex(real, imag);
        omega[0] = new Complex(1, 0);
        A0 = splitInput(A, 1);
        A1 = splitInput(A, 0);
        //appels récursifs
        y0 = FFT(A0, N/2);
        y1 = FFT(A1, N/2);
        for (int k = 0; k < ((N/2)-1); k++) {
            y[k] = y0[k].plus(omega[k].times(y1[k]));
            y[k+(N/2)] = y0[k].minus(omega[k].times(y1[k]));
            omega[0] = omega[0].times(omega[N-1]);
        }
        return y;
    }
}

Answer 1

Quelques pensées :

Chaque fois que je vois quelque chose répété aussi souvent que Math.ceil(N/2) ici, je pense que cela justifie d'avoir sa propre variable nommée. (Je sais que nommer les variables n'est pas toujours facile, mais je trouve que c'est essentiel pour la lisibilité.)

Complex A0[] = new Complex[((int) Math.ceil(N/2))];
Complex A1[] = new Complex[((int) Math.ceil(N/2))];

Remarquez que lorsque N == 1, le calcul aboutit à new Complex[0]. Je ne suis pas sûr de ce que cela fait, mais je pense que je placerais la vérification du cas de base N == 1 avant les allocations de mémoire.

Complex[] y0 = new Complex[((int) Math.ceil(N/2))];
Complex[] y1 = new Complex[((int) Math.ceil(N/2))];
/* ... */
    y0 = FFT(A0, N/2);
    y1 = FFT(A1, N/2);

Je pense que vous pouvez sauter les allocations new Complex[...] pour ces tableaux car vous ne stockez jamais rien dedans en fait.

Complex[] omega = new Complex[N];
/* ... */
        omega[0] = omega[0].times(omega[N]);

Je suis surpris que cela n'ait pas encore explosé -- omega[N] devrait déclencher une exception IndexOutOfBounds.

Answer 2

Problèmes qui sautent aux yeux :

1. (int) Math.ceil(N/2) Vous continuez à effectuer une division int, donc le Math.ceil() n'a aucun effet et vos tableaux divisés sont probablement incorrects pour un n impair
2. vous ne remplissez que omega[0] et omega[N-1], je m'attendrais à une NullPointerException lorsque vous essayez d'accéder à omega[1], ce qui se produirait lorsque N >= 6.
3. omega[N], comme mentionné également par sarnold
4. Vous allouez A0 et A1, et leur attribuez plus tard les résultats de splitInput