Выполнение FFT-вычислений с FFTW_MEASURE дает неверные результаты по сравнению с FFTW_ESTIMATE - PullRequest
Новая
       96

Выполнение FFT-вычислений с FFTW_MEASURE дает неверные результаты по сравнению с FFTW_ESTIMATE

0 голосов
/ 14 апреля 2020

Я хотел бы применить 1d-FFT к каждой строке моей сложной матрицы, созданной в Armadillo. Для этого я могу использовать встроенную функцию или использовать FFTW. Для этого я написал короткую тестовую программу: 

#include <iostream>
#include <arrayfire.h>
#include <armadillo>
#include <vector>
#include <chrono>
#include <complex>
#include <fftw3.h>

constexpr int matrix_size = 1024;
constexpr int bench_rounds = 10;

#define USE_ESTIMATE

arma::cx_mat arma_fft(const arma::cx_mat &in_mat){
    arma::cx_mat ret_mat(in_mat.n_cols, in_mat.n_rows, arma::fill::zeros);
    for(size_t i = 0; i < in_mat.n_rows; ++i)
        ret_mat.row(i) = arma::fft(in_mat.row(i));
    return ret_mat;
}

arma::cx_mat arma_ifft(const arma::cx_mat &in_mat){
    arma::cx_mat ret_mat(in_mat.n_cols, in_mat.n_rows, arma::fill::zeros);
    for(size_t i = 0; i < in_mat.n_rows; ++i)
        ret_mat.row(i) = arma::ifft(in_mat.row(i));
    return ret_mat;
}

#ifndef USE_ESTIMATE
arma::cx_mat fftw_fft(const arma::cx_mat &in_mat){
    fftw_plan p;
    arma::cx_mat ret_mat(in_mat.n_cols, in_mat.n_rows, arma::fill::zeros);
    arma::cx_rowvec temp_vec = in_mat.row(0);
    arma::cx_rowvec out_vec(temp_vec.size(), arma::fill::zeros);
    p = fftw_plan_dft_1d (matrix_size,
                          (fftw_complex*)temp_vec.memptr(),
                          (fftw_complex*)out_vec.memptr(),
                          FFTW_FORWARD,
                          FFTW_MEASURE);
    for(size_t i = 0; i < in_mat.n_rows; ++i){
        temp_vec = in_mat.row(i);
        out_vec = arma::cx_rowvec(temp_vec.size(), arma::fill::zeros);

        fftw_execute (p);

        ret_mat.row(i) = out_vec;
    }
    fftw_destroy_plan (p);
    return ret_mat;
}

arma::cx_mat fftw_ifft(const arma::cx_mat &in_mat){
    fftw_plan p;
    arma::cx_mat ret_mat(in_mat.n_cols, in_mat.n_rows, arma::fill::zeros);
    arma::cx_rowvec temp_vec = in_mat.row(0);
    arma::cx_rowvec out_vec(temp_vec.size(), arma::fill::zeros);
    p = fftw_plan_dft_1d (matrix_size,
                          (fftw_complex*)temp_vec.memptr(),
                          (fftw_complex*)out_vec.memptr(),
                          FFTW_BACKWARD,
                          FFTW_MEASURE);
    for(size_t i = 0; i < in_mat.n_rows; ++i){
        temp_vec = in_mat.row(i);
        out_vec = arma::cx_rowvec(temp_vec.size(), arma::fill::zeros);

        fftw_execute (p);

        ret_mat.row(i) = out_vec;
    }
    fftw_destroy_plan (p);
    return ret_mat / matrix_size;
}
#else
arma::cx_mat fftw_fft(const arma::cx_mat &in_mat){
    fftw_plan p;
    arma::cx_mat ret_mat(in_mat.n_cols, in_mat.n_rows, arma::fill::zeros);
    arma::cx_rowvec temp_vec = in_mat.row(0);
    arma::cx_rowvec out_vec(temp_vec.size(), arma::fill::zeros);

    for(size_t i = 0; i < in_mat.n_rows; ++i){
        temp_vec = in_mat.row(i);
        out_vec = arma::cx_rowvec(temp_vec.size(), arma::fill::zeros);

        p = fftw_plan_dft_1d (matrix_size,
                              (fftw_complex*)temp_vec.memptr(),
                              (fftw_complex*)out_vec.memptr(),
                              FFTW_FORWARD,
                              FFTW_ESTIMATE);

        fftw_execute (p);

        ret_mat.row(i) = out_vec;
    }
    fftw_destroy_plan (p);
    return ret_mat;
}

arma::cx_mat fftw_ifft(const arma::cx_mat &in_mat){
    fftw_plan p;
    arma::cx_mat ret_mat(in_mat.n_cols, in_mat.n_rows, arma::fill::zeros);
    arma::cx_rowvec temp_vec = in_mat.row(0);
    for(size_t i = 0; i < in_mat.n_rows; ++i){
        temp_vec = in_mat.row(i);
        arma::cx_rowvec out_vec = arma::cx_rowvec(temp_vec.size(), arma::fill::zeros);
        p = fftw_plan_dft_1d (matrix_size,
                              (fftw_complex*)temp_vec.memptr(),
                              (fftw_complex*)out_vec.memptr(),
                              FFTW_BACKWARD,
                              FFTW_ESTIMATE);
        fftw_execute (p);

        ret_mat.row(i) = out_vec;
    }
    fftw_destroy_plan (p);
    return ret_mat / matrix_size;
}
#endif
int main()
{
    //Check correctness
    arma::cx_mat test_mat = arma::randu<arma::cx_mat>(matrix_size, matrix_size);
    arma::cx_mat fft_test_mat = arma_fft(test_mat);
    arma::cx_mat fftw_test_mat = fftw_fft (test_mat);
    arma::cx_mat ifft_test_mat = arma_ifft(fft_test_mat);
    arma::cx_mat ifftw_test_mat = fftw_ifft (fftw_test_mat);
    std::cout << "Input and output for arma are equal: " << arma::approx_equal(test_mat, ifft_test_mat, "rel_tol", 0.01) << '\n';
    std::cout << "Forward output for arma and fftw are equal: " << arma::approx_equal(fft_test_mat, fftw_test_mat, "rel_tol", 0.01) << '\n';
    std::cout << "Output for arma and fftw are equal: " << arma::approx_equal(ifftw_test_mat, ifft_test_mat, "rel_tol", 0.01) << '\n';
}

При использовании FFTW_ESTIMATE мои результаты равны. При использовании FFTW_MEASURE для более высокой скорости обработки мои результаты неверны, и примерно половина полученной fftw-матрицы равна нулю. Почему? Я что-то здесь забыл?

...