void FourierTransform::FFT(Mat data, Direction direction)
{
    if (direction == Forward)
    {
        int N = data.total();
        int K = N/2;
        if (N <= 1)
        {
            return;
        }
        Mat even = Mat(1, N / 2, data.type(), Scalar(0,0));
        Mat odd = Mat(1, N / 2, data.type(), Scalar(0,0));
        separate(data, even, odd);

        FFT(even, Forward);
        FFT(odd, Forward);

        for (unsigned long u = 0; u < K; ++u) {
            vector<complex<double>> F_u_and_uK = F_u_and_F_uK_FFT(even, odd, u, K);
            data.at<Vec2d>(u)[0] = F_u_and_FuK[0].real();
            data.at<Vec2d>(u)[1] = F_u_and_F_uK[0].imag();

            data.at<Vec2d>(u+K)[0] = F_u_and_F_uK[1].real();
            data.at<Vec2d>(u+K)[1] = F_u_and_F_uK[1].imag();
        }

    }
    else if (direction == Backward)
    {
        multiply(data, Scalar(1, -1), data);
        FFT(data, Forward);
        multiply(data, Scalar(data.total(), data.total()), data);
    }
}

vector<complex<double>> FourierTransform::calculate_F_u_and_F_uK_FFT(Mat even, Mat odd, int u, int K)
{
    vector<complex<double>> result;
    complex<double> F_even(even.at<Vec2d>(u)[0], even.at<Vec2d>(u)[1]);
    complex<double> F_odd(odd.at<Vec2d>(u)[0], odd.at<Vec2d>(u)[1]);
    complex<double> W(cos( double(( -2 * M_PI * u)) / (2*K)), sin( double(( -2 * M_PI * u)) / (2*K)));
    complex<double> F_u(0, 0);
    complex<double> F_uK(0, 0);
    F_u = F_even + F_odd*W;
    F_uK = F_even - F_odd*W;
    result.push_back(F_u);
    result.push_back(F_uK);
    return result;
}