3d c2c fft с библиотекой fftw

Я пытаюсь выполнить 3D FFT с библиотекой FFTW, но у меня есть некоторые трудности с обратным преобразованием.

Сначала я делаю преобразование предисловия через:

fftwf_plan_dft_3d(_dimensions[0], _dimensions[1], _dimensions[2], (fftwf_complex*)_inputBuffer, (fftwf_complex*)_outputBuffer, FFTW_FORWARD, FFTW_ESTIMATE);

Хотя мои данные являются реальными данными, я использую преобразование от сложного к сложному, так как хочу заменить его позже opencl fft, который поддерживает только преобразования от сложного к сложному.

В трехмерном пространстве Фурье я использую очень простой фильтр нижних частот:

for all x, y, z:

// global position of the current bin
int gid = (y * w + x) + (z * w * h);

// position of the symmetric bin
vec3 conPos(M - x - 1, N - y - 1, L - z - 1);

// global position of the symmetric element
int conGid = (conPos.y * w + conPos.x) + (conPos.z * w * h);

if (sqrt(x * x + y * y + z * z) > 500)
{
    complex[gid].real = 0.0f;
    complex[gid].imag = 0.0f;
    complex[conGid].real = 0.0f;
    complex[conGid].imag = 0.0f;
}

Наконец, обратное преобразование:

fftwf_plan_dft_3d(_dimensions[0], _dimensions[1], _dimensions[2], (fftwf_complex*)_inputBuffer, (fftwf_complex*)_outputBuffer, FFTW_BACKWARD, FFTW_ESTIMATE);
// normalization ...

Результат не такой, как я ожидал. После обратного преобразования не все мнимые части равны нулю, как предполагалось.

Насколько я понимаю, после прямого преобразования реальных данных используется только половина общего размера буфера, а в другой половине нет сопряженных комплексных значений. (см.: c2c с реальными данными). мой собственный до обратного преобразования, но я не мог найти подсказку в документах fftw, какая половина вычисляется, а какая нет.

Я написал очень простой 2D-тест для просмотра этой симметрии в пространстве Фурье:

int w = 4;
int h = 4;
int size  = w * h;

cl_float rawImage[16] = ...; // loading image

fftwf_complex *complexImage = (fftwf_complex*) fftwf_malloc(sizeof(fftwf_complex) * size);
fftwf_complex *freqBuffer = (fftwf_complex*) fftwf_malloc(sizeof(fftwf_complex) * size);

for (int i = 0; i < size; i++)
{
    complexImage[i][0] = rawImage[i]; complexImage[i][1] = 0.0f;
}

fftwf_plan forward = fftwf_plan_dft_2d(w, h, complexImage, freqBuffer, FFTW_FORWARD, FFTW_ESTIMATE);

fftwf_execute(forward);

for (int y = 0; y < h; y++)
{
    for (int x = 0; x < w; x++)
    {
        int gid = y * w + x;
        qDebug() << gid << "real:" << freqBuffer[gid][0] << "imag:" << freqBuffer[gid][1];
    }
}

Это дает мне следующий результат:

gid
0    real 3060 imag 0 
1    real 510 imag 510 
2    real 0 imag 0 
3    real 510 imag -510 
4    real 510 imag 510 
5    real 0 imag -510 
6    real 0 imag 0 
7    real -510 imag 0 
8    real 0 imag 0 
9    real 0 imag 0 
10   real 0 imag 0 
11   real 0 imag 0 
12   real 510 imag -510 
13   real -510 imag 0 
14   real 0 imag 0 
15   real 0 imag 510

Насколько я понимаю, симметричных значений нет. Почему?

Было бы неплохо, если бы кто-нибудь подсказал мне.

Привет

Волк


3d fft fftw
person DerHandwerk    schedule 29.04.2012    source источник

Ответы (2)


arrow_upward
1
arrow_downward

Если вам нужен строго реальный результат после обратного БПФ (за вычетом обычного числового шума от использования арифметики конечного размера), вы должны убедиться, что входные данные, которые вы вводите в полный IFFT, полностью сопряженно-симметричны (последняя половина вектора - это зеркальное комплексное сопряжение первой половины). Не похоже, чтобы вы заставили свои данные быть такими.

person hotpaw2    schedule 29.04.2012
comment
Я обновил свой код (см. выше), чтобы сохранить симметрию, но он тоже не работает должным образом. - person DerHandwerk; 30.04.2012

arrow_upward
1
arrow_downward

Эта ссылка вводит в заблуждение. ДПФ реального сигнала (в общем случае) не приводит к тому, что половина выходных отсчетов равна нулю. Он просто накладывает на них (сопряженную) симметрию.

Таким образом, в вашем коде фильтра вы манипулируете только половиной выходных значений, которые должны быть. Каждый раз, когда вы манипулируете выходным бином n, вам также необходимо манипулировать бином Nn (где N — длина ДПФ), чтобы поддерживайте симметрию, которая даст вам реальный результат, когда вы примените обратное ДПФ.

Я бы посоветовал сначала решить гораздо более простую задачу — 1D-фильтр. Как только вы это сделаете, будет легко масштабироваться до 3D.

person Oliver Charlesworth    schedule 29.04.2012
comment
Это звучит разумно. Означает ли это, что в случае 1D FFT мне нужно получить только от n = 0 до N/2, поскольку на каждом этапе я фильтрую комплекс [n] и комплекс [N-n]? - person DerHandwerk; 29.04.2012
comment
@DerHandwerk: Не уверен, что понимаю, о чем вы спрашиваете. Одномерное ДПФ длины N будет иметь только N/2+1 независимых выходных отсчетов; соотношение y[n] = y[Nn]* (где * обозначает комплексно-сопряженное соединение) . - person Oliver Charlesworth; 29.04.2012
comment
Хорошо, насколько я понимаю, для 1D FFT мой нижний проход должен делать следующее: if amplitude > value then y[n] = 0; y[N - n - 1] = 0; end для n = 0 до N - 1 - person DerHandwerk; 30.04.2012