IQ signal FM descrimination using phase derivative
I simply use unwrapped atan2(IQ(i)) - atan2(IQ(i-1)) to estimate a discrete derivative, then low pass filter to below 15 kHz.  Although with a shallow slope, the 1st order approximation to atan() given by Boschen will work just as well.  

Your noise might be due to not unwrapping the phase delta, or to not low pass filtering after doing the phase discrimination.  Not sure whether your angle() function is the same as atan2().