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

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

Учитывая простую (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

`

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

В: (+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): "Во многих практических приложениях входные данные в [i] - чисто действительные числа, и в этом случае выход DFT удовлетворяет "эрмитовой" избыточности: out [i] является сопряжением out [ni]. Этими обстоятельствами можно воспользоваться, чтобы примерно в два раза улучшить как скорость, так и использование памяти ".

Таким образом, матрица 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 ответа

Решение

У вас есть реальная функция 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')

Но в этих обстоятельствах я знаю, что, например, только из коэффициентов на верхней диагонали я могу восстановить 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 завершено для обратного преобразования в реальное.

Источник

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

#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;
}
Источник
Другие вопросы по тегам