diff --git a/NOTES.txt b/NOTES.txt index 3b14908..dead55b 100644 --- a/NOTES.txt +++ b/NOTES.txt @@ -2,15 +2,23 @@ This file contains random notitions ----------------------------------- + +Valid sample rates +------------------ + Sample rates between 300001 Hz and 900000 Hz (inclusive) are not supported. They cause an invalid configuration of the RTL chip. rsamp_ratio = 28.8 MHz * 2**22 / sample_rate If bit 27 and bit 28 of rsamp_ratio are different, the RTL chip malfunctions. -The RTL chip has a configurable 32-tap FIR filter. + +Behaviour of RTL and Elonics tuner +---------------------------------- + +The RTL chip has a configurable 32-tap FIR filter running at 28.8 MS/s. RTL-SDR currently configures it for cutoff at 1.2 MHz (2.4 MS/s). -Casual test of ADC errors: +Casual test of ADC mismatch: * DC offset in order of 1 code step * I/Q gain mismatch in order of 4% * I/Q phase mismatch in order of 1% of sample interval @@ -32,18 +40,54 @@ mode causes a brief level spike, while manually rewriting the same IF gain in AGC mode does not have any effect). It seems more likely that AGC is a digital gain in the downsampling filter. -Default settings in librtlsdr: - Elonics LNA gain: when auto tuner gain: autonomous control with slow update - otherwise gain as configured via rtlsdr_set_tuner_gain - Elonics mixer gain: autonomous control disabled, - gain depending on rtlsdr_set_tuner_gain - Elonics IF linearity: optimize sensitivity (default), auto switch disabled - Elonics IF gain: +6, +0, +0, +0, +9, +9 (non-standard mode) - Elonics IF filters: matched to sample rate (note this may not be optimal) - RTL AGC mode off + +Default settings in librtlsdr +----------------------------- + +Elonics LNA gain: when auto tuner gain: autonomous control with slow update + otherwise gain as configured via rtlsdr_set_tuner_gain +Elonics mixer gain: autonomous control disabled, + gain depending on rtlsdr_set_tuner_gain +Elonics IF linearity: optimize sensitivity (default), auto switch disabled +Elonics IF gain: +6, +0, +0, +0, +9, +9 (non-standard mode) +Elonics IF filters: matched to sample rate (note this may not be optimal) +RTL AGC mode off -Local radio stations: +Effect of settings on baseband SNR +---------------------------------- + +STATION SRATE LNA IF GAIN AGC SOFT BW IF LEVEL GUARD/PILOT + +radio3 1 MS/s 24 dB default off 150 kHz 0.19 -62.6 dB/Hz +radio3 1.5 MS 24 dB default off 150 kHz 0.19 -62.7 dB/Hz +radio3 2 MS/s 24 dB default off 150 kHz 0.18 -62.7 dB/Hz + +radio3 2 MS/s 34 dB default off 150 kHz 0.46 -64.0 dB/Hz +radio3 2 MS/s 34 dB default off 80 kHz -64.0 dB/Hz +radio3 2 MS/s 34 dB default off 150 kHz adccal -64.0 dB/Hz + +radio4 2 MS/s 24 dB default off 150 kHz 0.04 -41.1 dB/Hz + +radio4 1 MS/s 34 dB default off 150 kHz 0.06 -43.3 dB/Hz +radio4 1 MS/s 34 dB default off 80 kHz -51.2 dB/Hz + +radio4 2 MS/s 34 dB default off 150 kHz 0.10 -42.4 dB/Hz +radio4 2 MS/s 34 dB default off 80 kHz -48.2 dB/Hz +radio4 2 MS/s 34 dB default off 150 kHz adccal -42.4 dB/Hz + +Note: all measurements 10 second duration +Note: all measurements have LO frequency set to station + 250 kHz + +Conclusion: Sample rate (1 MS/s to 2 MS/s) has little effect on quality. +Conclusion: LNA gain has little effect on quality. +Conclusion: Narrow IF filter improves quality of weak station. +Conclusion: ADC gain/offset calibration has no effect on quality. + + +Local radio stations +-------------------- + radio2 92600000 (good) radio3 96800000 (good) radio4 94300000 (bad) @@ -52,3 +96,4 @@ radio538 102100000 (medium) radio10 103800000 (bad) radio west 89300000 (medium) +-- diff --git a/TODO.txt b/TODO.txt index 7ee9a50..1d56db0 100644 --- a/TODO.txt +++ b/TODO.txt @@ -1,19 +1,16 @@ +* (experiment) make nice plot of baseband distortion due to IF filtering +* (experiment) consider downsampling IF signal before FM detection * (experiment) measure effect of IF gain on baseband SNR * (experiment) measure effect of IF gain linearity on baseband SNR * (experiment) measure effect of RTL AGC mode on baseband SNR -* (experiment) measure effect of ADC calibration on baseband SNR -* (experiment) measure effect of IF bandwidth on baseband SNR -* (experiment) measure effect of IF sample rate on baseband SNR * (experiment) try if RTL AGC mode improves FM decoding * (feature) support 'M' 'k' suffixes for sample rates and tuning frequency * (feature) implement off-line FM decoder in Python for experimentation * (feature) implement stereo pilot pulse-per-second -* (quality) consider DC offset calibration * (speedup) maybe replace high-order FIR downsampling filter with 2nd order butterworth followed by lower order FIR filter * figure out why we sometimes lose stereo lock * it looks like IF level sometimes varies so much that it saturates the receiver; perhaps this can be solved by dynamically managing the hardware gain in response to level measurements -* (quality) figure out if I/Q balance can improve weak stations * (quality) figure out if hardware gain settings can improve weak stations * (feature) implement RDS decoding * (quality) consider FM demodulation with PLL instead of phase discriminator diff --git a/pyfm.py b/pyfm.py index d2e9c30..dfa2069 100644 --- a/pyfm.py +++ b/pyfm.py @@ -6,6 +6,8 @@ import sys import types import numpy import numpy.fft +import numpy.linalg +import numpy.random import scipy.signal @@ -78,24 +80,37 @@ def firFilter(d, coeff): return scipy.signal.lfilter(coeff, 1, d) -def quadratureDetector(d): - """FM frequency detector based on quadrature demodulation.""" +def quadratureDetector(d, fs): + """FM frequency detector based on quadrature demodulation. + Return an array of real-valued numbers, representing frequencies in Hz.""" + + k = fs / (2 * numpy.pi) # lazy version def g(d): prev = None for b in d: - if prev is None: - yield numpy.angle(b[1:] * b[:-1].conj()) - else: - x = numpy.concatenate((prev[-1:], b[:-1])) - yield numpy.angle(b * x.conj()) + if prev is not None: + x = numpy.concatenate((prev[1:], b[:1])) + yield numpy.angle(x * prev.conj()) * k prev = b + yield numpy.angle(prev[1:] * prev[:-1].conj()) * k if isinstance(d, types.GeneratorType): return g(d) else: - return numpy.angle(d[1:] * d[:-1].conj()) + return numpy.angle(d[1:] * d[:-1].conj()) * k + + +def modulateFm(sig, fs, fcenter=0): + """Create an FM modulated IQ signal. + + sig :: modulation signal, values in Hz + fs :: sample rate in Hz + fcenter :: center frequency in Hz + """ + + return numpy.exp(2j * numpy.pi * (sig + fcenter).cumsum() / fs) def spectrum(d, fs=1, nfft=None, sortfreq=False): @@ -232,11 +247,12 @@ def pll(d, centerfreq, bandwidth): return y, phasei, phaseq, phaseerr, freq, phase -def pilotLevel(d, fs, freqshift, bw=150.0e3): +def pilotLevel(d, fs, freqshift, nfft=None, bw=150.0e3): """Calculate level of the 19 kHz pilot vs noise floor in the guard band. d :: block of raw I/Q samples or lazy I/Q sample stream fs :: sample frequency in Hz + nfft :: FFT length freqshift :: frequency offset in Hz bw :: half-bandwidth of IF signal in Hz @@ -255,10 +271,10 @@ def pilotLevel(d, fs, freqshift, bw=150.0e3): d = firFilter(d, b) # Demodulate FM. - d = quadratureDetector(d) + d = quadratureDetector(d, fs) # Power spectral density. - f, q = spectrum(d, fs=fs, sortfreq=False) + f, q = spectrum(d, fs=fs, nfft=nfft, sortfreq=False) # Locate 19 kHz bin. k19 = int(19.0e3 * len(q) / fs) @@ -266,7 +282,7 @@ def pilotLevel(d, fs, freqshift, bw=150.0e3): k19 = k19 - kw + numpy.argmax(q[k19-kw:k19+kw]) # Calculate pilot power. - p19 = numpy.sum(q[k19-1:k19+2]) * fs * 0.75 / len(q) + p19 = numpy.sum(q[k19-1:k19+2]) * fs * 1.5 / len(q) # Calculate noise floor in guard band. k17 = int(17.0e3 * len(q) / fs) @@ -278,3 +294,66 @@ def pilotLevel(d, fs, freqshift, bw=150.0e3): return (p19db, guarddb, guarddb - p19db) + +def modulateAndReconstruct(sigfreq, sigampl, nsampl, fs, noisebw=None, ifbw=None, ifnoise=0): + """Create a pure sine wave, modulate to FM, add noise, filter, demodulate. + + sigfreq :: frequency of sine wave in Hz + sigampl :: amplitude of sine wave in Hz (carrier swing) + nsampl :: number of samples + fs :: sample rate in Hz + noisebw :: calculate noise after demodulation over this bandwidth + ifbw :: IF filter bandwidth in Hz, or None for no filtering + ifnoise :: IF noise level + + Return (ampl, phase, noise) + where ampl is the amplitude of the reconstructed sine wave (~ sigampl) + phase is the phase shift after reconstruction + noise is the standard deviation of noise in the reconstructed signal + """ + + # Make sine wave. + sig0 = sigampl * numpy.sin(2*numpy.pi*sigfreq/fs * numpy.arange(nsampl)) + + # Modulate to IF. + fm = modulateFm(sig0, fs=fs, fcenter=0) + + # Add noise. + if ifnoise: + fm += numpy.sqrt(0.5) * numpy.random.normal(0, ifnoise, nsampl) + fm += 1j * numpy.sqrt(0.5) * numpy.random.normal(0, ifnoise, nsampl) + + # Filter IF. + if ifbw is not None: + b = scipy.signal.firwin(61, 2.0 * ifbw / fs, window='nuttall') + fm = scipy.signal.lfilter(b, 1, fm) + fm = fm[61:] + + # Demodulate. + sig1 = quadratureDetector(fm, fs=fs) + + # Fit original sine wave. + k = len(sig1) + m = numpy.zeros((k, 3)) + m[:,0] = numpy.sin(2*numpy.pi*sigfreq/fs * (numpy.arange(k) + nsampl - k)) + m[:,1] = numpy.cos(2*numpy.pi*sigfreq/fs * (numpy.arange(k) + nsampl - k)) + m[:,2] = 1 + fit = numpy.linalg.lstsq(m, sig1) + csin, ccos, coffset = fit[0] + del fit + + # Calculate amplitude, phase. + ampl1 = numpy.sqrt(csin**2 + ccos**2) + phase1 = numpy.arctan2(-ccos, csin) + + # Calculate residual noise. + res1 = sig1 - m[:,0] * csin - m[:,1] * ccos + + if noisebw is not None: + b = scipy.signal.firwin(61, 2.0 * noisebw / fs, window='nuttall') + res1 = scipy.signal.lfilter(b, 1, res1) + + noise1 = numpy.sqrt(numpy.mean(res1 ** 2)) + + return ampl1, phase1, noise1 +