Странный, но близкий FFT и ifft изображения в переполнении стека

Я написал программу, которая загружает, сохраняет и выполняет FFT и IFFT для черно-белых изображений PNG. После большой отладочной головной боли я наконец-то получил внятный вывод и обнаружил, что он искажает исходное изображение.

вход:

FFT:

IFFT:

Насколько я тестировал, данные пикселей в каждом массиве сохраняются и преобразуются правильно. Пиксели хранятся в двух массивах: «данные», которые содержат ч / б значение каждого пикселя, и «комплексные данные», которые в два раза длиннее «данных» и хранят реальное ч / б значение и мнимые части каждого пикселя в чередующихся индексах. Мой алгоритм fft работает с массивом, структурированным как ‘complex_data’. После кода для чтения команд от пользователя вот такой код:

if (cmd == "fft")
{
if (height > width) size = height;
else size = width;

N = (int)pow(2.0, ceil(log((double)size)/log(2.0)));

temp_data = (double*) malloc(sizeof(double) * width * 2); //array to hold each row of the image for processing in FFT()

for (i = 0; i < (int) height; i++)
{
for (j = 0; j < (int) width; j++)
{
temp_data[j*2] = complex_data[(i*width*2)+(j*2)];
temp_data[j*2+1] = complex_data[(i*width*2)+(j*2)+1];
}
FFT(temp_data, N, 1);
for (j = 0; j < (int) width; j++)
{
complex_data[(i*width*2)+(j*2)] = temp_data[j*2];
complex_data[(i*width*2)+(j*2)+1] = temp_data[j*2+1];
}
}
transpose(complex_data, width, height); //tested
free(temp_data);
temp_data = (double*) malloc(sizeof(double) * height * 2);
for (i = 0; i < (int) width; i++)
{
for (j = 0; j < (int) height; j++)
{
temp_data[j*2] = complex_data[(i*height*2)+(j*2)];
temp_data[j*2+1] = complex_data[(i*height*2)+(j*2)+1];
}
FFT(temp_data, N, 1);
for (j = 0; j < (int) height; j++)
{
complex_data[(i*height*2)+(j*2)] = temp_data[j*2];
complex_data[(i*height*2)+(j*2)+1] = temp_data[j*2+1];
}
}
transpose(complex_data, height, width);

free(temp_data);
free(data);

data = complex_to_real(complex_data, image.size()/4); //tested
image = bw_data_to_vector(data, image.size()/4); //tested
cout << "*** fft success ***" << endl << endl;void FFT(double* data, unsigned long nn, int f_or_b){ // f_or_b is 1 for fft, -1 for ifft

unsigned long n, mmax, m, j, istep, i;
double wtemp, w_real, wp_real, wp_imaginary, w_imaginary, theta;
double temp_real, temp_imaginary;

// reverse-binary reindexing to separate even and odd indices
// and to allow us to compute the FFT in place

n = nn<<1;
j = 1;
for (i = 1; i < n; i += 2) {
if (j > i) {
swap(data[j-1], data[i-1]);
swap(data[j], data[i]);
}
m = nn;
while (m >= 2 && j > m) {
j -= m;
m >>= 1;
}
j += m;
};

// here begins the Danielson-Lanczos section

mmax = 2;
while (n > mmax) {
istep = mmax<<1;
theta = f_or_b * (2 * M_PI/mmax);
wtemp = sin(0.5 * theta);
wp_real = -2.0 * wtemp * wtemp;
wp_imaginary = sin(theta);
w_real = 1.0;
w_imaginary = 0.0;
for (m = 1; m < mmax; m += 2) {
for (i = m; i <= n; i += istep) {
j = i + mmax;
temp_real = w_real * data[j-1] - w_imaginary * data[j];
temp_imaginary = w_real * data[j] + w_imaginary * data[j-1];

data[j-1] = data[i-1] - temp_real;
data[j] = data[i] - temp_imaginary;
data[i-1] += temp_real;
data[i] += temp_imaginary;
}
wtemp = w_real;
w_real += w_real * wp_real - w_imaginary * wp_imaginary;
w_imaginary += w_imaginary * wp_real + wtemp * wp_imaginary;
}
mmax=istep;
}}

Мой ifft такой же, только с f_or_b, установленным в -1 вместо 1. Моя программа вызывает FFT () в каждой строке, транспонирует изображение, снова вызывает FFT () в каждой строке, а затем транспонирует обратно. Может быть, ошибка с моей индексацией?