Understanding Integration Scale in Multi-Scale Harris Corner Detector
I think the first smoothing, by $\sigma_D$, is only done to get more stable derivatives whereas in the second step the convolution by a Gaussian with $\sigma_I$ is done to establish the 'scale-space' in which the operator is applied. 

Ignoring $\sigma_D$, this looks like this in Matlab:

dx = [-1 0 1; -1 0 1; -1 0 1];  % Simple mask for derivative 
Ix = conv2(im, dx, 'same');     % Convolve against image
Iy = conv2(im, dx', 'same');    % Again for y-direction with transposed mask

% Gaussian filter, as defined by Kovesi
% www.csse.uwa.edu.au/~pk/Research/MatlabFns/index.html
sigma = 1;
g = fspecial('gaussian', max(1, fix(6*sigma)), sigma);  

% Determine smoothed squared image derivatives via convolution
Ix2 = conv2(Ix.^2, g, 'same'); 
Iy2 = conv2(Iy.^2, g, 'same');
Ixy = conv2(Ix.*Iy, g, 'same');

% Metric to define the corner score
cim = (Ix2.*Iy2 - Ixy.^2) - k*(Ix2 + Iy2).^2;    % Original measure
cim = (Ix2.*Iy2 - Ixy.^2)./(Ix2 + Iy2 + eps);    % Kovesi measure