Использование fftw_plan_many_dft вместе с броненосной матрицей приводит к неверным данным для FFT, но корректным данным для IFFT - PullRequest
Новая
       104

Использование fftw_plan_many_dft вместе с броненосной матрицей приводит к неверным данным для FFT, но корректным данным для IFFT

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

Я хотел бы применить n независимых FFT к сложной 2d-матрице, созданной броненосцем, то есть по одному FFT для каждой строки. Я придерживаюсь предложенного здесь подхода: Взаимодействуя с броненосцем 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;

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;
}

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;
}

arma::cx_mat multi_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);
    int NCOL=in_mat.n_cols, NROW=in_mat.n_rows;
    int frank=1;
    int howmany=NCOL;
    int n[]={NROW};
    int idist=NROW, odist=NROW;
    int istride=1, ostride=1;
    int *inembed=n, *onembed=n;
    p = fftw_plan_many_dft (frank,
                            n,
                            howmany,
                            (fftw_complex*)in_mat.memptr(),
                            inembed,
                            istride,
                            idist,
                            (fftw_complex*)ret_mat.memptr(),
                            onembed,
                            ostride,
                            odist,
                            FFTW_FORWARD,
                            FFTW_ESTIMATE);
    fftw_execute (p);
    fftw_destroy_plan (p);
    return ret_mat;
}

arma::cx_mat multi_ifftw_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);
    int NCOL=in_mat.n_cols, NROW=in_mat.n_rows;
    int frank=1;
    int howmany=NCOL;
    int n[]={NROW};
    int idist=NROW, odist=NROW;
    int istride=1, ostride=1;
    int *inembed=n, *onembed=n;
    p = fftw_plan_many_dft (frank,
                            n,
                            howmany,
                            (fftw_complex*)in_mat.memptr(),
                            inembed,
                            istride,
                            idist,
                            (fftw_complex*)ret_mat.memptr(),
                            onembed,
                            ostride,
                            odist,
                            FFTW_BACKWARD,
                            FFTW_ESTIMATE);
    fftw_execute (p);
    fftw_destroy_plan (p);
    return ret_mat / matrix_size;
}

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 fftw_multi_mat = multi_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);
    arma::cx_mat ifftw_multi_mat = multi_ifftw_fft (fftw_multi_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';
    std::cout << "Multi-Mat-output is comparable for forward-operation: " << arma::approx_equal(fftw_test_mat, fftw_multi_mat, "rel_tol", 0.01)
              << " and backward operation: " << arma::approx_equal(ifftw_test_mat, ifftw_multi_mat, "rel_tol", 0.01)
              << '\n';
}

Мой вывод: 

Input and output for arma are equal: 1
Forward output for arma and fftw are equal: 1
Output for arma and fftw are equal: 1
Multi-Mat-output is comparable for forward-operation: 0 and backward operation: 1

т.е. Функция dft делает что-то еще с моими данными по сравнению с другими функциями, где я перебираю каждую строку. Тем не менее, при преобразовании обратно я получаю правильный результат.
Я также проверил, перебирал ли я столбцы вместо строк, и транспонировал входную матрицу, но это тоже не помогло. Есть что-то еще, что я пропустил?

...