Reorganize notes; measure effect of LNA gain, ADC calibration.

Python code to experiment with IF filtering.

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
+Default settings in librtlsdr
-                    otherwise gain as configured via rtlsdr_set_tuner_gain
+-----------------------------
-  Elonics mixer gain: autonomous control disabled,
+
-                      gain depending on rtlsdr_set_tuner_gain
+Elonics LNA gain: when auto tuner gain: autonomous control with slow update
-  Elonics IF linearity: optimize sensitivity (default), auto switch disabled
+                  otherwise gain as configured via rtlsdr_set_tuner_gain
-  Elonics IF gain: +6, +0, +0, +0, +9, +9 (non-standard mode)
+Elonics mixer gain: autonomous control disabled,
-  Elonics IF filters: matched to sample rate (note this may not be optimal)
+                    gain depending on rtlsdr_set_tuner_gain
-  RTL AGC mode off
+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)
 
+--

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

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):
+def quadratureDetector(d, fs):
-    """FM frequency detector based on quadrature demodulation."""
+    """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:
+            if prev is not None:
-                yield numpy.angle(b[1:] * b[:-1].conj())
+                x = numpy.concatenate((prev[1:], b[:1]))
-            else:
+                yield numpy.angle(x * prev.conj()) * k
-                x = numpy.concatenate((prev[-1:], b[:-1]))
-                yield numpy.angle(b * x.conj())
             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
+