function [f, df] = discreteGaussian(sig)
    x = floor(-3*sig):ceil(3*sig);
    f = exp(-0.5*x.^2/sig^2);
    f = f / sum(f);
    df = -x.*f/sig^2;