Armadillo 2D FFTW от реальной до сложной DFT симметрии Герметия - PullRequest
Новая
       111

Armadillo 2D FFTW от реальной до сложной DFT симметрии Герметия

0 голосов
/ 02 мая 2018

Я пытаюсь выполнить 2-мерное преобразование Фурье некоторых реальных данных. Для этого я использую библиотеку FFTW, так как она значительно быстрее библиотеки Армадилло.

Учитывая простую (4x4) стартовую матрицу: AAA: 

    0        0        0        0
    0   1.0000   2.0000   3.0000
    0   2.0000   4.0000   6.0000
    0   3.0000   6.0000   9.0000

`

Если я использую встроенный БПФ в броненосце, вывод будет выглядеть следующим образом:

BBB: (+3.600e+01,+0.000e+00) (-1.200e+01,+1.200e+01) (-1.200e+01,+0.000e+00) (-1.200e+01,-1.200e+01) (-1.200e+01,+1.200e+01) (+0.000e+00,-8.000e+00) (+4.000e+00,-4.000e+00) (+8.000e+00,+0.000e+00) (-1.200e+01,+0.000e+00) (+4.000e+00,-4.000e+00) (+4.000e+00,+0.000e+00) (+4.000e+00,+4.000e+00) (-1.200e+01,-1.200e+01) (+8.000e+00,+0.000e+00) (+4.000e+00,+4.000e+00) (+0.000e+00,+8.000e+00)

Но если я использую FFTW, я получаю:

CCC: (+3.600e+01,+0.000e+00) (+0.000e+00,-8.000e+00) (+4.000e+00,+0.000e+00) (0,0) (-1.200e+01,+1.200e+01) (+4.000e+00,-4.000e+00) (-1.200e+01,-1.200e+01) (0,0) (-1.200e+01,+0.000e+00) (-1.200e+01,+0.000e+00) (+8.000e+00,+0.000e+00) (0,0) (-1.200e+01,+1.200e+01) (+4.000e+00,-4.000e+00) (+4.000e+00,+4.000e+00) (0,0

Выполнение соответствующего IFFT на матрице BBB и CCC дает в точности начальную матрицу AAA.

Согласно документации: (http://www.fftw.org/fftw3_doc/One_002dDimensional-DFTs-of-Real-Data.html#One_002dDimensional-DFTs-of-Real-Data): «Во многих практических применениях входные данные в [i] являются чисто действительными числами, и в этом случае выходной сигнал DFT удовлетворяет« эрмитовой »избыточности: out [i] является сопряжением out [n-i]. Этими обстоятельствами можно воспользоваться, чтобы примерно в два раза улучшить как скорость, так и использование памяти ».

Следовательно, матрица CCC требует какой-то операции для извлечения герметичности Гермета, но я слишком математически нуб, чтобы понять, что это за операция. Может ли кто-нибудь помочь мне с этим?

Кроме того, Armadillo хранит данные в формате col Major и FFTW в формате основных строк, в соответствии с документацией это не должно иметь значения, если вы передаете размеры строки / столбца в обратном порядке в функцию плана?

Спасибо за внимание.

Прилагается мой код: 

#include <iostream>
#include <fftw3.h>
#include "armadillo"

using namespace arma;
using namespace std;


int main(int argc, char** argv)
{

    mat AAA=zeros(4,4);
    mat IBB=zeros(4,4);
    cx_mat BBB(4,4);



    for (int xx=0;xx<=3;xx++){
        for ( int yy=0;yy<=3;yy++){

    AAA(xx,yy)= xx*yy;
        }
    }


    cx_mat CCC (4,4);
    cx_mat CCCC(4,4);
    mat ICC =zeros(4,4);



     fftw_plan plan=fftw_plan_dft_r2c_2d(4, 4,(double(*))&AAA(0,0), (double(*)[2])&CCC(0,0), FFTW_ESTIMATE);
     fftw_plan plan2=fftw_plan_dft_c2r_2d(4, 4,(double(*)[2])&CCCC(0,0), (double(*))&ICC(0,0), FFTW_ESTIMATE);

     //Perform Armadillo FFT (Correct output)
     BBB=fft2(AAA);
     //Perform armadillo IFFT
     IBB=real(ifft2(BBB));




    //Perform FFTW- FFT
    fftw_execute(plan);
    //Allocate fourier array to another array as imput array is destroyed
    CCCC=CCC;

    //Perform FFTW- IFFT on newly allocated array
    fftw_execute(plan2);
    //Must re-normalise the array by the number of elements
    ICC=ICC/(4*4);

    //myst rescale by the number of elements in the array

BBB.print("BBB:");
CCC.print("CCC:");  

IBB.print("IBB:");
ICC.print("ICC:");




    return 0;
}
`

Ответы [ 2 ]

0 голосов
/ 08 мая 2018

Просто чтобы добавить к ответу Каве (принятый ответ), чтобы полностью восстановить избыточность Гермета, необходимо было выполнить ДПФ, как указано в его ответе, а затем выбрать подматрицу, игнорируя нулевые частоты. Выполнение переворота слева направо, переворот вверх и вниз, а затем комплексное сопряжение результирующей матрицы. Надеюсь, что это помогает другим. Вот код 

#include <iostream>
#include <fftw3.h>
#include "armadillo"

using namespace arma;
using namespace std;


int main(int argc, char** argv)
{

    mat AAA=zeros(6,6);
    mat IBB=zeros(6,6);
    cx_mat ccgin(6,6);
    cx_mat ccgout(6,6);
    cx_mat BBB(6,6);



    for (int xx=0;xx<=5;xx++){
        for ( int yy=0;yy<=5;yy++){

    AAA(xx,yy)= xx*(xx+yy);
        }
    }


    cx_mat CCC (4,6);
    cx_mat CCCC(4,6);
    mat ICC =zeros(6,6);

    cx_mat con(3,5);



     fftw_plan plan=fftw_plan_dft_r2c_2d(6, 6,(double(*))&AAA(0,0), (double(*)[2])&CCC(0,0), FFTW_ESTIMATE);
     fftw_plan plan2=fftw_plan_dft_c2r_2d(6, 6,(double(*)[2])&CCCC(0,0), (double(*))&ICC(0,0), FFTW_ESTIMATE);

     //Perform Armadillo FFT (Correct output)
     BBB=fft2(AAA);

     //Perform armadillo IFFT
     IBB=real(ifft2(BBB));




    //Perform FFTW- FFT
    fftw_execute(plan);
    //Allocate fourier array to another array as imput array is destroyed
    CCCC=CCC;

    //Perform FFTW- IFFT on newly allocated array
    fftw_execute(plan2);
    //Must re-normalise the array by the number of elements
    ICC=ICC/(6*6);

    //must rescale by the number of elements in the array

BBB.print("BBB:");
CCC.print("CCC:");  

IBB.print("IBB:");
ICC.print("ICC:");


//Recover Hermetian redundancy
con=fliplr(flipud(conj(CCC(span(1,3),span(1,5)))));
con.print("fliplr(flipud(conj(CCC(span(1,3),span(1,5)))));:");


    return 0;
}
0 голосов
/ 03 мая 2018

У вас есть реальная функция A: 

0        0        0        0
0   1.0000   2.0000   3.0000
0   2.0000   4.0000   6.0000
0   3.0000   6.0000   9.0000

Преобразование Фурье вещественной функции является эрмитовым. Это означает, что действительная часть спектра четна X(iw) = X(-iw), а мнимая часть спектра нечетна X(iw)=-X(-iw). Другими словами 

imag(BBB) == -imag(BBB') && real(BBB) == real(BBB')

Но в этих условиях я знаю, что, например, только из коэффициентов на верхней диагонали я могу восстановить ВВВ для обратного преобразования.

fftw_plan_dft_r2c_2d также объясняет, что именно поэтому CCC должен быть равен nd / 2 + 1 x nd (1 столбец заполнения) для сохранения выходных данных. Таким образом, вы можете объявить CCC и CCCC следующим образом: 

cx_mat CCC (4,3);
cx_mat CCCC(4,3);

Но так как вы упоминаете, что FFTW отличается от майора строк Armadillo, вы должны также отразить последствия этого в своем коде: 

cx_mat CCC (3,4);
cx_mat CCCC(3,4);

И вдруг ваш результат выглядит совершенно иначе: 

BBB:
(+3.600e+01,+0.000e+00) (-1.200e+01,+1.200e+01) (-1.200e+01,+0.000e+00) (-1.200e+01,-1.200e+01)
(-1.200e+01,+1.200e+01) (+0.000e+00,-8.000e+00) (+4.000e+00,-4.000e+00) (+8.000e+00,+0.000e+00)
(-1.200e+01,+0.000e+00) (+4.000e+00,-4.000e+00) (+4.000e+00,+0.000e+00) (+4.000e+00,+4.000e+00)
(-1.200e+01,-1.200e+01) (+8.000e+00,+0.000e+00) (+4.000e+00,+4.000e+00) (+0.000e+00,+8.000e+00)
CCC:
(+3.600e+01,+0.000e+00) (-1.200e+01,+1.200e+01) (-1.200e+01,+0.000e+00) (-1.200e+01,-1.200e+01)
(-1.200e+01,+1.200e+01) (+0.000e+00,-8.000e+00) (+4.000e+00,-4.000e+00) (+8.000e+00,+0.000e+00)
(-1.200e+01,+0.000e+00) (+4.000e+00,-4.000e+00) (+4.000e+00,+0.000e+00) (+4.000e+00,+4.000e+00)
IBB:
    0        0        0        0
    0   1.0000   2.0000   3.0000
    0   2.0000   4.0000   6.0000
    0   3.0000   6.0000   9.0000
ICC:
    0        0        0        0
    0   1.0000   2.0000   3.0000
    0   2.0000   4.0000   6.0000
    0   3.0000   6.0000   9.0000

Из того, что находится в CCC, вы можете восстановить оставшиеся коэффициенты, и ваше обратное преобразование будет правильным. Ударь меня, если что-то останется неясным.

Реконструкция может быть выполнена, например, следующим образом: 

(+3.6e+01,+0.0e+00) (-1.2e+01,+1.2e+01) (-1.2e+01,+0.0e+00 (-1.2e+01,-1.2e+01)
(-1.2e+01,+1.2e+01) (+0.0e+00,-8.0e+00) (+4.0e+00,-4.0e+00 (+8.0e+00,+0.0e+00)
(-1.2e+01,+0.0e+00) (+4.0e+00,-4.0e+00) (+4.0e+00,+0.0e+00 (+4.0e+00,+4.0e+00)
conj(CCC(1,2))      conj(CCC(4,2))      conj(CCC(4,3))     conj(CCC(2,2))

И CCCC завершено для обратного преобразования в реальное.

...