Real-time FM demod results in clicking and strangely-shaped output

Question

hardware engineer here dipping his toes into software DSP and running into some issues.

I'm attempting to demodulate in realtime a simple narrowband FM signal with samples from an RTL-SDR, and play the resulting audio through my speakers.

My Sampling Approach

The way my code is set up currently, I have a Python sounddevice running at an audio sample rate of 24khz which makes a call to my getSamples() function every time it needs new audio samples to play. My SDR is sampling at 240khz which gives me a nice even decimation factor of 10 to get to my audio rate. I'm processing 10,000 RF samples at a time in this callback, which means I'm giving the audio device 1000 samples at a time.

The getSamples() function is actually pretty large and contains some calls to UI objects to things like RSSI bars, but below is a slimmed-down version containing only the relevant demodulation code to get from baseband to NFM audio. The FM demod algorithm was borrowed from this post.

(sig is the scipy.signal module and np is numpy)

#########################################
#   Main Sample Processing Routine
#########################################

# Read samples
samples = self.sdr.read_samples(RadioSamplesPer)

# Decimate samples before audio processing
baseband = sig.decimate(samples,int(Decimation))

# FM Demod #
if self.mode == "NFM":

    # Filter RF baseband (lowpass at 5khz)
    filtered = Filters.LP_NFM_RF(baseband, AudioSampleRate)

    # Add in the last sample from the previous block for "history"
    baseband = np.insert(baseband, 0, self.lastSample)

    # Compute phi from the complex sample
    phi = np.arctan2(np.imag(baseband), np.real(baseband))
    
    # Demodulate by taking the difference between each sample
    outputSamples = np.diff(np.unwrap(phi)/(2*np.pi*0.1))

    # Get the new last sample
    self.lastSample = baseband[AudioSamplesPer]

    # Filter Audio (Bandpass at 300-3000khz)
    outputSamples = Filters.BP_NFM_AUD(outputSamples, AudioSampleRate)

    # Deemphasis Filter (RC filter with R=50 and C=4uF)
    outputSamples = Filters.LP_RC(outputSamples, 50, 4e-6, AudioSamplesPer)

    # Clip to (-1,1)
    outputSamples = np.clip(outputSamples,-1,1)

# Return the samples
return outputSamples

This seems to be generally working. I can get real-time demodulated audio, and it sounds correct, but it's not perfect. I can hear an audible click every time a new sample block is processed, which is both annoying and a sign that things aren't working quite right. Below is an annotated capture of a demodulated 1kHz tone coming from my demodulation routine (click for full size).

Issue 1 (red): Sample Jumps Between Blocks

Every other sample block, I'm seeing discontinuities at the transition point between blocks. Originally, I was seeing this issue between every block, but adding in the "history" for the FM demodulator seen in the above code cut that in half.

I'm assuming this issue stems from the demodulator not knowing what is coming in the next sample block. The fact that it got better when adding "history" seems to indicate I'm on the right track. But my non-DSP background isn't helping me to figure out where to go next.

Issue 2 (blue): Strange Output Envelope

This is the other strange issue causing the clicking between sample blocks. For some reason each sample block has a decaying output envelope instead of a constant envelope like I would expect with a solid 1kHz tone. Perhaps one of my filters is causing this unintentionally? It's a time-domain artifact, so I wouldn't think that would be the case.

Answer 1

I'm about 99.44% sure that you're filtering wrong, in that you're not carrying over any previous filter state. You've kind of fixed this in your differentiation step by saving the previous sample, because the filter that np.diff (effectively) implements uses the previous sample as it's current 'state'.

In your filter-and-decimate stage (I assume that's what Filters.LP_NFM_RF is doing) you need to properly initialize the filter. If it's a FIR filter you need to do a correct overlap-and-add (search on that term). If it's an IIR filter you need to retain the state of the filter from running the last set of data through it and use that as the state for the next set of data. If Filters is a package you didn't write, then hopefully they've implemented something to correctly splice filtered segments together; you need to figure out how they want you to use that -- and use it.

Answer 2

So Tim Wescott's answer was most likely a part of the issue, but during my attempts to correct my IIR filters I realized something even simpler was going on. My LP_RC filter class was calculating its output entirely incorrectly, and in a spectacularly bad way.

For some reason, my sleep-deprived mind tried to implement an RC filter by simply multiplying the input sample array by the filter's time constant.

def LP_RC(samples,R,C,length):
    f = np.arange(0,length,1)
    w = 2.0*np.pi*f
    return samples * 1/(w*R*C+1)

I've since corrected the my RC implenetation to something sensible, and added in initial/final conditions for my RF & audio filtering as well. My output is now buttery-smooth with no strange artifacting.