Computing blocking on digital images/ videos

I am trying to implement a visual quality assessment metric which calculates blocking artifact for digital images/ videos. The metric is proposed in "A Blind Reference-Free Blockiness Measure" for Chunhua Chen and Jeffrey A. Bloom (plz find the paper attached).

Denote the given pixel array (an image or a video frame) as I(x,y), where 0 <= x <= Sx-1, 0 <= y <= Sy-1, Sx is the height and Sy is the width of the given image pixel array.

The proposed metric computes the vertical blockiness BMv as follow: Step 1: The absolute difference between horizontally adjacent pixels is firstly calculated. This operation results in a gradient image Dc(x,y), which is given by:

Dc(x,y) = I(x,y) - I(x+1m y), where 0 <= x <= Sx-2, 0 <= y <= Sy-1.

function [horFourier, verFourier] = DFT(hProfile, vProfile)

horFourier = abs(fft(hProfile));
verFourier = abs(fft(vProfile));

end

Step 2: Each element in the gradient image is normalized.

My Matlab Code:

function [horDiffNorm, verDiffNorm] = normalization(H, V)

      H = double(H);
      V = double(V);

      [rowsH, colsH, depthH] = size(H);
      [rowsV, colsV, depthV] = size(V);

      for jH = 1:1:rowsH
          for iH = 1:1:colsH
              sumH = 0;
              for zH = 1:1:colsH
                  if (zH ~= iH)
                      sumH = sumH + power(H(jH, zH), 2);
                  end
              end
              horDiffNorm(jH, iH) = H(jH, iH)/ sqrt(double(1/(2*(colsH-1))*sumH));
           end
      end

Steps 3: An average operation along each column is then applied, resulting in a 1-D array pf length Sx-1 called horizontal profile.

function [horProfile, verProfile] = calculateProfile(H, V)

      [rowsH, colsH] = size(H);
      [rowsV, colsV] = size(V);

      for iH=1:1:colsH
          sumH = 0;
          for jH=1:1:rowsH
              sumH = sumH + H(jH, iH);
          end
          horProfile(iH) = sumH/ rowsH; 
      end

The result: </matlabcentral/answers/uploaded_files/138027/signal1.png>

Steps 4: Take 1-D DFT to PH and just consider the magnitude FPH of the DFT coefficients.

function [horFourier, verFourier] = DFT(hProfile, vProfile)

horFourier = abs(fft(hProfile));
verFourier = abs(fft(vProfile));

end

Step 5: Computing vertical blockiness measure

    function [verticalBlockiness] = verticalBlockinessMeasure(hFourier)
        global maxHorBS
        for k = 2:1:maxHorBS
            sum = 0;
            for i = 1:1:k-1
                sum = sum + power(hFourier(i*(length(hFourier)-1)/k-1), 2);
            end  
            BMv(k-1) = 1/hFourier(0)*(sqrt(1/(k-1)*sum));
        end
        verticalBlockiness = max(BMv);
    end

When I run the code, I get below error: