function h = compute_entropy(f, lim)
% computes the entropy for 2D probability density function f
% lim specifies the bounds of the square integration domain
g = @(x, y) -f(x,y).* log2(f(x,y));
h = integral2(g,-lim ,lim ,-lim ,lim, 'AbsTol', 1e-3);