javascript - FFT 实现错误(Nayuki vs Octave)

Question

这个问题可能部分是 FFT 知识和部分编程知识，但我想我会把它贴在这里看看你的想法。我正在尝试使用 Project Nayuki's code 在 JavaScript 中实现渐变过滤器并且不能完全模仿我在 C++ (FFTW) 和 Octave/MATLAB 中已经完成的工作。我将 672 的初始数据数组零填充到 2048，并在空间域中创建斜坡滤波器。下面是斜坡滤波器前后的数据图像，使用 Octave 的 FFT:

这是 Octave 代码:

% 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));

解决了我所有的问题。