FFT执行错误（Nayuki vs Octave）

Question

这个问题可能是FFT知识和部分编程知识的一部分，但想到我会在这里发布它以查看您的想法。我试图在JavaScript中使用Project Nayuki's code实现一个斜坡过滤器，并且不能完全模仿我已经在C++（FFTW）和Octave/MATLAB中完成的工作。我将672到2048的初始数据数组填零，并在空间域中创建斜坡过滤器。下面是数据的图像之前和斜坡过滤器后，利用倍频的FFT：

而这里的八度代码：

% This script checks my FBP reconstruction code 
clear 

BaseFolder = "/home/steven/C++/TestSolutio/"; 
a = textread([BaseFolder "proj.txt"],"%f"); 
b = textread([BaseFolder "norm.txt"],"%f"); 
p = zeros(size(a)); 
for n = 0:499 
    p((672*n+1):(672*n+672)) = -log(a((672*n+1):(672*n+672)) ./ b); 
end 

dfan = (2.0*0.0625)/(80.0); 
FilterSize = (2*(672-1)) + 1; 
Np = 2048; 
FilterPadding = (Np-FilterSize-1)/2; 
FilterOriginal = zeros(FilterSize, 1); 
for f = 1:FilterSize 
    nf = (-672+1) + f - 1; 
    if(nf == 0) 
    FilterOriginal(f) = 1.0/(8.0 * dfan^2); 
    else 
    if(mod(nf,2) == 0) FilterOriginal(f) = 0; 
    else FilterOriginal(f) = -0.5/(pi*sin(nf*dfan))^2; 
    endif 
    endif 
end 
RampFilter = zeros(Np, 1); 
for f = 1:Np 
    if(f <= FilterPadding || f > (FilterSize+FilterPadding)) RampFilter(f) = 0; 
    else RampFilter(f) = FilterOriginal(f-FilterPadding); 
    endif 
end 
Filter = abs(fft(RampFilter)); 

proj_id = 0; 
ProjBuffer = zeros(Np,1); 
ProjPadding = (Np-672)/2; 
for f = 1:Np 
    if(f <= ProjPadding || f > (672+ProjPadding)) ProjBuffer(f) = 0; 
    else ProjBuffer(f) = p(672*proj_id+f-ProjPadding); 
    endif 
end 

ProjFilter = fft(ProjBuffer); 
ProjFilter = ProjFilter .* Filter; 
Proj = ifft(ProjFilter); 
ProjFinal = Proj((ProjPadding+1):(ProjPadding+672)); 

plot(1:672, p((672*proj_id+1):(672*proj_id+672))) 
axis([1 672 -5 10]) 

figure 
plot(1:Np, Filter) 

figure 
plot(1:672, ProjFinal) 

当我尝试这样做使用JavaScript，它看起来好像有一半的信号被翻转并添加到另一半，但我不知道真的发生了什么：

这里的JS功能：

function filterProj(proj){ 
    // Initialization variables 
    var padded_size = 2048; 
    var n_channels = 672; 
    var d_fan = (2.0*0.0625)/80.0; 
    // Create ramp filter 
    var filter_size = (2*(n_channels-1))+1; 
    var filter_padding = (padded_size - filter_size - 1)/2; 
    var ramp_filter = new Array(); 
    var nf; 
    for(f = 0; f < filter_size; f++){ 
    nf = (-n_channels+1) + f; 
    if(nf == 0) ramp_filter.push(1.0/(8.0*Math.pow(d_fan,2.0))); 
    else { 
     if(nf % 2 == 0) ramp_filter.push(0.0); 
     else ramp_filter.push(-0.5/Math.pow((Math.PI*Math.sin(nf*d_fan)),2.0)); 
    } 
    } 
    // Pad filter with zeros & transform 
    var filter_real = new Array(); 
    var filter_img = new Array(); 
    var filter = new Array(); 
    for(f = 0; f < padded_size; f++){ 
    if(f < filter_padding || f > (filter_size+filter_padding-1)){ 
     filter_real.push(0.0); 
    } 
    else { 
     filter_real.push(ramp_filter[(f-filter_padding)]); 
    } 
    filter_img.push(0.0); 
    } 
    transform(filter_real, filter_img); 
    for(f = 0; f < padded_size; f++){ 
    filter_real[f] = Math.abs(filter_real[f]); 
    } 
    // For each projection: 
    // Pad with zeros, take FFT, multiply by filter, take inverse FFT, and remove padding 
    var proj_padding = (padded_size - n_channels)/2; 
    for(n = 0; n < 500; n++){ 
    var proj_real = new Array(); 
    var proj_img = new Array(); 

    for(f = 0; f < padded_size; f++){ 
     if(f < proj_padding || f >= (n_channels+proj_padding)){ 
     proj_real.push(0.0); 
     } 
     else { 
     proj_real.push(proj[(n_channels*n + (f-proj_padding))]); 
     } 
     proj_img.push(0.0); 
    } 
    transform(proj_real, proj_img); 
    for(f = 0; f < padded_size; f++){ 
     proj_real[f] *= filter_real[f]; 
    } 
    inverseTransform(proj_real, proj_img); 

    for(f = 0; f < n_channels; f++){ 
     proj[(n_channels*n+f)] = (d_fan*proj_real[(proj_padding+f)])/padded_size; 
    } 
    } 
} 

任何帮助/建议将不胜感激！

更新

作为一个例子，这里的FFT后的斜坡过滤器，使用相同的空间域斜坡过滤器输入：

Answer 1

一些时间和调查，我发现后，那么我的错误不是FFT相关的，而是一个简单的复杂数学错误。在C++/Matlab/Octave中，复杂数据类型会重载abs（）函数来计算复杂的幅度。但是，Nayuki代码将输入/输出数据分解为其真实和复杂的部分，因此需要手动计算复杂的数值。总之，将下面一行：

filter_real[f] = Math.abs(filter_real[f]); 

这一个：

filter_real[f] = Math.sqrt(Math.pow(filter_real[f],2) + Math.pow(filter_img[f],2)); 

解决我的所有问题。