逆FFT实现问题：JS生成的水平带倾斜，与Mathematica结果不符

我基于Radix-2实现了逆快速傅里叶变换（IFFT）用于图像处理：将尺寸为N×N（N为平方数）的图像转换为长度为N×N的向量。当输入矩阵仅在中心上方和对称下方各有一个1时，逆变换结果为垂直带，频率与1的位置匹配，这部分是正确的。但当1位于中心左侧和右侧时，得到的水平带存在微小倾斜；而使用Mathematica内置的InverseFourier函数执行相同操作，得到的水平带完全水平。请问以下JavaScript实现是否存在问题？

实现代码 #

IFFT核心函数 #

function invFFT_radix2(out, start, input, N, offset, s ) {
    if (N === 1) {
        out[start] = input[offset];
    } else {
        invFFT_radix2(out, start,     input, N/2, offset,   2*s);
        invFFT_radix2(out, start+N/2, input, N/2, offset+s, 2*s);
        for (var k = 0; k < N/2; k++) {
            var twiddle = cisExp(2*Math.PI*k/N);
            var t = out[start+k];
            out[start+k]     = t.plus ( twiddle.times(out[start+k+N/2]) );
            out[start+k+N/2] = t.minus( twiddle.times(out[start+k+N/2]) );
        }
    }
}

调用方式 #

invFFT_radix2(out, 0, input, input.length, 0, 1 )

辅助函数与说明 #

function cisExp(x) { // e^ix = cos x + i*sin x
    return new Complex(Math.cos(x), Math.sin(x));
}

• input为输入向量，out用于存储变换结果，所有元素均为Complex类型（包含plus、minus、times等基础运算方法）。

结果对比 #

• JS实现生成的图像顶部放大图：
• Mathematica生成的图像顶部放大图：

JS生成的图像左侧为白线，右侧变为黑线。

你的Radix-2 IFFT实现存在两个关键问题，导致水平带出现倾斜：

1. 旋转因子符号错误
   逆FFT的旋转因子应为e^(-2πik/N)，但你当前传入cisExp的参数是2*Math.PI*k/N，符号完全相反。这个错误会导致频率分量的相位偏移，直接引发水平带的倾斜。修正后的旋转因子计算应为：

   var twiddle = cisExp(-2*Math.PI*k/N);
   

2. 缺少归一化步骤
   Mathematica的InverseFourier默认会对结果做归一化（除以总长度N），你的实现没有这一步。虽然这不会直接导致倾斜，但会影响输出幅值，建议在IFFT完成后对out中的每个元素除以input.length。

另外，需注意图像的行列遍历顺序：处理2D图像FFT时，通常需先对每一行做1D FFT，再对每一列做1D FFT（或反之）。如果输入向量是按行优先展开的，要确保IFFT的维度顺序与Mathematica一致，否则也可能出现方向或相位偏差。

修正后的核心函数示例：

function invFFT_radix2(out, start, input, N, offset, s ) {
    if (N === 1) {
        out[start] = input[offset];
    } else {
        invFFT_radix2(out, start,     input, N/2, offset,   2*s);
        invFFT_radix2(out, start+N/2, input, N/2, offset+s, 2*s);
        for (var k = 0; k < N/2; k++) {
            // 修正旋转因子符号
            var twiddle = cisExp(-2*Math.PI*k/N);
            var t = out[start+k];
            var temp = twiddle.times(out[start+k+N/2]);
            out[start+k]     = t.plus(temp);
            out[start+k+N/2] = t.minus(temp);
        }
    }
}

// 调用后添加归一化步骤
for (let i = 0; i < out.length; i++) {
    out[i] = out[i].divide(out.length);
}

内容的提问来源于stack exchange，提问作者Denis Cousineau