Computing group delay directly

In an attempt to better understand group delay, I tried writing simple code to compute it based directly on it's definition, rather than using the existing grpdelay function (with it's fancy efficient DFT based algorithm). My results weren't even close to those found with grpdelay, and I'm trying to figure where I went wrong.

Group delay is defined as the (negative) derivative of phase response, so given the filter coefficients a & b, and sampling rate fs I tried approximating the derivative like so:

[h f] = freqz(b, a, N, fs);
g1 = -diff(unwrap(arg(h)))./diff(f)/(2*pi);

Since arg(h) has units of radians and f has units of Hz, I would expect this to give the group delay in seconds, at each of the frequencies in f. To get the same using built-in matlab commands, I computed

g2 = grpdelay(b, a, f, fs)/fs;

The results of these two are not at all similar. What am I doing wrong?

[SAM Arani]

I happen to have same problem, Have you find the problem yet?