I've realised that I didn't implement the decimation step correctly when processing the output of my first receiver test.

In particular, I didn't low-pass filter the demodulated signal prior to downsampling, meaning that other signals present in the broadcast FM baseband signal (e.g. the stereo audio subcarrier) were aliased in to the mono audio channel I extracted.

I've fixed that, and reprocessed the IQ capture used in my previous post. The resulting audio is noticeably improved.

Here's the updated code:

import sys
import numpy as np
from scipy.signal import firwin
from scipy.io import wavfile

def pk(arr):
    return np.max(np.abs(arr))

def norm(chan):
    chan0 = chan - np.mean(chan)
    out = chan0 / pk(chan0)
    return out

def demod(xt):
    return np.diff(np.unwrap(np.angle(xt)))

def main(infile, outfile):
    channel = np.load(infile)

    sample_rate = 250e3
    M = 6 # decimate by this

    iq = norm(channel[0]) + 1j*norm(channel[1])
    x = demod(iq)
    # low pass filter
    taps = firwin(numtaps=101, cutoff=15e3, fs=sample_rate)
    x = np.convolve(x, taps, 'valid')

    x = x[::M] # decimate
    x /= pk(x) # normalize
    x *= 32767
    x = x.astype(np.int16)

    wavfile.write(outfile, int(sample_rate / M), x)

if __name__ == '__main__':
    if len(sys.argv) < 3:
        print(f'usage: {sys.argv[0]}  ')
    else:
        main(sys.argv[1], sys.argv[2])

As an aside, I note that this fix doesn't help with the two additional captures from which I'm yet to extract recognisable audio.