Filtered backprojection in MATLAB and designing filter

I'm trying to write my own MATLAB code to compute the inverse radon transform (iradon) and thus far I have managed to successfully reconstruct an image using a ramp filter, a hamming window, and also using convolution of the 1D projections in the spatial domain with a window h in my code based on the textbook by Kak and Shakey. However, I think that if I take the FFT of the window h and multiply that by the FFT of the 1D projections, I should be able to obtain the same reconstruction. Unfortunately, what I do get is a mess.

function [out] = myfbp4(arg2);


ph = phantom();
sino = radon(phantom,0:0.1:179);

rho2 = 183; % max rho
rho1 = -183; % min rho;
step = 1;
npos = 367;
dtheta = 0.1;
angles = deg2rad(0:dtheta:179);
nproj = length(angles);



n1 = ceil(log(npos)/log(2)); 
N = 2^n1;                  % smallest power of 2 greater than npos (for fft efficiency)
N = 1024;   % for zero padding
fs = 1; % sampling frequency
T = 1;  % sample spacing

% grid for reconstructed image
x1 = rho1:step:rho2;
y1 = rho1:step:rho2;

[yyy, xxx] = ndgrid(-y1, x1);


% build filter h according to Kak and Shakey 
h = -floor(npos/2):T:floor(npos/2);

for i = 1:length(h)
    if (mod(h(i),2) == 1) 
        h(i) = -1./(h(i)^2 * pi^2 * T^2);
    else
        h(i) = 0;
    end    
    h(ceil(npos/2)) = 1/(4*T^2);
end

%%%%%%%%%%%   RAMP FILTER   %%%%%%%%%%%%%%%%
%%%%%%%%%% this is un needed when using h %%
%%%%%%%%%% see below... %%%%%%%%%%%%%%%%%%%%
rampFilterNum = [0:1:N/2-1 -N/2:1:-1];
rampFilterAbs = abs(rampFilterNum);
rampFilter = rampFilterAbs *fs/N; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% positions made to correspond to N point fft
drho = (rho2 - rho1) / (npos-1);
rho  = rho1 + (0:(npos-1)) * drho;
positions = zeros(nproj, length(rho));
for i = 1:nproj,
  positions(i, :) = rho;
end

if (strcmp(arg2,'filter')) 
    % compute FT of h and multiply by fft projections
    fftProj = fft(sino, N);
    hfft = fft(h,N);
    fftProjFiltered = bsxfun(@times, hfft', fftProj);
    ifftProj = real(ifft(fftProjFiltered));
    filteredProjections = ifftProj;
end

if (strcmp(arg2,'conv')) 
    % make image my convolution of projections with h
    for iproj = 1:nproj
        sino(:, iproj) = conv(sino(:,iproj), h, 'same');
    end
    filteredProjections = sino;
end


% display the image through backprojection
fdata = zeros(size(xxx));
for iproj = 1:nproj
    theta = angles(iproj);
    rho1 = xxx*cos(theta) + yyy*sin(theta); % rotate coordinate system by theta
    %r = x1;
    r = positions(iproj,:);

    fdata1 = filteredProjections(1:npos,iproj); % filtered projections
    %fdata1  interp1(
    fdata2 = interp1(r, fdata1, rho1, 'linear', 0); 
    fdata = fdata + deg2rad(dtheta) * fdata2; %theta*fdata2;
end

out = fdata;
end

Just running out = myfbp4('conv') or myfbp4('filter') will show you the different results. It seems that the convolution works fine, but the filtering approach isn't working as I'd hoped.

Can anyone see the problem? (apologies if there is any redundant code, I tried to cut out most of it... I should also mentioned that this code was borrowed from somewhere and modified a bit, but I don't remember where I found it).

Thanks in advance

EDIT: Problem solved. The issue was that I was not taking the absolute value of the fourier transform of the window h to obtain the frequency window. For those who find this, hfft = abs(fft(h,N)) should replace hfft = fft(h, N).

Problem solved. The issue was that I was not taking the absolute value of the fourier transform of the window h to obtain the frequency window. For those who find this, hfft = abs(fft(h,N)) should replace hfft = fft(h, N).