stevend12
stevend12

Reputation: 75

FFT implementation error (Nayuki vs Octave)

This question is probably part FFT knowledge and part programming knowledge, but figured I'd post it here to see what you think. I'm trying to implement a ramp filter in JavaScript using Project Nayuki's code and can't quite imitate what I've already done in C++ (FFTW) and Octave/MATLAB. I am zero-padding the initial data array of 672 to 2048 and creating the ramp filter in the spatial domain. Here's images of the data before and after the ramp filter, using Octave's FFT:

Before FFT (Octave) After FFT (Octave)

And here's the Octave code:

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

When I try to do this using JavaScript, it looks as if half the signal is flipped and added to the other half, but I don't know what's really happening:

Before FFT (JS) After FFT (JS)

Here's the JS function:

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;
    }
  }
}

Any help/advice would be greatly appreciated!

Update

As an example, here's the ramp filters after FFT, using the same spatial-domain ramp filter as input:

Filters

Upvotes: 3

Views: 346

Answers (1)

stevend12
stevend12

Reputation: 75

Well after some time and investigation I found that my error was not FFT related but a simple complex math error. In C++/Matlab/Octave, the complex data types overload the abs() function to calculate the complex magnitude. However, the Nayuki code splits the input/output data into its real and complex parts, so the complex magnitude needs to be calculated manually. In summary, replacing this line:

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

with this one:

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

solved all my problems.

Upvotes: 1

Related Questions