Использование FFTW в MATLAB против C/C++: разные результаты ⇐ C++

Программы на C++. Форум разработчиков
Ответить Пред. темаСлед. тема
Гость
 Использование FFTW в MATLAB против C/C++: разные результаты

Сообщение Гость »


Я изучаю различия между fftw в MATLAB и C/C++. Я думаю, что это может быть простой вопрос, но по какой-то причине я его не вижу.
Предположим следующий код MATLAB:

Код: Выделить всё

Nx = 8;
Ny = 8;
Lx =16;
%-----------
xi_x =  (2*pi)/Lx;
xi  = ((0:Nx-1)/Nx)*(2*pi);
x =  xi/xi_x;
ylow = 0;
yupp =6;
ygl  = cos ( pi * ( 0 : Ny ) / Ny )';
ygl = (1/2)*(((yupp-ylow)*ygl) + (yupp+ylow));
[X,Y]   = meshgrid(x,ygl);
%define function u:
u = (Y-ylow) .* (Y-yupp) .* sin( (2*pi / Lx) * X);
uh = fft(u,[],2);
%pad uh with zeros
uhp=[uh(:,1:Nx/2) zeros(Ny,(m-Nx)) uh(:,Nx/2+1:Nx)];
%take ifft
up=real(ifft(uhp,[],2));
u1 = real(ifft(uh,[],2));  % то же, что и вы
Мой вопрос выше, поскольку оба

Код: Выделить всё

uh
and

Код: Выделить всё

uhp
are the same values with extra zero padding, why is the inverse of fft of these two different? Meaning if

Код: Выделить всё

uh
is:

Код: Выделить всё

0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
0.0000 + 0.0000i   0.0000 + 5.2721i  -0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.0000 - 0.0000i  -0.0000 - 0.0000i   0.0000 - 5.2721i
-0.0000 + 0.0000i   0.0000 +18.0000i  -0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.0000 - 0.0000i  -0.0000 - 0.0000i   0.0000 -18.0000i
-0.0000 + 0.0000i   0.0000 +30.7279i  -0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.0000 - 0.0000i  -0.0000 - 0.0000i   0.0000 -30.7279i
-0.0000 + 0.0000i   0.0000 +36.0000i  -0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.0000 - 0.0000i  -0.0000 - 0.0000i   0.0000 -36.0000i
-0.0000 + 0.0000i   0.0000 +30.7279i  -0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.0000 - 0.0000i  -0.0000 - 0.0000i   0.0000 -30.7279i
-0.0000 + 0.0000i   0.0000 +18.0000i  -0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.0000 - 0.0000i  -0.0000 - 0.0000i   0.0000 -18.0000i
0.0000 + 0.0000i   0.0000 + 5.2721i  -0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.0000 - 0.0000i  -0.0000 - 0.0000i   0.0000 - 5.2721i
0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
and

Код: Выделить всё

uhp
is:

Код: Выделить всё

0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
0.0000 + 0.0000i   0.0000 + 5.2721i  -0.0000 + 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.0000 - 0.0000i  -0.0000 - 0.0000i   0.0000 - 5.2721i
-0.0000 + 0.0000i   0.0000 +18.0000i  -0.0000 + 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.0000 - 0.0000i  -0.0000 - 0.0000i   0.0000 -18.0000i
-0.0000 + 0.0000i   0.0000 +30.7279i  -0.0000 + 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.0000 - 0.0000i  -0.0000 - 0.0000i   0.0000 -30.7279i
-0.0000 + 0.0000i   0.0000 +36.0000i  -0.0000 + 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.0000 - 0.0000i  -0.0000 - 0.0000i   0.0000 -36.0000i
-0.0000 + 0.0000i   0.0000 +30.7279i  -0.0000 + 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.0000 - 0.0000i  -0.0000 - 0.0000i   0.0000 -30.7279i
-0.0000 + 0.0000i   0.0000 +18.0000i  -0.0000 + 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.0000 - 0.0000i  -0.0000 - 0.0000i   0.0000 -18.0000i
0.0000 + 0.0000i   0.0000 + 5.2721i  -0.0000 + 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.0000 - 0.0000i  -0.0000 - 0.0000i   0.0000 - 5.2721i
0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
Then why the

Код: Выделить всё

 real(ifft(uh,[],2))
equals:

Код: Выделить всё

    0         0         0         0         0         0         0         0
0   -0.9320   -1.3180   -0.9320   -0.0000    0.9320    1.3180    0.9320
-0.0000   -3.1820   -4.5000   -3.1820   -0.0000    3.1820    4.5000    3.1820
0   -5.4320   -7.6820   -5.4320   -0.0000    5.4320    7.6820    5.4320
-0.0000   -6.3640   -9.0000   -6.3640   -0.0000    6.3640    9.0000    6.3640
0   -5.4320   -7.6820   -5.4320   -0.0000    5.4320    7.6820    5.4320
-0.0000   -3.1820   -4.5000   -3.1820   -0.0000    3.1820    4.5000    3.1820
0   -0.9320   -1.3180   -0.9320   -0.0000    0.9320    1.3180    0.9320
0         0         0         0         0         0         0         0
and

Код: Выделить всё

up=ifft(uhp,[],2);
:

Код: Выделить всё

 0         0         0         0         0         0         0         0         0         0         0         0
0   -0.4393   -0.7610   -0.8787   -0.7610   -0.4393   -0.0000    0.4393    0.7610    0.8787    0.7610    0.4393
-0.0000   -1.5000   -2.5981   -3.0000   -2.5981   -1.5000   -0.0000    1.5000    2.5981    3.0000    2.5981    1.5000
0   -2.5607   -4.4352   -5.1213   -4.4352   -2.5607   -0.0000    2.5607    4.4352    5.1213    4.4352    2.5607
-0.0000   -3.0000   -5.1962   -6.0000   -5.1962   -3.0000   -0.0000    3.0000    5.1962    6.0000    5.1962    3.0000
0   -2.5607   -4.4352   -5.1213   -4.4352   -2.5607   -0.0000    2.5607    4.4352    5.1213    4.4352    2.5607
-0.0000   -1.5000   -2.5981   -3.0000   -2.5981   -1.5000   -0.0000    1.5000    2.5981    3.0000    2.5981    1.5000
0   -0.4393   -0.7610   -0.8787   -0.7610   -0.4393   -0.0000    0.4393    0.7610    0.8787    0.7610    0.4393
0         0         0         0         0         0         0         0         0         0         0         0
The reason I ask b/c in C/C++ when padding my arrays similarly and taking the ifft my code returns the same output, so

Код: Выделить всё

real(ifft(uh,[],2))
=

Код: Выделить всё

real(ifft(uhp,[],2))
!
Here's my C/C++ code:

Код: Выделить всё

static const int nx = 8;
static const int ny = 8;
static const int ncomp = 2;
static const int nxk = nx/2 + 1;
static const int Mx = 3*nx/2, My= 3*ny/2;
double makeSpatialMesh2DTrans(double dx, double dy, double ylow, double yupp,double xarr[], double yarr[]);
int xy(int x, int y, int number_of_cols);
class MyFftwClass {
public:
MyFftwClass(void);
~MyFftwClass(void);
void r2cexecute(double rArr[], double cArr[][ncomp]);

private:
static fftw_plan s_plan;
double *m_buffer_in;
fftw_complex *m_buffer_out;

};
class MyiFftwClass { //inverse
public:
MyiFftwClass(void);
~MyiFftwClass(void);
void c2rexecute(double cArr[][ncomp],double rArr[]);

private:
static fftw_plan s_plan1;
fftw_complex *m_buffer_in1;
double *m_buffer_out1;
};
MyFftwClass r2c;
MyiFftwClass c2r;
int main(){
double Lx = 16.;
double ylow = 0.;
double yupp = 6.;
double Ly = (yupp-ylow);
double dx = Lx / nx;
double dy = Ly / ny;
double *XX;
XX = (double*) fftw_malloc((nx*(ny+1))*sizeof(double));
memset(XX, 42, (nx*(ny+1)) * sizeof(double));
double *YY;
YY = (double*) fftw_malloc((nx*(ny+1))*sizeof(double));
memset(YY, 42, (nx*(ny+1)) * sizeof(double));

double kmax = makeSpatialMesh2DTrans(dx, dy, ylow, yupp, XX, YY);
double *uFun;
uFun = (double*) fftw_malloc((nx*(ny+1))*sizeof(double));
double A = 2*M_PI/Lx;
for(int y = 0; y< ny+1; y++){
for(int x = 0; x< nx; x++){
uFun[xy(x, y, nx)] = ( (YY[xy(x, y, nx)]-ylow) * (YY[xy(x, y, nx)]-yupp) * sin(A*XX[xy(x, y, nx)]));
}
}
}

fftw_complex *uFunk;
uFunk = (fftw_complex*) fftw_malloc((nx*(ny + 1)) * sizeof(fftw_complex));
memset(uFunk, 42, (nx*(ny + 1))* sizeof(fftw_complex));
r2c.r2cexecute(uFun,uFunk);

double N[2] = {nx,ny};
double M[2] = {Mx,My};

fftw_complex *uFunk_pad;
uFunk_pad = (fftw_complex*) fftw_malloc((Mx*My)* sizeof(fftw_complex));

//Add zero padding to the function u
std::cout 

Источник: [url]https://stackoverflow.com/questions/78144549/taking-fftw-in-matlab-vs-c-c-different-outputs[/url]
Вернуться к началу
Реклама
Ответить Пред. темаСлед. тема

Быстрый ответ

Изменение регистра текста: 
Смайлики
:) :( :oops: :roll: :wink: :muza: :clever: :sorry: :angel: :read: *x)
Ещё смайлики…
   
К этому ответу прикреплено по крайней мере одно вложение.

Если вы не хотите добавлять вложения, оставьте поля пустыми.

Максимально разрешённый размер вложения: 15 МБ.

  • Похожие темы
    Ответы
    Просмотры
    Последнее сообщение

Вернуться в «C++»