From d4e12c75cefc450e2f0480dfb730e28d83137088 Mon Sep 17 00:00:00 2001
From: Joris van Rantwijk <joris@jorisvr.nl>
Date: Sun, 5 Jan 2014 16:33:05 +0100
Subject: [PATCH] Add notes. Start building Python module with useful
 algorithms.

---
 NOTES.txt  |  17 +++++
 TODO.txt   |  12 ++++
 plltest.py | 103 -----------------------------
 pyfm.py    | 186 +++++++++++++++++++++++++++++++++++++++++++++++++++++
 4 files changed, 215 insertions(+), 103 deletions(-)
 create mode 100644 NOTES.txt
 delete mode 100644 plltest.py
 create mode 100644 pyfm.py

diff --git a/NOTES.txt b/NOTES.txt
new file mode 100644
index 0000000..78eeef5
--- /dev/null
+++ b/NOTES.txt
@@ -0,0 +1,17 @@
+
+This file contains random notitions
+-----------------------------------
+
+Sample rates between 301000 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.
+RTL-SDR currently configures it for cutoff at 1.2 MHz (2.4 MS/s).
+
+Casual test of ADC errors:
+ * 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
+
diff --git a/TODO.txt b/TODO.txt
index 8ab4322..e3af8f5 100644
--- a/TODO.txt
+++ b/TODO.txt
@@ -1,3 +1,15 @@
+* (experiment) measure raw signal for radio3, radio4 for ~ 1 minute
+   * different gain settings
+   * with/without RTL AGC mode
+   * with several IF gain settings
+* confirm theories about effect of gain, IF gain, AGC
+* look for effect of gain on baseband SNR
+* look for effect of ADC calibration on baseband SNR
+* look for effect of IF bandwidth on baseband SNR
+
+* (experiment) try if RTL2832 FIR filter can be optimized
+* (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
diff --git a/plltest.py b/plltest.py
deleted file mode 100644
index ecfd3d5..0000000
--- a/plltest.py
+++ /dev/null
@@ -1,103 +0,0 @@
-"""
-Test PLL algorithm.
-
-Usage: testpll.py baseband.dat centerfreq bandwidth > output.dat
-
-  baseband.dat   Raw 16-bit signed little-endian sample stream
-  centerfreq     Center frequency relative to sample frequency (0.5 = Nyquist)
-  bandwidth      Approximate bandwidth of PLL relative to sample frequency
-  output.dat     ASCII file with space-separated data
-"""
-
-import sys
-import numpy
-
-
-def pll(d, centerfreq, bandwidth):
-
-    minfreq = (centerfreq - bandwidth) * 2 * numpy.pi
-    maxfreq = (centerfreq + bandwidth) * 2 * numpy.pi
-
-    w = bandwidth * 2 * numpy.pi
-    phasor_a = numpy.poly([ numpy.exp(-1.146*w), numpy.exp(-5.331*w) ])
-    phasor_b = numpy.array([ sum(phasor_a) ])
-
-    loopfilter_b = numpy.poly([ numpy.exp(-0.1153*w) ])
-    loopfilter_b *= 0.62 * w
-
-    n = len(d)
-    y = numpy.zeros(n)
-    phasei = numpy.zeros(n)
-    phaseq = numpy.zeros(n)
-    phaseerr = numpy.zeros(n)
-    freq = numpy.zeros(n)
-    phase = numpy.zeros(n)
-    freq[0] = centerfreq * 2 * numpy.pi
-
-    phasor_i1 = phasor_i2 = 0
-    phasor_q1 = phasor_q2 = 0
-    loopfilter_x1 = 0
-
-    for i in xrange(n):
-
-        psin = numpy.sin(phase[i])
-        pcos = numpy.cos(phase[i])
-        y[i] = pcos
-
-        pi = pcos * d[i]
-        pq = psin * d[i]
-
-        pi = phasor_b[0] * pi - phasor_a[1] * phasor_i1 - phasor_a[2] * phasor_i2
-        pq = phasor_b[0] * pq - phasor_a[1] * phasor_q1 - phasor_a[2] * phasor_q2
-        phasor_i2 = phasor_i1
-        phasor_i1 = pi
-        phasor_q2 = phasor_q1
-        phasor_q1 = pq
-
-        phasei[i] = pi
-        phaseq[i] = pq
-
-        if pi > abs(pq):
-            perr = pq / pi
-        elif pq > 0:
-            perr = 1
-        else:
-            perr = -1
-        phaseerr[i] = perr
-
-        dfreq = loopfilter_b[0] * perr + loopfilter_b[1] * loopfilter_x1
-        loopfilter_x1 = perr
-
-        if i + 1 < n:
-            freq[i+1] = min(maxfreq, max(minfreq, freq[i] - dfreq))
-            p = phase[i] + freq[i+1]
-            if p > 2 * numpy.pi:  p -= 2 * numpy.pi
-            if p < -2 * numpy.pi: p += 2 * numpy.pi
-            phase[i+1] = p
-
-    return y, phasei, phaseq, phaseerr, freq, phase
-
-
-def main():
-
-    if len(sys.argv) != 4:
-        print >>sys.stderr, __doc__
-        sys.exit(1)
-
-    infile      = sys.argv[1]
-    centerfreq  = float(sys.argv[2])
-    bandwidth   = float(sys.argv[3])
-
-    d = numpy.fromfile(infile, '<i2')
-    d = d.astype(numpy.float64) / 32767.0
-
-    (y, phasei, phaseq, phaseerr, freq, phase) = pll(d, centerfreq, bandwidth)
-
-    print '#output phasei, phaseq, phaseerr freq phase'
-    for i in xrange(len(y)):
-        print y[i], phasei[i], phaseq[i], phaseerr[i], freq[i], phase[i]
-
-
-if __name__ == '__main__':
-    main()
-
diff --git a/pyfm.py b/pyfm.py
new file mode 100644
index 0000000..831be66
--- /dev/null
+++ b/pyfm.py
@@ -0,0 +1,186 @@
+"""
+Test lab for FM decoding algorithms.
+"""
+
+import sys
+import types
+import numpy
+import numpy.fft
+
+
+def readRawSamples(fname):
+    """Read raw sample file from rtl_sdr."""
+
+    d = numpy.fromfile(fname, dtype=numpy.uint8)
+    d = d.astype(numpy.float64)
+    d = (d - 128) / 128.0
+
+    return d[::2] - 1j * d[1::2]
+
+
+def lazyRawSamples(fname, blocklen):
+    """Return generator over blocks of raw samples."""
+
+    f = file(fname, 'rb')
+    while 1:
+        d = f.read(2 * blocklen)
+        if len(d) < 2 * blocklen:
+            break
+        d = numpy.fromstring(d, dtype=numpy.uint8)
+        d = d.astype(numpy.float64)
+        d = (d - 128) / 128.0
+        yield d[::2] - 1j * d[1::2]
+
+
+def freqShiftIQ(d, freqshift):
+    """Shift frequency by multiplication with complex phasor."""
+
+    def g(d, freqshift):
+        p = 0
+        for b in d:
+            n = len(b)
+            w = numpy.exp((numpy.arange(n) + p) * (2j * numpy.pi * freqshift))
+            p += n
+            yield b * w
+
+    if isinstance(d, types.GeneratorType):
+        return g(d, freqshift)
+    else:
+        n = len(d)
+        w = numpy.exp(numpy.arange(n) * (2j * numpy.pi * freqshift))
+        return d * w
+
+
+def spectrum(d, fs=1, nfft=None, sortfreq=False):
+    """Calculate Welch-style power spectral density.
+
+    fs       :: sample rate, default to 1
+    nfft     :: FFT length, default to block length
+    sortfreq :: True to put negative freqs in front of positive freqs
+
+    Use Hann window with 50% overlap.
+
+    Return (freq, Pxx)."""
+
+    if not isinstance(d, types.GeneratorType):
+        d = [ d ]
+
+    prev = None
+
+    if nfft is not None:
+        assert nfft > 0
+        w = numpy.hanning(nfft)
+        q = numpy.zeros(nfft)
+
+    pos = 0
+    i = 0
+    for b in d:
+
+        if nfft is None:
+            nfft = len(b)
+            assert nfft > 0
+            w = numpy.hanning(nfft)
+            q = numpy.zeros(nfft)
+
+        while pos + nfft <= len(b):
+
+            if pos < 0:
+                t = numpy.concatenate((prev[pos:], b[:pos+nfft]))
+            else:
+                t = b[pos:pos+nfft]
+
+            t *= w
+            tq = numpy.fft.fft(t)
+            tq *= numpy.conj(tq)
+            q += numpy.real(tq)
+
+            del t
+            del tq
+
+            pos += (nfft+(i%2)) // 2
+            i += 1
+
+        pos -= len(b)
+        if pos + len(b) > 0:
+            prev = b
+        else:
+            prev = numpy.concatenate((prev[pos+len(b):], b))
+
+    if i > 0:
+        q /= (i * numpy.sum(numpy.square(w)) * fs)
+
+    f = numpy.arange(nfft) * (fs / float(nfft))
+    f[nfft//2:] -= fs
+
+    if sortfreq:
+        f = numpy.concatenate((f[nfft//2:], f[:nfft//2]))
+        q = numpy.concatenate((q[nfft//2:], q[:nfft//2]))
+
+    return f, q
+
+
+def pll(d, centerfreq, bandwidth):
+
+    minfreq = (centerfreq - bandwidth) * 2 * numpy.pi
+    maxfreq = (centerfreq + bandwidth) * 2 * numpy.pi
+
+    w = bandwidth * 2 * numpy.pi
+    phasor_a = numpy.poly([ numpy.exp(-1.146*w), numpy.exp(-5.331*w) ])
+    phasor_b = numpy.array([ sum(phasor_a) ])
+
+    loopfilter_b = numpy.poly([ numpy.exp(-0.1153*w) ])
+    loopfilter_b *= 0.62 * w
+
+    n = len(d)
+    y = numpy.zeros(n)
+    phasei = numpy.zeros(n)
+    phaseq = numpy.zeros(n)
+    phaseerr = numpy.zeros(n)
+    freq = numpy.zeros(n)
+    phase = numpy.zeros(n)
+    freq[0] = centerfreq * 2 * numpy.pi
+
+    phasor_i1 = phasor_i2 = 0
+    phasor_q1 = phasor_q2 = 0
+    loopfilter_x1 = 0
+
+    for i in xrange(n):
+
+        psin = numpy.sin(phase[i])
+        pcos = numpy.cos(phase[i])
+        y[i] = pcos
+
+        pi = pcos * d[i]
+        pq = psin * d[i]
+
+        pi = phasor_b[0] * pi - phasor_a[1] * phasor_i1 - phasor_a[2] * phasor_i2
+        pq = phasor_b[0] * pq - phasor_a[1] * phasor_q1 - phasor_a[2] * phasor_q2
+        phasor_i2 = phasor_i1
+        phasor_i1 = pi
+        phasor_q2 = phasor_q1
+        phasor_q1 = pq
+
+        phasei[i] = pi
+        phaseq[i] = pq
+
+        if pi > abs(pq):
+            perr = pq / pi
+        elif pq > 0:
+            perr = 1
+        else:
+            perr = -1
+        phaseerr[i] = perr
+
+        dfreq = loopfilter_b[0] * perr + loopfilter_b[1] * loopfilter_x1
+        loopfilter_x1 = perr
+
+        if i + 1 < n:
+            freq[i+1] = min(maxfreq, max(minfreq, freq[i] - dfreq))
+            p = phase[i] + freq[i+1]
+            if p > 2 * numpy.pi:  p -= 2 * numpy.pi
+            if p < -2 * numpy.pi: p += 2 * numpy.pi
+            phase[i+1] = p
+
+    return y, phasei, phaseq, phaseerr, freq, phase
+
+