I've done the same code on both MATLAB and Python, but ifft2 returns different values?

Question

I've been trying to implement an homomorphic filter in frequency domain on both MATLAB and Python using OpenCV2 and NumPy, the MATLAB code gives the expected answer but the Python does not, the resulting image is very weird. I've tested all variables and came to the conclusion the only point there is a difference is the IFFT. On MATLAB, the results can be applied normally to the exp function and return the filtered original image expected, but the values of Python ifft are very different. I happened to see other posts with similar problems, but no satisfactory answer (perhaps i'm just bad at searching too...). The MATLAB code

function [ img_r ] = homomorphic( img, D0, n )

[N, M] = size(img);
img_bk = double(img);
img_bk = log(img_bk+1);

img_freq = fftshift(fft2(img_bk));
magA = uint8(10*log(1+abs(img_freq)));

cu = M/2;
cv = N/2;
Hf = zeros(N,M);
for v = 1:N
    dv = v - cv;
    for u = 1:M
        du = u - cu;
        D = sqrt(du*du + dv*dv);

        num = 1;
        if D > 0
            den = 1+((D0/D)^(2*n));
        else
            den = 0; %to replace +inf
        end

        if den ~= 0
            H = num/den;            
        else
            H = 0;
        end
        img_freq(v,u) = H*img_freq(v,u);

    end
end

magB = uint8(10*log(1+abs(img_freq)));
img_r = (ifft2(ifftshift(img_freq)));
img_r = exp(img_r);
img_r = uint8(img_r);

and the Python code (might have some bugs but overall works)

import numpy as np
import cv2

def homomorphic(img, D0, n=2):
    [N,M] = img.shape
    img_bk = np.log(1 + np.float64(img))
    img_freq = np.fft.fftshift(np.fft.fft2(img_bk))
    cu = M/2.0
    cv = N/2.0    

    for v in range(N):
        dv = v - cv
        for u in range(M):
            du = u - cu
            D = np.sqrt(du*du + dv*dv)
            if D != 0:
                a = 1.0 + (D0/D)**(2*n)
                H = 1/a
            else:
                print D
                H = 0
            img_freq[v][u] = H*img_freq[v][u]

    img_r = np.abs(np.fft.ifft2(np.fft.ifftshift(img_freq)))
    eimg = np.exp(img_r)
    eimg = np.uint8(eimg)
    return eimg

I really don't get it, why the results are so different? Does anyone have any idea?