DCT using FFT results in complex result

I'm trying to implement a DCT10 accourding to this paper https://www.researchgate.net/publication/330405662_PittPack_Open-Source_FFT-Based_Poisson%27s_Equation_Solver_for_Computing_With_Accelerators (section "Neumann Boundary Condition"). However I have the problem that after performing the FFT and half-sample shifting, the result is not purely real (which i think it should be, right ?) Therefore when truncating the imaginary part, the mentioned reverse transform will not result in my original values.

Here is my Matlab code (DCT in first dimension):

function X_dct = dct_type2(x_sig)
    N = size(x_sig);

    % shuffle to prepare for FFT
    x_hat = zeros(N);
    for m = 1 : N(2)
        for n = 1 : (N(1) / 2)
            x_hat(n, m) = x_sig((2 * n) - 1, m);
            x_hat(N(1) - n + 1, m) = x_sig(2 * n, m);
        end
    end

    % perform FFT
    X_hat_dft = fft(x_hat, N(1), 1);

    % apply shifting by half-sample
    X_dct = zeros(N);
    for m = 1 : N(2)
        for k = 1 : N(1)
            X_dct(k, m) = 2 * exp(-1i * (pi * (k-1)) / (2 * N(1))) * X_hat_dft(k, m);
        end
    end
end

Can somebody explain what is the problem here ? Or is my assumption wrong that the result should be purely real ?

So it turns out that it is correct to drop the non-zero imaginary part using this technique, even though it intuitively appeared wrong to me. That the reverse transform did't recorver the original values was merely a scaling issue of the frequency components.