IQ signal FM descrimination using phase derivative I'm trying to demodulate an FM signal recorded as IQ samples from RTL-SDR at 2.4Msps with RTL-SDR centered on an FM station's frequency.

So the input signal is 8 bit IQ baseband at 2.4MSps.

The data is available here.

The FFT of the input signal looks fine and matches what I see in SDR-Sharp so I think the data is good.

My understanding is that I should be able to discriminate the FM in Octave by:


Removing the DC offset 
Low pass filtering it to 240KHz and decimating
Computing discrete derivative of the phase component


The FFT looks good after the low pass filtering and decimation.
My assumption is that applying a FIR low pass filter shouldn't distort phase (not sure if I'm correct)

I'm then computing discrete derivative of the phase of the signal by convoluting the phase signal with the discrete derivative's impulse response of [1,-1].

I'm expecting that the spectrum of the derivative signal would contain the composite FM signal with 0-15KHz occupied by the mono L+R component the 19KHz pilot tone and so on.
But I'm getting white noise in the FFT of the phase derivative (see below).

My Octave code is as follows

% 2.4MSps IQ baseband signal of an FM station from RTL-SDR recorded in SDR#
  samples = audioread("SDRSharp_20191229_133205Z_102102000Hz_IQ.wav", "native");
  csamples = arrayfun(@complex,samples(:,1),samples(:,2));

% 1) - remove DC offset
  csamples_ac = csamples - mean(csamples);   

% 2) - low pass filter and decimate to get complex 240KSps (240KHz Widband FM)
  csamples_filtered = decimate(csamples_ac, 10, "fir");

% 3) - descrete derivative of phase
  filtered_derivative = conv(angle(csamples_filtered), [1,-1]);

  plot(abs(fftshift(fft(filtered_derivative, 32768))));
  axis("tight");




Edit: Thanks to hotpaw2 and Dan Boschen I now know what my mistakes were and have a working Octave script. The main mistake was not unwrapping the phase, also when discriminating by convolution with [1,-1] one needs to discard the first and last samples from the result because of possible artifacts of the convolution which give huge values at the ends of the sequence, the third mistake was not converting to int16 and not getting rid of DC in the sound samples before using Octave's audioplayer function.

In case anyone ever needs such a script a working version is below.
(Note that the phase-based discrimination suffers from a low-pass filtering 1/f rolloff effect, for more details see this answer)

pkg load signal;

#reading SDR Sharp recording (2.4Msps,RAW,8bit,baseband)
samples = audioread("SDRSharp_20191229_183956Z_99500000Hz_IQ.wav", "native");
#one channel real, other channel - imaginary components
unshifted_signal = double(samples(:,1)) .+ double(samples(:,2))*j;

#removing DC, this is important, 
unshifted_signal = unshifted_signal - mean(unshifted_signal);

#a 1550 Hz frequency shift to move the real carrier center frequency as close to 0 as possible (de-trending the phase)
#the frequency shift is not absolutely necessary, small diviations don't affect the demodulation
#this would be station specific, selected to minimized the slope on plot(unwrp_phase)
freq_shift = exp(-j*1550*2*pi*[1:length(unshifted_signal)]/2400000)';
#freq_shift = 1;
signal = unshifted_signal.*freq_shift;

#low-pass filter and decimate down to 240KHz
rcv_240 = decimate(signal, 10);
#calculating the phase, and unwrapping it (unwrapping is important)
unwrp_phase = unwrap(angle(rcv_240));

#doing the actual demodulation (discrimination)
#not including the ends because of possible convolution artifacts
phase_drv = conv(unwrp_phase, [1,-1])(2:end-1);

#extract L+R submodulated at 0-15KHz by low-pass filtering and decimation
monoLplusR_unnormalized = decimate(phase_drv, 8);

#normalize the range for 16bits for audioplayer, it is important to have no DC here (removed earlier)
monoLplusR = int16(monoLplusR_unnormalized/max(abs(monoLplusR_unnormalized))*32766);

#play at 30Ksps (0-15KHz)
player = audioplayer(monoLplusR, 30000, 16);
play(player);
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().