# 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:

1. Removing the DC offset

2. Low pass filtering it to 240KHz and decimating

3. 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#
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;

#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);

• Someone just liked this post (perhaps you) and I realized in seeing it that is is very relevant to what you are now doing so wanted to make sure you saw it: dsp.stackexchange.com/questions/38601/… Dec 29 '19 at 20:42
• @DanBoschen, I was trying to apply small frequency shifts as you suggested for fine-tuning (by multiplying the complex filtered 240KHz signal by exp(-j*2*pinfshift/fsampling)) and I observed fading amplitude when deviating further from fc. But more interestingly there was a "dead spot" near (1800Hz off) my original center frequency. The amplitude of the unwrapped phase was much ~100 times lower at that point. I wonder why this happened, just something wrong with my octave code or something well-known.
– axk
Dec 29 '19 at 22:06
• Are you only looking at the real (or imaginary) output? That is the only way you could see amplitude fading by simply multiplying by the exponential frequency -or do you mean you see amplitude fading in your demodulated signal? Dec 29 '19 at 22:09
• Derive and plot how your derivative response would actually be for your implementation mapping instantaneous frequency versus time to a magnitude versus time quantity. For example look at the plot I gave below where the discrimator response is linear and useful at the origin, but as you shift away from center with a frequency offset the discriminator gets non linear (when it is flat at the top of that sinusoidal shape there would be no output for a smaller FM deviation signal). Map what you are doing specifically to see if you have such a non-linearity. Dec 29 '19 at 23:05
• Also remember to hard limit your waveform prior to demod -- any AM will translated to AM of the demodulated signal but given this is FM that would be an undesired interference. Dec 29 '19 at 23:36

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().

• @axk This may be the better (and if true then the one to mark as correct) answer. He makes an excellent point that I didn't consider is your phase isn't unwrapped. If you are processing in Octave you can use the unwrap command to test this unwrap(angle(x)). Otherwise the difference is producing a large (noise) impulse everytime you cross the $2\pi$ boundary. The points I gave are still valid considerations but I think hotpaw2 may have gotten to the bottom of the noise issue specifically. Dec 29 '19 at 17:51
• That indicates you have a frequency offset error (your estimate of what the actual carrier frequency or center frequency of the FM modulation is off by that slope. Just detrend that and amplify the result and you should be there. Dec 31 '19 at 22:17
• (I describe that particular aspect in more detail in the link I pointed to you in the comment under your question) Dec 31 '19 at 22:49
• The SNR you would get is effected by the magnitude to frequency slope of your discriminator - you want to have the highest slope within its usable range —- have you looked at that yet to see that you maximized it? Just use some test sinusoidal FM tones if you are not confident with the math involved or to confirm your math Jan 1 '20 at 1:06
• You might want to try de-trending in the time domain: estimate the slope of the phase and subtract it out (subtract a line with that same slope). Jan 1 '20 at 1:52

At 240 KHz the phase increment per unit sample is very small so your derivative is a high pass filter with a cut-off well above 15KHz. (Review freqz([1 -1]) to see what I mean.) You have a frequency discriminator but its phase slope is very small in converting frequency variation to amplitude variation for your signals of interest.

You can try increasing the number of samples between differences [1 0 0 0 0 -1] for examples or low pass filter and then decimate your phase domain waveform to continue with your approach if you really only want what is in the 0 to 15KHz band.

Also consider the implementation using the imaginary result of the complex conjugate multiplication of the waveform with a delayed version of itself directly without having to compute the angle (as another option).

This would be $$I[n]Q[n-m]-I[n-m]Q[n]$$

Where m is increased to increase the slope (gain) of the discriminator. Note how with m=3 below the slope of output magnitude versus frequency variation on the input is increased. So I am suspecting your slope is way too low.

• Yes see my comment under his. It isn't a contradition, with the shallow slope you will have less output SNR. Ultimately you want to match your slope to the modulation index which is driven by the frequency deviation in the modulated signal. So review what the frequency deviation is for broadcast FM and then derive your discriminator as a linear slope of Frequency In to Magnitude Out and see if your slope is optimally set for that deviation. Dec 29 '19 at 17:56
• Regarding strongly attenuated near DC--- consider what a "DC" phase is; 0 frequency meaning unmodulated carrier. You will always have a high pass response on the phase given frequency is the derivative of phase. This is why pre-emphasis is employed in FM broadcast to equalize for this effect, reducing hissing at high frequencies (pre-emphasize the higher frequencies at transmit and then de-emphasize them in the receiver). FM is direct modulation based on frequency but SNR is directly proportional to the phase. (so FM has a natural 1/f roll-off on the modulated signal SNR) Dec 29 '19 at 18:03
• (So if you do a delay and complex conjugate multiply frequency discriminator you won't have the phase wrapping issue.....but I think @hotpaw2 's response was insightful so if that does turn out to be the case, his would be the "correct" answer). Dec 29 '19 at 18:06
• Trying to make sense of the FFT I'm getting for the derivative with the impulse response of [1, 0, -1] vs just [1,-1]. For [-1,0,-1] I'm getting a curve similar to the frequency response of [1, 0, -1] with the signal superimposed on top of it i.imgur.com/4cvdoL2.png . As if the time domain was a delta function, but the impulse response starts in the middle of the frequency range.
– axk
Jan 4 '20 at 20:14
• Have you tried a test signal of a simple FM tone sampled at your same rate and then passed through your discriminator? You can do that to easily derive your discriminator function of output magnitude vs frequency in --- that should give you a lot of insight into both the maximum deviation your discrimnator will support and that the math is making sense (could have scaling issues etc). Jan 4 '20 at 20:29