Working SDRPlay audio out - committing before I break something

demo.py

@@ -0,0 +1,83 @@
+#! /usr/bin/env python3
+
+#
+#   Example script for interfacing with SDRPlay radio.
+#
+# Modified from: https://github.com/pothosware/SoapySDR/wiki/PythonSupport
+#
+
+import SoapySDR
+from SoapySDR import * #SOAPY_SDR_ constants
+import numpy as np
+import sys
+import struct
+
+#enumerate devices
+devices = [dict(device) for device in SoapySDR.Device.enumerate()]
+
+if len(devices) == 0:
+    print('No SDR devices available.')
+    sys.exit(1)
+
+#create device instance
+sdr = SoapySDR.Device(devices[0])
+
+#apply settings
+sdr.setSampleRate(SOAPY_SDR_RX, 0, 1e6)
+sdr.setFrequency(SOAPY_SDR_RX, 0, 105500000)
+
+#setup a stream (complex floats)
+samples = 10000
+buff = np.array([0]*samples*2, np.int16)
+rxStream = sdr.setupStream(SOAPY_SDR_RX, SOAPY_SDR_CS16)
+sdr.activateStream(rxStream)
+
+
+#receive some samples
+for i in range(1000):
+    sr = sdr.readStream(rxStream, [buff], samples*2)
+
+    print(sr.ret, sr.flags, sr.timeNs, file=sys.stderr)
+
+    sys.stdout.buffer.write(struct.pack('=%dh' % len(buff), *buff))
+
+#shutdown the stream
+sdr.deactivateStream(rxStream) #stop streaming
+sdr.closeStream(rxStream)
+
+
+
+
+
+'''
+#create a re-usable buffer for rx samples
+# https://stackoverflow.com/a/32877245
+
+    #while True:
+    sr = sdr.readStream(rxStream, [buff], samples)
+    #print(sr.ret) #num samples or error code
+    #print(sr.flags) #flags set by receive operation
+    #print(sr.timeNs) #timestamp for receive buffer
+
+    # we now have a np buffer of [(re_16:1, im_16:1), re_16:2, im_16:2), ...]
+    # that needs to become a buffer of [re_8:1, im_8:1, re_8:2, im_8:2, ...]
+    #print(buff)
+    #output = buff.flatten('F')
+    #print(output)
+
+    #output = np.empty(len(buff) * 2, dtype=np.float32)
+    #output[0::2] = output_left
+    #output[1::2] = output_right
+    #output = output.astype(int)
+    #print(output)
+
+    #scale_factor = 127 / np.max(buff)
+
+    # scale the values using numpy.interp
+    #uint8_buffer = np.interp(buff, (buff.min(), buff.max()), (-127, 127)).astype(np.int8)
+
+    output = buff
+    #output /= (np.max(np.int16) / np.max(np.int8))
+    #output /= (32768 / 127)
+    #output = np.divide(output, 1).astype(np.int8)
+'''

example.py

@@ -10,11 +10,10 @@ import SoapySDR
 from SoapySDR import * #SOAPY_SDR_ constants
 import numpy #use numpy for buffers
 import sys
+import struct
 
 #enumerate devices
 devices = [dict(device) for device in SoapySDR.Device.enumerate()]
-for device in devices:
-    print(device)
 
 if len(devices) == 0:
     print('No SDR devices available.')
@@ -23,12 +22,6 @@ if len(devices) == 0:
 #create device instance
 sdr = SoapySDR.Device(devices[0])
 
-#query device info
-print(sdr.listAntennas(SOAPY_SDR_RX, 0))
-print(sdr.listGains(SOAPY_SDR_RX, 0))
-freqs = sdr.getFrequencyRange(SOAPY_SDR_RX, 0)
-for freqRange in freqs: print(freqRange)
-
 #apply settings
 sdr.setSampleRate(SOAPY_SDR_RX, 0, 1e6)
 sdr.setFrequency(SOAPY_SDR_RX, 0, 912.3e6)
@@ -38,14 +31,14 @@ rxStream = sdr.setupStream(SOAPY_SDR_RX, SOAPY_SDR_CF32)
 sdr.activateStream(rxStream) #start streaming
 
 #create a re-usable buffer for rx samples
-buff = numpy.array([0]*1024, numpy.complex64)
+samples = 1000
+buff = numpy.array([0]*samples*2, numpy.float32)
 
 #receive some samples
 for i in range(10):
-    sr = sdr.readStream(rxStream, [buff], len(buff))
-    print(sr.ret) #num samples or error code
-    print(sr.flags) #flags set by receive operation
-    print(sr.timeNs) #timestamp for receive buffer
+    sr = sdr.readStream(rxStream, [buff], samples)
+
+    sys.stdout.buffer.write(struct.pack('=%df' % len(buff), *buff))
 
 #shutdown the stream
 sdr.deactivateStream(rxStream) #stop streaming

fm1s.py

@@ -18,10 +18,8 @@ if optimized:
 	import fastfm # Cython
 
 MAX_DEVIATION = 200000.0 # Hz
-INPUT_RATE = 256000
-OUTPUT_RATE = 32000
-#INPUT_RATE = 1e6
-#OUTPUT_RATE = 50000
+INPUT_RATE = 1000000
+OUTPUT_RATE = 50000
 
 if debug_mode:
 	OUTPUT_RATE=256000
@@ -67,7 +65,7 @@ while True:
 	data = sys.stdin.buffer.read((INPUT_RATE * 2) // 10)
 	if not data:
 		break
-	data = remaining_data + data
+	#data = remaining_data + data
 
 	if len(data) < 4:
 		remaining_data = data
@@ -85,12 +83,16 @@ while True:
 
 	# Find angle (phase) of I/Q pairs
 	iqdata = numpy.frombuffer(data, dtype=numpy.uint8)
+	print(iqdata.dtype, iqdata.shape, iqdata, file=sys.stderr)
 	iqdata = iqdata - 127.5
 	iqdata = iqdata / 128.0
+	print(iqdata.dtype, iqdata.shape, iqdata, file=sys.stderr)
 	iqdata = iqdata.view(complex)
+	print(iqdata.dtype, iqdata.shape, iqdata, file=sys.stderr)
 
 	angles = numpy.angle(iqdata)
 
+
 	# Determine phase rotation between samples
 	# (Output one element less, that's we always save last sample
 	# in remaining_data)

fm1s_int16.py

@@ -0,0 +1,225 @@
+#!/usr/bin/env python3
+# FM demodulator based on I/Q (quadrature) samples
+
+#
+# This code belongs to https://github.com/elvis-epx/sdr
+# No licence information is provided.
+#
+
+import struct, math, random, sys, numpy, filters, time
+
+optimized = "-o" in sys.argv
+debug_mode = "-d" in sys.argv
+disable_pll = "--disable-pll" in sys.argv
+if disable_pll:
+	optimized = False
+
+if optimized:
+	import fastfm # Cython
+
+MAX_DEVIATION = 200000.0 # Hz
+INPUT_RATE = 1000000
+OUTPUT_RATE = 50000
+
+if debug_mode:
+	OUTPUT_RATE=256000
+
+DECIMATION = INPUT_RATE / OUTPUT_RATE
+assert DECIMATION == math.floor(DECIMATION)
+
+FM_BANDWIDTH = 15000 # Hz
+STEREO_CARRIER = 38000 # Hz
+DEVIATION_X_SIGNAL = 0.999 / (math.pi * MAX_DEVIATION / (INPUT_RATE / 2))
+
+pll = math.pi - random.random() * 2 * math.pi
+last_pilot = 0.0
+deviation_avg = math.pi - random.random() * 2 * math.pi
+last_deviation_avg = deviation_avg
+tau = 2 * math.pi
+
+# Downsample mono audio
+decimate1 = filters.decimator(DECIMATION)
+
+# Deemph + Low-pass filter for mono (L+R) audio
+lo = filters.deemphasis(INPUT_RATE, 75, FM_BANDWIDTH, 120)
+
+# Downsample jstereo audio
+decimate2 = filters.decimator(DECIMATION)
+
+# Deemph + Low-pass filter for joint-stereo demodulated audio (L-R)
+lo_r = filters.deemphasis(INPUT_RATE, 75, FM_BANDWIDTH, 120)
+
+# Band-pass filter for stereo (L-R) modulated audio
+hi = filters.band_pass(INPUT_RATE,
+	STEREO_CARRIER - FM_BANDWIDTH, STEREO_CARRIER + FM_BANDWIDTH, 120)
+
+# Filter to extract pilot signal
+pilot = filters.band_pass(INPUT_RATE,
+	STEREO_CARRIER / 2 - 100, STEREO_CARRIER / 2 + 100, 120)
+
+last_angle = 0.0
+remaining_data = b''
+
+while True:
+	# Ingest 0.1s worth of data
+	bytes_per_sample = 1
+	data = sys.stdin.buffer.read((INPUT_RATE * bytes_per_sample * 2) // 10)
+	if not data:
+		break
+	#data = remaining_data + data
+
+	if len(data) < 4:
+		remaining_data = data
+		continue
+
+	# Save one sample to next batch, and the odd byte if exists
+	if len(data) % 2 == 1:
+		print("Odd byte, that's odd", file=sys.stderr)
+		remaining_data = data[-3:]
+		data = data[:-1]
+	else:
+		remaining_data = data[-2:]
+
+	samples = len(data) // (2 * bytes_per_sample)
+
+	# Find angle (phase) of I/Q pairs
+	iqdata = numpy.frombuffer(data, dtype=numpy.int16)
+	print(iqdata.dtype, iqdata.shape, iqdata, file=sys.stderr)
+	iqdata = iqdata / (2**15)
+	print(iqdata.dtype, iqdata.shape, iqdata, file=sys.stderr)
+	iqdata = iqdata.view(complex)
+	print(iqdata.dtype, iqdata.shape, iqdata, file=sys.stderr)
+
+	angles = numpy.angle(iqdata)
+
+
+	# Determine phase rotation between samples
+	# (Output one element less, that's we always save last sample
+	# in remaining_data)
+	rotations = numpy.ediff1d(angles)
+
+	# Wrap rotations >= +/-180º
+	rotations = (rotations + numpy.pi) % (2 * numpy.pi) - numpy.pi
+
+	# Convert rotations to baseband signal
+	output_raw = numpy.multiply(rotations, DEVIATION_X_SIGNAL)
+	output_raw = numpy.clip(output_raw, -0.999, +0.999)
+
+	# At this point, output_raw contains two audio signals:
+	# L+R (mono-compatible) and L-R (joint-stereo) modulated in AM-SC,
+	# carrier 38kHz
+
+	# Downsample and low-pass L+R (mono) signal
+	output_mono = lo.feed(output_raw)
+	output_mono = decimate1.feed(output_mono)
+
+	# Filter pilot tone
+	detected_pilot = pilot.feed(output_raw)
+
+	# Separate ultrasonic L-R signal by high-pass filtering
+	output_jstereo_mod = hi.feed(output_raw)
+
+	# Demodulate L-R, which is AM-SC with 53kHz carrier
+
+	if optimized:
+		output_jstereo, pll, STEREO_CARRIER, \
+		last_pilot, deviation_avg, last_deviation_avg = \
+			fastfm.demod_stereo(output_jstereo_mod,
+						pll,
+						STEREO_CARRIER,
+						INPUT_RATE,
+						detected_pilot,
+						last_pilot,
+						deviation_avg,
+						last_deviation_avg)
+	else:
+		output_jstereo = []
+
+		for n in range(0, len(output_jstereo_mod)):
+			# Advance carrier
+			pll = (pll + tau * STEREO_CARRIER / INPUT_RATE) % tau
+
+			# Standard demodulation
+			output_jstereo.append(math.cos(pll) * output_jstereo_mod[n])
+
+			if disable_pll:
+				continue
+
+			############ Carrier PLL #################
+
+			# Detect pilot zero-crossing
+			cur_pilot = detected_pilot[n]
+			zero_crossed = (cur_pilot * last_pilot) <= 0
+			last_pilot = cur_pilot
+			if not zero_crossed:
+				continue
+
+			# When pilot is at 90º or 270º, carrier should be around 180º
+			# t=0    => cos(t) = 1,  cos(2t) = 1
+			# t=π/2  => cos(t) = 0,  cos(2t) = -1
+			# t=π    => cos(t) = -1, cos(2t) = 1
+			# t=-π/2 => cos(t) = 0,  cos(2t) = -1
+			ideal = math.pi
+			deviation = pll - ideal
+			if deviation > math.pi:
+				# 350º => -10º
+				deviation -= tau
+			deviation_avg = 0.99 * deviation_avg + 0.01 * deviation
+			rotation = deviation_avg - last_deviation_avg
+			last_deviation_avg = deviation_avg
+
+			if abs(deviation_avg) > math.pi / 8:
+				# big phase deviation, reset PLL
+				# print("Resetting PLL", file=sys.stderr)
+				pll = ideal
+				pll = (pll + tau * STEREO_CARRIER / INPUT_RATE) % tau
+				deviation_avg = 0.0
+				last_deviation_avg = 0.0
+
+			# Translate rotation to frequency deviation e.g.
+			# cos(tau + 3.6º) = cos(1.01 * tau)
+			# cos(tau - 9º) = cos(tau * 0.975)
+                        #
+			# Overcorrect by 5% to (try to) sync phase,
+                        # otherwise only the frequency would be synced.
+
+			STEREO_CARRIER /= (1 + (rotation * 1.05) / tau)
+
+			'''
+			print("%d deviationavg=%f rotation=%f freq=%f" %
+				(n,
+				deviation_avg * 180 / math.pi,
+				rotation * 180 / math.pi,
+				STEREO_CARRIER),
+				file=sys.stderr)
+			time.sleep(0.05)
+			'''
+
+	# Downsample, Low-pass/deemphasis demodulated L-R
+	output_jstereo = lo_r.feed(output_jstereo)
+	output_jstereo = decimate2.feed(output_jstereo)
+
+	assert len(output_jstereo) == len(output_mono)
+
+	# Scale to 16-bit and divide by 2 for channel sum
+	output_mono = numpy.multiply(output_mono, 32767 / 2.0)
+	output_jstereo = numpy.multiply(output_jstereo, 32767 / 2.0)
+
+	# Output stereo by adding or subtracting joint-stereo to mono
+	output_left = output_mono + output_jstereo
+	output_right = output_mono - output_jstereo
+
+	if not debug_mode:
+		# Interleave L and R samples using NumPy trickery
+		output = numpy.empty(len(output_mono) * 2, dtype=output_mono.dtype)
+		output[0::2] = output_left
+		output[1::2] = output_right
+		output = output.astype(int)
+	else:
+		output = numpy.empty(len(output_mono) * 3, dtype=output_mono.dtype)
+		output[0::3] = output_mono
+		output[1::3] = output_jstereo
+		output[2::3] = numpy.multiply(detected_pilot, 32767)
+		output = output.astype(int)
+
+	sys.stdout.buffer.write(struct.pack('<%dh' % len(output), *output))

fm_rt_stereo

@@ -1,6 +1,4 @@
-#!/bin/sh -x
-
-# Decode and play stereo broadcast FM in realtime.
+#!/bin/sh
 
 #
 # This code belongs to https://github.com/elvis-epx/sdr
@@ -8,4 +6,16 @@
 #
 
 
-rtl_sdr -f 106.3M -s 256k - | ./fm1s.py | sox -t raw -r 32000 -b 16 -c 2 -L -e signed-integer - -d
+#./example.py | ./fm1s.py | sox -t raw -r 50000 -b 16 -c 2 -L -e signed-integer - -d
+
+#./example.py > buffer.iq
+#./demo.py > sdrplay_buffer.iq
+./fm1s_2.py < sdrplay_buffer.iq
+
+
+# Run the radio direct to audio. Something about 1M samples and 50000 output rate sounds ass, but it works
+#rtl_sdr -f 105.5M -s 1M - | ./fm1s.py | sox -t raw -r 50000 -b 16 -c 2 -L -e signed-integer - -d
+
+# Get data from the nooelec device
+#rtl_sdr -f 105.5M -s 1M -n 10000 nooelec_buffer.iq
+./fm1s.py < nooelec_buffer.iq

nooelec

@@ -0,0 +1,7 @@
+#! /bin/bash
+
+if ! [ -f nooelec_buffer.iq ]; then
+rtl_sdr -f 105.5M -s1M -n 1000000 nooelec_buffer.iq
+fi
+./fm1s.py < nooelec_buffer.iq > nooelec_buffer.raw
+./play.sh < nooelec_buffer.raw

play.sh

@@ -0,0 +1,3 @@
+#! /bin/bash
+
+sox -t raw -r 50000 -b 16 -c 2 -L -e signed-integer - -d

sdrplay

@@ -0,0 +1,7 @@
+#! /bin/bash
+
+if ! [ -f sdrplay_int16.iq ]; then
+./demo.py > sdrplay_int16.iq
+fi
+./fm1s_int16.py < sdrplay_int16.iq > sdrplay_int16.raw
+./play.sh < sdrplay_int16.raw

setup

@@ -1,5 +1,5 @@
 #! /bin/bash
 
-sudo apt install rtl-sdr portaudio19-dev python3-all-dev
+sudo apt install rtl-sdr sox portaudio19-dev python3-all-dev
 
 pip install pyaudio filters