Erreur d'implémentation FFT (Nayuki vs Octave)

Question

Cette question fait probablement partie de la connaissance FFT et de la connaissance de la programmation partielle, mais j'ai pensé que je l'afficherais ici pour voir ce que vous en pensez. J'essaye d'implémenter un filtre de rampe en JavaScript en utilisant Project Nayuki's code et je n'arrive pas à imiter ce que j'ai déjà fait en C++ (FFTW) et Octave/MATLAB. Je suis le remplissage zéro du tableau de données initial de 672 à 2048 et la création du filtre de rampe dans le domaine spatial. Voici les images des données avant et après le filtre à rampe, en utilisant la FFT d'Octave:

Et voici le code 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) 

Lorsque je tente de le faire en utilisant JavaScript, il On dirait que la moitié du signal est inversée et ajoutée à l'autre moitié, mais je ne sais pas ce qui se passe réellement:

est ici la fonction 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; 
    } 
    } 
} 

Toute aide/conseils serait grandement apprécié!

Mise à jour

À titre d'exemple, voici les filtres de rampe après FFT, en utilisant le même filtre de rampe spatiale domaine en entrée:

Answer 1

bien après un certain temps et d'enquête, je trouve que mon l'erreur n'était pas liée à la FFT mais à une simple erreur mathématique complexe. En C++/Matlab/Octave, les types de données complexes surchargent la fonction abs() pour calculer l'amplitude complexe. Cependant, le code de Nayuki divise les données d'entrée/sortie en ses parties réelles et complexes, ainsi l'ampleur complexe doit être calculée manuellement. En résumé, en remplaçant cette ligne:

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

avec celui-ci:

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

résolu tous mes problèmes.