2016-09-05 3 views
0

Я хочу вычислить преобразование БПФ и обратно, чтобы проверить, работает ли оно одинаково. У меня есть приложение больших 3D матрицы в моем коде, я попытался проверить его с 4 * 4 * 4 матрицы и вот мой кодВычисление 3D FFT и обратного БПФ в C

`

#include <stdio.h> 
#include <stdlib.h> 
#include <complex.h> 
#include <time.h> 
#include <math.h> 
#include <fftw3.h> 


int main() 
{ 
int N = 4; //Dimension of matrix 
unsigned int seed = 1; 
double *in = (double*)malloc(sizeof(double)*N*N*N); 
fftw_complex *out = fftw_malloc(sizeof(fftw_complex)*N*N*N); 
double *out1 = (double*)malloc(sizeof(double)*N*N*N); 

fftw_plan plan_backward; 
fftw_plan plan_forward; 

srand (seed); 

for(int i = 0; i < N; i++) 
{ 
    for(int j = 0; j < N; j++) 
    { 
     for(int k = 0; k < N; k++) 
     { 
      in[i*(N*N) + j*N + k] = rand (); 
     } 
    } 
} 

printf(" Given matrix in\n"); 


for(int i = 0; i < N; i++) 
{ 
    for(int j = 0; j < N; j++) 
    { 
     for(int k = 0; k < N; k++) 
     { 
      printf("%f\n", in[i*(N*N) + j*N + k]); 
     } 
    } 
} 


printf("\n"); 

plan_backward = fftw_plan_dft_r2c_3d (N, N, N, in, out, FFTW_ESTIMATE); 

fftw_execute (plan_backward); 

fftw_destroy_plan (plan_backward); 

printf("out matrix\n"); 

for(int i = 0; i < N; i++) 
{ 
    for(int j = 0; j < N; j++) 
    { 
     for(int k = 0; k < N; k++) 
     { 
      printf("%f\t%f\n", creal(out[i*(N*N) + j*N + k]), cimag(out[i*(N*N) + j*N + k])); 
     } 
    } 
} 

printf("\n"); 

plan_forward = fftw_plan_dft_c2r_3d (N, N, N, out, out1, FFTW_ESTIMATE); 

fftw_execute (plan_forward); 

fftw_destroy_plan (plan_forward); 

printf("out1 matrix\n"); 


for(int i = 0; i < N; i++) 
{ 
    for(int j = 0; j < N; j++) 
    { 
     for(int k = 0; k < N; k++) 
     { 
      printf("%f\n", out1[i*(N*N) + j*N + k]); 
     } 
    } 
} 

fftw_free(in); 
free(out); 
fftw_free(out1); 

return 0; 

}` 

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

ответ

0

Ваш FFT не нормализуется. Между вашим входом и выходом существует постоянный фактор.

Посмотрите here

Эти преобразования ненормализованное, так что R2C следует c2r преобразование (или наоборот) приведет к исходным данным масштабируется по числу элементов-то реальные данные, произведение (логических) измерений реальных данных.

Таким образом, коэффициент должен быть N * N * N. Просто разделите данные на этот коэффициент, и вы должны вернуть те же данные, что и ваш ввод.

+0

Получил это Спасибо. Просто разъяснение мне нужно выполнить эту нормализацию только для данных, которые я инвертирую? скажем, я применяю сначала c2r для сложных данных и возвращаюсь, я делаю r2c/(N * N * N)? – verito

+0

Преобразование Фурье линейно, поэтому вы можете выполнять нормализацию везде, где хотите: 'c2r -> 1/N^3 -> r2c',' 1/N^3 -> c2r -> r2c' и т. Д. –

0

Как описано в руководстве по эксплуатации FFTW3:

Этих преобразования являются ненормированными, так что R2C с последующим c2r преобразование (или наоборот) приведет к исходным данным, масштабированных по количеству элементов реальных данных Т. Е. Произведение (логических) размеров реальных данных.

Количество реальных размеров в вашем случае 4^3 и разделив первый результат 115474520512 на это число возвращает первый вход 115474520512/(4^3) = 1804289383.

+0

Спасибо. У меня есть ответ :) – verito