How do I quantize in JPEG compression?

I am studying the JPEG compression algorithm. I followed some simple versioning instruction to implement it in MatLab and I followed the quantization process.

So, in JPEG, I use 8x8 blocks as a unit to perform a direct transformation, and then quantize each block based on a quantization matrix (or divide by N for simplicity in the instructions).

I have implemented DCT myself and it works the same as inline dct2, so I think there is no problem in my DCT code.

function transformed_matrix = dctTransform(block, inverse)
m_size = size(block, 1);
A = zeros(m_size);

for i = 0:m_size - 1
    for j = 0:m_size - 1
        if i == 0
            a = sqrt(1 / m_size);
        else
            a = sqrt(2 / m_size);
        end            
        A(i + 1, j + 1) = a * cos(pi * (j + 0.5) * i / m_size);
    end
end

if inverse == true
    transformed_matrix = A' * block * A;
else
    transformed_matrix = A * block * A';
end

end

      

        
        
        
      

    

Then I start implementing the quantization, I made a simple version like below (Only for grayscale now ..):

function quantized_matrix = quantize(block, quality, inverse, mode)
m_size = size(block, 1);
N = 16;
DEFAULT_QUANTIZATION_MATRIX = ...
    [16 11 10 16 24 40 51 61
     12 12 14 19 26 58 60 55
     14 13 16 24 40 57 69 56
     14 17 22 29 51 87 80 62
     18 22 37 56 68 109 103 77
     24 35 55 64 81 104 113 92
     49 64 78 87 103 121 120 101
     72 92 95 98 112 100 103 99] * quality;

% check for input size and mode
if strcmp(mode, 'default') && m_size == 8
    if inverse == true
        quantized_matrix = block .* DEFAULT_QUANTIZATION_MATRIX;
    else
        quantized_matrix = round(block ./ DEFAULT_QUANTIZATION_MATRIX);
    end
else
    if inverse == true
        quantized_matrix = block * N;
    else
        quantized_matrix = round(block / N);
    end
end

end

      

        
        
        
      

    

My main program code

I = im2double(imread('../images/lena.bmp'));
block_size = 8;
fun = @(block_struct) quantize(dctTransform(block_struct.data, false), 1, false, 'defualt')
fun2 = @(block_struct) dctTransform(block_struct.data, false)
fun3 = @(block_struct) dct2(block_struct.data)
I2 = blockproc(I, [block_size block_size], fun2);
I3 = blockproc(I, [block_size block_size], fun3);
I4 = blockproc(I, [block_size block_size], fun);
subplot(2,2,1), imshow(I, []), title('The Original Image');
subplot(2,2,2), imshow(I2, []), title('The DCT Image');
subplot(2,2,3), imshow(I3, []), title('The builtin DCT Image');
subplot(2,2,4), imshow(I4, []), title('The Quantized Image');

      

        
        
        
      

    

There is no difference between my DCT and the built-in DCT implementation, so I think there must be something wrong with my quantization implementation. I checked the DCT calculation result, most of the numbers in the matrix are very small and so I finally have a black image (all rounded to 0). Is there a misunderstanding of JPEG compression from my implementation?

Any help is appreciated.