2011-03-23 2 views
1

fftwpp를 사용하여 내 데이터와 회선 커널을 푸리에 공간으로 변환하고 스칼라 곱과 같이 곱한 다음 다시 실제 공간으로 변환합니다. 프로그램을 처음 실행하면 완전히 0으로 채워지는 배열이 만들어집니다. 다시 실행하면 원하는 결과를 얻을 수 있습니다.fftwpp는 컨볼 루션 후에 0으로 채워진 배열을 만듭니다.

실행하는 동안 wisdom3.txt이 생성됩니다. 삭제하면 프로그램이 0으로 채워진 배열을 다시 만드는 데 많은 시간이 걸립니다.

내 코드가 잘못되었습니다.

// sx, sy and sz are the dimensions of my data 

int szp = sz/2 + 1; 
size_t align = sizeof(Complex); 

// creates arrays to store the data in, the double one is for the real data 
// the Complex one for the fourier data 
array3<double> f(sx, sy, sz, align); 
array3<Complex> g(sx, sy, szp, align); 

// copying data into double array 
for(int k = 0; k < sz; k++) 
    for(int j = 0; j < sy; j++) 
     for(int i = 0; i < sx; i++) 
      f(i, j, k) = data[i + sx * j + sx * sy * k]; 

// transforming data into fourier space 
rcfft3d Forward3(sz, f, g); 
Forward3.fft(f, g); 


// generate the kernel 
array3<double> kernel(sx, sy, sz); 
array3<Complex> kernel2(sx, sy, szp, align); 
// more code to create the kernel left out ... 


// transform the kernel into the fourier space 
rcfft3d ForwardKernel3(sz, kernel, kernel2); 
ForwardKernel3.fft(kernel, kernel2); 


// multiplying data and kernel in fourier space together 
for(int k = 0; k < szp; k++) 
    for(int j = 0; j < sy; j++) 
     for(int i = 0; i < sx; i++) 
      g(i, j, k) = g(i, j, k) * kernel2(i, j, k); 


// transform back to normal space 
crfft3d Backward3(sz, g, f); 
Backward3.fftNormalized(g, f); 


// putting everything in the results array and normalize 
for(int k = 0; k < sz; k++) 
    for(int j = 0; j < sy; j++) 
     for(int i = 0; i < sx; i++) 
      result[i + sx * j + sx * sy * k] = 
       (f(i, j, k) >= thresholdValue ? f(i, j, k) : 0); 

답변

3

FFTW ++는 FFTW에 대한 래퍼입니다. FFTW는 데이터를 효율적으로 처리하기위한 실행 계획을 생성해야합니다. 계획이 생성되면 it can be reused을 사용하여 다양한 데이터 집합을 처리합니다. wisdom3.txt이 표시된 파일은 프로그램이 처음 실행될 때 생성되는 계획에 대한 정보입니다. 일단 존재하면 연속 실행시로드되어 프로그램을 빠르게 실행할 수 있습니다. 삭제하면 FFTW가 다시 생성하여 프로그램이 느리게 실행됩니다.

첫 번째 실행시 출력이 0 인 이유는이 계획 생성 단계 때문입니다. 그것은 FFTW FAQ에 설명되어 있습니다.

FFTW ++를 사용하지 않았기 때문에 FT를 수행하기 전에 계획이 생성되었는지 확인하기 위해 호출해야하는 정확한 방법을 모르겠습니다. 하지만 입출력 배열을 정의하고 데이터로 초기화하기 전에 즉시 호출해야합니다. 코드를 여러 번 실행해야하는 경우 지혜 파일을 보관해야합니다.

+0

감사합니다. 배열의 선언 바로 뒤에서 변환의 선언을 옮겼습니다. 이제 첫 번째 실행에서 올바르게 작동합니다. –

관련 문제