Программа 2d DCT не работает

Я должен сделать это 2d DCT изображения для моего проекта. Я перевел формулу прямо в код. Логически все выглядит хорошо, но не дает нужного результата. Я сопоставил это с функцией Matlab, чтобы проверить результаты для матрицы 3x3, но они неверны.

Кроме того, то, что и как я закодировал, дает огромное количество циклов, так что фактическая операция с изображением выполняется часами.

Любое предложение по уменьшению количества циклов и указаний на программную ошибку было бы замечательным. Спасибо.

Это мой код.

    double alpha_p, alpha_q;
    double pi = Math.atan(1.0) * 4.0;
    //dct begins
    System.out.println("it begins");
    for (int p = 0; p < M; p++) {
        for (int q = 0; q < N; q++) {
            if (p == 0)
                alpha_p = 1 / sqrt(M);
            else
                alpha_p = sqrt(2 / M);
            if (q == 0)
                alpha_q = 1 / sqrt(N);
            else
                alpha_q = sqrt(2 / N);
            double toreturn = 0;
            for (int m = 0; m < M; m++) {
                for (int n = 0; n < N; n++) {
                    toreturn = toreturn + img[m][n]
                            * cos(((2 * m + 1) * p * pi) / 2 * M)
                            * cos(((2 * n + 1) * q * pi) / 2 * N);
                }
            }
            dctimg[p][q] = alpha_p * alpha_q * toreturn;
            System.out.println("euta");
        }
    }
    // dct over
    System.out.println("its over");

    //inverse dct begins
    for (int m = 0; m < M; m++) {
        for (int n = 0; n < N; n++) {
            double toreturn = 0;
            for (int p = 0; p < M; p++) {
                for (int q = 0; q < N; q++) {
                    if (p == 0)
                        alpha_p = 1 / sqrt(M);
                    else
                        alpha_p = sqrt(2 / M);
                    if (q == 0)
                        alpha_q = 1 / sqrt(N);
                    else
                        alpha_q = sqrt(2 / N);
                    toreturn = toreturn + alpha_p * alpha_q * dctimg[p][q]
                                          * cos(((2 * m + 1) * p * pi) / 2 * M)
                                          * cos(((2 * n + 1) * q * pi) / 2 * N);
                }
            }
            finalimg[m][n] = toreturn;
        }
    }
    //inverse dct over