17. Ejemplo Punto a Punto

En este cap铆tulo reunimos muchos de los conceptos que aprendimos anteriormente y analizamos un ejemplo completo de recepci贸n y decodificaci贸n de una se帽al digital real. Analizaremos el Radio Data System (RDS), que es un protocolo de comunicaciones para incorporar peque帽as cantidades de informaci贸n en transmisiones de radio FM, como la estaci贸n y el nombre de la canci贸n. Tendremos que demodular FM, cambiar de frecuencia, filtrar, diezmar, remuestrear, sincronizar, decodificar y analizar los bytes. Se proporciona un archivo IQ de ejemplo para fines de prueba o si no tiene un SDR a mano.

Introducci贸n a la radio FM y RDS

Para entender el RDS primero debemos repasar las emisiones de radio FM y c贸mo estan estructuradas estas se帽ales. Probablemente est茅 familiarizado con la parte de audio de las se帽ales de FM, que son simplemente se帽ales de audio moduladas en frecuencia y transmitidas en frecuencias centrales correspondientes al nombre de la estaci贸n, por ejemplo, 鈥淲PGC 95.5 FM鈥 est谩 centrada exactamente en 95,5 MHz. Adem谩s de la parte de audio, cada transmisi贸n de FM contiene otros componentes que se modulan en frecuencia junto con el audio. En lugar de simplemente buscar en Google la estructura de la se帽al, echemos un vistazo a la densidad espectral de potencia (PSD) de una se帽al de FM de ejemplo, despu茅s de la demodulaci贸n de FM. Solo vemos la parte positiva porque la salida de la demodulaci贸n de FM es una se帽al real, aunque la entrada sea compleja (veremos el c贸digo para realizar esta demodulaci贸n en breve).

Power spectral density (PSD) of an FM radio signal after the FM demodulation, showing RDS

Al observar la se帽al en el dominio de la frecuencia, notamos las siguientes se帽ales individuales:

  1. Una se帽al de alta potencia entre 0 - 17 kHz

  2. Un tono a 19 kHz

  3. Centrada en 38 kHz y aproximadamente 30 kHz de ancho, vemos una se帽al sim茅trica de aspecto interesante.

  4. Se帽al en forma de doble l贸bulo centrada en 57 kHz

  5. Se帽al en forma de un solo l贸bulo centrada en 67 kHz

Esto es esencialmente todo lo que podemos determinar con solo mirar el PSD, y recuerde que esto es despu茅s de la demodulaci贸n de FM. El PSD antes de la demodulaci贸n de FM se parece a lo siguiente, lo que realmente no nos dice mucho.

Power spectral density (PSD) of an FM radio signal before any demodulation

Dicho esto, es importante entender que cuando modulas una se帽al en FM, una frecuencia m谩s alta en la se帽al de datos conducir谩 a una frecuencia m谩s alta en la se帽al de FM resultante. Por lo tanto, la presencia de una se帽al centrada en 67 kHz aumenta el ancho de banda total ocupado por la se帽al de FM transmitida, ya que el componente de frecuencia m谩xima ahora es de alrededor de 75 kHz, como se muestra en el primer PSD anterior. La reglas del ancho de banda de Carson aplicado a FM nos dice que las estaciones de FM ocupan aproximadamente 250 kHz de espectro, raz贸n por la cual generalmente tomamos muestras a 250 kHz (recuerde que cuando usamos muestreo en cuadratura/IQ, el ancho de banda recibido es igual a su frecuencia de muestreo).

Como comentario breve, algunos lectores pueden estar familiarizados con observar la banda de FM usando un SDR o un analizador de espectro y ver su espectrograma, y pensar que las se帽ales de bloque y adyacentes a algunas de las estaciones de FM son RDS.

Spectrogram of the FM band

Resulta que esas se帽ales en bloque son en realidad HD Radio, una versi贸n digital de la misma se帽al de radio FM (mismo contenido de audio). Esta versi贸n digital genera una se帽al de audio de mayor calidad en el receptor porque la FM anal贸gica siempre incluir谩 algo de ruido despu茅s de la demodulaci贸n, ya que es un esquema anal贸gico, pero la se帽al digital se puede demodular/decodificar con cero ruido, suponiendo que no haya errores de bits.

Volvamos a las cinco se帽ales que descubrimos en nuestro PSD; el siguiente diagrama se etiqueta para qu茅 se utiliza cada se帽al.

Components within an FM radio signal, including mono and stereo audio, RDS, and DirectBand signals

Repasando cada una de estas se帽ales sin ning煤n orden en particular:

Las se帽ales de audio mono y est茅reo simplemente transportan la se帽al de audio, en un patr贸n en el que sumarlas y restarlas da como resultado los canales izquierdo y derecho.

El tono piloto de 19 kHz se utiliza para demodular el audio est茅reo. Si duplicas el tono act煤a como referencia de frecuencia y fase, ya que la se帽al de audio est茅reo est谩 centrada en 38 kHz. Se puede duplicar el tono simplemente elevando al cuadrado las muestras; recuerde la propiedad de Fourier de cambio de frecuencia que aprendimos en el capitulo Dominio de la Frecuencia .

DirectBand era una red de transmisi贸n de datos inal谩mbrica de Am茅rica del Norte propiedad de Microsoft y operada por ella, tambi茅n llamada 鈥淢SN Direct鈥 en los mercados de consumo. DirectBand transmit铆a informaci贸n a dispositivos como receptores GPS port谩tiles, relojes de pulsera y estaciones meteorol贸gicas dom茅sticas. Incluso permit铆a a los usuarios recibir mensajes cortos desde Windows Live Messenger. Una de las aplicaciones m谩s exitosas de DirectBand fueron los datos de tr谩fico local en tiempo real mostrados en los receptores GPS Garmin, que eran utilizados por millones de personas antes de que los tel茅fonos inteligentes se volvieran omnipresentes. El servicio DirectBand se cerr贸 en enero de 2012, lo que plantea la pregunta: 驴por qu茅 lo vemos en nuestra se帽al de FM grabada despu茅s de 2012? Mi 煤nica suposici贸n es que la mayor铆a de los transmisores de FM fueron dise帽ados y construidos mucho antes de 2012, e incluso sin ninguna 鈥渟uministro鈥 activo por DirectBand, todav铆a transmite algo, tal vez s铆mbolos piloto.

Por 煤ltimo, llegamos a RDS, que es el foco del resto de este cap铆tulo. Como podemos ver en nuestro primer PSD, RDS tiene aproximadamente 4 kHz de ancho de banda (antes de ser modulado en FM) y se encuentra entre el audio est茅reo y la se帽al DirectBand. Es un protocolo de comunicaciones digitales de baja velocidad de datos que permite a las estaciones de FM incluir identificaci贸n de la estaci贸n, informaci贸n del programa, hora y otra informaci贸n diversa junto con el audio. El est谩ndar RDS se publica como est谩ndar IEC 62106 y se puede encontrar aqui.

La se帽al RDS

En este cap铆tulo usaremos Python para recibir RDS, pero para comprender mejor c贸mo recibirlo, primero debemos aprender c贸mo se forma y transmite la se帽al.

Lado de transmisi贸n

La informaci贸n RDS que transmitir谩 la estaci贸n de FM (por ejemplo, nombre de pista, etc.) est谩 codificada en conjuntos de 8 bytes. Cada conjunto de 8 bytes, que corresponde a 64 bits, se combina con 40 鈥渂its de verificaci贸n鈥 para formar un 煤nico 鈥済rupo鈥. Estos 104 bits se transmiten juntos, aunque no hay un intervalo de tiempo entre los grupos, por lo que desde la perspectiva del receptor recibe estos bits sin parar y debe determinar el l铆mite entre los grupos de 104 bits. Veremos m谩s detalles sobre la codificaci贸n y la estructura del mensaje una vez que nos sumerjamos en el lado de la recepci贸n.

Para transmitir estos bits de forma inal谩mbrica, RDS utiliza BPSK, que como aprendimos en el capitulo Modulaci贸n Digital es un esquema de modulaci贸n digital simple que se utiliza para asignar unos y ceros a la fase de una portadora. Como muchos protocolos basados en BPSK, RDS utiliza codificaci贸n diferencial, lo que simplemente significa que los 1 y 0 de los datos se codifican en cambios de 1 y 0, lo que le permite ya no preocuparse si est谩 desfasado 180 grados (m谩s sobre esto m谩s adelante). Los s铆mbolos BPSK se transmiten a 1187,5 s铆mbolos por segundo y, debido a que BPSK transporta un bit por s铆mbolo, eso significa que RDS tiene una velocidad de datos sin procesar de aproximadamente 1,2 kbps (incluida la sobrecarga). RDS no contiene ninguna codificaci贸n de canal (tambi茅n conocida como correcci贸n de errores directa), aunque los paquetes de datos contienen una verificaci贸n de redundancia c铆clica (CRC) para saber cu谩ndo ocurri贸 un error.

Luego, la se帽al BPSK final se cambia de frecuencia hasta 57 kHz y se agrega a todos los dem谩s componentes de la se帽al de FM, antes de ser modulada en FM y transmitida por aire en la frecuencia de la estaci贸n. Las se帽ales de radio FM se transmiten a una potencia extremadamente alta en comparaci贸n con la mayor铆a de las dem谩s comunicaciones inal谩mbricas, 隆hasta 80 kW! Esta es la raz贸n por la que muchos usuarios de SDR tienen un filtro de rechazo de FM (es decir, un filtro de eliminaci贸n de banda) en l铆nea con su antena; por lo que FM no a帽ade interferencias a lo que est谩n intentando recibir.

Si bien esto fue solo una breve descripci贸n general del lado de la transmisi贸n, profundizaremos en m谩s detalles cuando hablemos de la recepci贸n de RDS.

La de recepci贸n

Para demodular y decodificar RDS, realizaremos los siguientes pasos, muchos de los cuales son pasos del lado de transmisi贸n a la inversa (no es necesario memorizar esta lista, recorreremos cada paso individualmente a continuaci贸n):

  1. Reciba una se帽al de radio FM centrada en la frecuencia de la estaci贸n (o lea en una grabaci贸n de IQ), generalmente a una frecuencia de muestreo de 250 kHz

  2. Demodular la FM usando lo que se llama 鈥渄emodulaci贸n en cuadratura鈥

  3. Cambio de frecuencia de 57 kHz para que la se帽al RDS est茅 centrada en 0 Hz

  4. Filtro de paso bajo, para filtrar todo excepto RDS (tambi茅n act煤a como filtro combinado)

  5. Diezmar por 10 para que podamos trabajar con una frecuencia de muestreo m谩s baja, ya que de todos modos filtramos las frecuencias m谩s altas.

  6. Remuestrear a 19 kHz lo que nos dar谩 un n煤mero entero de muestras por s铆mbolo

  7. Sincronizaci贸n de tiempo a nivel de s铆mbolo, usando Mueller y Muller en este ejemplo

  8. Sincronizaci贸n fina de frecuencia mediante un bucle Costas

  9. Demodular el BPSK a 1 y 0

  10. Decodificaci贸n diferencial, para deshacer la codificaci贸n diferencial que se aplic贸

  11. Decodificaci贸n de los 1 y 0 en grupos de bytes

  12. An谩lisis de los grupos de bytes en nuestro resultado final.

Si bien esto puede parecer muchos pasos, RDS es en realidad uno de los protocolos de comunicaciones digitales inal谩mbricas m谩s simples que existen. Un protocolo inal谩mbrico moderno como WiFi o 5G requiere un libro de texto completo para cubrir solo la informaci贸n de la capa PHY/MAC de alto nivel.

Ahora profundizaremos en el c贸digo Python utilizado para recibir RDS. Este c贸digo ha sido probado para funcionar usando una Grabaci贸n de radio FM que puedes encontrar aqu铆., aunque deber铆a poder transmitir su propia se帽al siempre que se reciba con una SNR lo suficientemente alta, simplemente sintonice la frecuencia central de la estaci贸n y muestree a una velocidad de 250 kHz. Para maximizar la potencia de la se帽al recibida (por ejemplo, si est谩 en interiores), es 煤til utilizar una antena dipolo de media onda de la longitud correcta (~1,5 metros), no las antenas de 2,4 GHz que vienen con Pluto. Dicho esto, FM es una se帽al muy fuerte y, si est谩s cerca de una ventana o afuera, las antenas de 2,4 GHz probablemente ser谩n suficientes para captar las estaciones de radio m谩s potentes.

En esta secci贸n presentaremos peque帽as porciones del c贸digo individualmente, con discusi贸n, pero el mismo c贸digo se proporciona al final de este cap铆tulo en un bloque grande. Cada secci贸n presentar谩 un bloque de c贸digo y luego explicar谩 lo que est谩 haciendo.

Adquirir una se帽al

import numpy as np
from scipy.signal import resample_poly, firwin, bilinear, lfilter
import matplotlib.pyplot as plt

# Read in signal
x = np.fromfile('/home/marc/Downloads/fm_rds_250k_1Msamples.iq', dtype=np.complex64)
sample_rate = 250e3
center_freq = 99.5e6

Leemos en nuestra grabaci贸n de prueba, que fue muestreada a 250 kHz y centrada en una estaci贸n de FM recibida con una SNR alta. Aseg煤rese de actualizar la ruta del archivo para reflejar su sistema y d贸nde guard贸 la grabaci贸n. Si ya tiene un SDR configurado y funcionando desde Python, no dude en recibir una se帽al en vivo, aunque es 煤til haber probado primero todo el c贸digo con un known-to-work IQ recording. A lo largo de este c贸digo usaremos x para almacenar la se帽al actual que se est谩 manipulando.

Demodulaci贸n FM

# Quadrature Demod
x = 0.5 * np.angle(x[0:-1] * np.conj(x[1:])) # see https://wiki.gnuradio.org/index.php/Quadrature_Demod

Como se analiz贸 al principio de este cap铆tulo, varias se帽ales individuales se combinan en frecuencia y se modulan en FM para crear lo que realmente se transmite a trav茅s del aire. Entonces el primer paso es deshacer esa modulaci贸n FM. Otra forma de pensarlo es que la informaci贸n se almacena en la variaci贸n de frecuencia de la se帽al que recibimos, y queremos demodularla para que la informaci贸n ahora est茅 en amplitud, no en frecuencia. Tenga en cuenta que la salida de esta demodulaci贸n es una se帽al real, aunque hayamos introducido una se帽al compleja.

Lo que hace esta 煤nica l铆nea de Python es primero calcular el producto de nuestra se帽al con una versi贸n retardada y conjugada de nuestra se帽al. A continuaci贸n, encuentra la fase de cada muestra en ese resultado, que es el momento en el que pasa de complejo a real. Para demostrarnos que esto nos da la informaci贸n contenida en las variaciones de frecuencia, consideremos un tono en la frecuencia f con alguna fase arbitraria \phi, que podemos representar como e^ {j2 \pi (ft + \phi)}. Cuando se trata de tiempo discreto, que utiliza un n煤mero entero n en lugar de t, esto se convierte en e^{j2 \pi (f n + \phi)}. La versi贸n conjugada y retrasada es e^{-j2 \pi (f (n-1) + \phi)}. Multiplicar estos dos lleva a e^{j2 \pi f}, lo cual es genial porque \phi desapareci贸, y cuando calculamos la fase de esa expresi贸n nos queda solo f.

Un efecto secundario conveniente de la modulaci贸n FM es que las variaciones de amplitud de la se帽al recibida en realidad no cambian el volumen del audio, a diferencia de la radio AM.

Dezplazamiento en frecuencia

# Freq shift
N = len(x)
f_o = -57e3 # amount we need to shift by
t = np.arange(N)/sample_rate # time vector
x = x * np.exp(2j*np.pi*f_o*t) # down shift

A continuaci贸n bajamos la frecuencia en 57 kHz, usando el truco e^{j2 \pi f_ot} que aprendimos en el cap铆tulo Sincronizaci贸n donde f_o es el cambio de frecuencia en Hz y t es solo un vector de tiempo, el hecho de que comience en 0 no es importante, lo que importa es que use el per铆odo de muestreo correcto (que es inverso a la frecuencia de muestreo). Adem谩s, debido a que se trata de una se帽al real, en realidad no importa si usas -57 o +57 kHz porque las frecuencias negativas coinciden con las positivas, por lo que de cualquier manera cambiaremos nuestro RDS a 0. Hz.

Filtrar para aislar RDS

# Low-Pass Filter
taps = firwin(numtaps=101, cutoff=7.5e3, fs=sample_rate)
x = np.convolve(x, taps, 'valid')

Ahora debemos filtrar todo excepto RDS. Como tenemos RDS centrado en 0 Hz, eso significa que lo que queremos es un filtro paso bajo. Usamos firwin() para dise帽ar un filtro FIR (es decir, encontrar los taps), que solo necesita saber cu谩ntos taps queremos que tenga el filtro y la frecuencia de corte. Tambi茅n se debe proporcionar la frecuencia de muestreo o, de lo contrario, la frecuencia de corte no tiene sentido para firwin. El resultado es un filtro paso bajo sim茅trico, por lo que sabemos que las derivaciones ser谩n n煤meros reales y podemos aplicar el filtro a nuestra se帽al mediante una convoluci贸n. Elegimos 'valid' para deshacernos de los efectos de borde de hacer convoluci贸n, aunque en este caso realmente no importa porque estamos alimentando una se帽al tan larga que algunas muestras extra帽as en cada borde son innecesarias. No voy a desperdiciar nada.

Nota al margen: en alg煤n momento actualizar茅 el filtro anterior para usar un filtro coincidente adecuado (creo que el coseno elevado de ra铆z es lo que usa RDS), por razones conceptuales, pero obtuve las mismas tasas de error usando el enfoque firwin() que un apropiado filtro de acoplamiento de GNU Radio, por lo que claramente no es un requisito estricto.

Decimate por 10

# Decimate by 10, now that we filtered and there wont be aliasing
x = x[::10]
sample_rate = 25e3

Cada vez que filtra hasta una peque帽a fracci贸n de su ancho de banda (por ejemplo, comenzamos con 125 kHz de ancho de banda real y ahorramos solo 7,5 kHz de eso), tiene sentido diezmar. Recuerde el comienzo del cap铆tulo Muestreo IQ donde aprendimos sobre la frecuencia Nyquist y c贸mo almacenar completamente informaci贸n de banda limitada siempre que muestreemos al doble de la frecuencia m谩s alta. Bueno, ahora que usamos nuestro filtro paso bajo, nuestra frecuencia m谩s alta es de aproximadamente 7,5 kHz, por lo que solo necesitamos una frecuencia de muestreo de 15 kHz. S贸lo para estar seguros agregaremos algo de margen y usaremos una nueva frecuencia de muestreo de 25 kHz (esto terminar谩 funcionando bien matem谩ticamente m谩s adelante).

Realizamos la diezma simplemente descartando 9 de cada 10 muestras, ya que anteriormente ten铆amos una frecuencia de muestreo de 250 kHz y ahora queremos que sea de 25 kHz. Esto puede parecer confuso al principio, porque descartar el 90% de las muestras parece como si estuvieras descartando informaci贸n, pero si revisas el cap铆tulo Muestreo IQ ver谩s por qu茅 en realidad no estamos perdiendo nada, porque Se filtr贸 correctamente (que actu贸 como nuestro filtro anti-aliasing) y redujo nuestra frecuencia m谩xima y, por lo tanto, el ancho de banda de la se帽al.

Desde la perspectiva del c贸digo, este es probablemente el paso m谩s simple de todos, pero aseg煤rese de actualizar su variable sample_rate para reflejar la nueva frecuencia de muestreo.

Remuestreo a 19 kHz

# Resample to 19kHz
x = resample_poly(x, 19, 25) # up, down
sample_rate = 19e3

En el cap铆tulo Formador de Pulso afianzamos el concepto de 鈥渕uestras por s铆mbolo鈥 y aprendimos la conveniencia de tener un n煤mero entero de muestras por s铆mbolo (un valor fraccionario es v谩lido, pero no conveniente). Como se mencion贸 anteriormente, RDS utiliza BPSK y transmite 1187,5 s铆mbolos por segundo. Si continuamos usando nuestra se帽al tal como est谩, muestreada a 25 kHz, tendremos 21.052631579 muestras por s铆mbolo (haga una pausa y piense en los c谩lculos si eso no tiene sentido). Entonces, lo que realmente queremos es una frecuencia de muestreo que sea un m煤ltiplo entero de 1187,5 Hz, pero no podemos bajarla demasiado o no podremos 鈥渁lmacenar鈥 el ancho de banda completo de nuestra se帽al. En la subsecci贸n anterior hablamos de que necesitamos una frecuencia de muestreo de 15 kHz o superior, y elegimos 25 kHz s贸lo para darnos algo de margen.

Encontrar la mejor frecuencia de muestreo para remuestrear se reduce a cu谩ntas muestras por s铆mbolo queremos, y podemos trabajar hacia atr谩s. Hipot茅ticamente, consideremos apuntar a 10 muestras por s铆mbolo. La velocidad de s铆mbolo RDS de 1187,5 multiplicada por 10 nos dar铆a una frecuencia de muestreo de 11,875 kHz, que lamentablemente no es lo suficientemente alta para Nyquist. 驴Qu茅 tal 13 muestras por s铆mbolo? 1187,5 multiplicado por 13 nos da 15437,5 Hz, que est谩 por encima de 15 kHz, pero es un n煤mero bastante impar. 驴Qu茅 tal la siguiente potencia de 2, es decir, 16 muestras por s铆mbolo? 隆1187,5 multiplicado por 16 es exactamente 19 kHz! El n煤mero par es menos una coincidencia y m谩s una elecci贸n de dise帽o de protocolo.

Para remuestrear de 25 kHz a 19 kHz, usamos resample_poly() que aumenta la muestra con un valor entero, filtra y luego reduce la muestra con un valor entero. Esto es conveniente porque en lugar de ingresar 25000 y 19000 podemos usar 25 y 19. Si hubi茅ramos usado 13 muestras por s铆mbolo usando una frecuencia de muestreo de 15437,5 Hz, no podr铆amos usar resample_poly() y el proceso de remuestreo ser铆a mucho m谩s complicado.

Una vez m谩s, recuerde siempre actualizar su variable sample_rate cuando realice una operaci贸n que la cambie.

Sincronizaci贸n en Tiempo (S铆mbolo-Nivel)

# Symbol sync, using what we did in sync chapter
samples = x # for the sake of matching the sync chapter
samples_interpolated = resample_poly(samples, 32, 1) # we'll use 32 as the interpolation factor, arbitrarily chosen, seems to work better than 16
sps = 16
mu = 0.01 # initial estimate of phase of sample
out = np.zeros(len(samples) + 10, dtype=np.complex64)
out_rail = np.zeros(len(samples) + 10, dtype=np.complex64) # stores values, each iteration we need the previous 2 values plus current value
i_in = 0 # input samples index
i_out = 2 # output index (let first two outputs be 0)
while i_out < len(samples) and i_in+32 < len(samples):
    out[i_out] = samples_interpolated[i_in*32 + int(mu*32)] # grab what we think is the "best" sample
    out_rail[i_out] = int(np.real(out[i_out]) > 0) + 1j*int(np.imag(out[i_out]) > 0)
    x = (out_rail[i_out] - out_rail[i_out-2]) * np.conj(out[i_out-1])
    y = (out[i_out] - out[i_out-2]) * np.conj(out_rail[i_out-1])
    mm_val = np.real(y - x)
    mu += sps + 0.01*mm_val
    i_in += int(np.floor(mu)) # round down to nearest int since we are using it as an index
    mu = mu - np.floor(mu) # remove the integer part of mu
    i_out += 1 # increment output index
x = out[2:i_out] # remove the first two, and anything after i_out (that was never filled out)

Finalmente estamos listos para nuestra sincronizaci贸n de s铆mbolo/tiempo, aqu铆 usaremos exactamente el mismo c贸digo de sincronizaci贸n de reloj de Mueller y Muller del cap铆tulo Sincronizaci贸n, cons煤ltelo si desea obtener m谩s informaci贸n sobre c贸mo funciona. Establecemos la muestra por s铆mbolo (sps) en 16 como se analiz贸 anteriormente. Mediante experimentaci贸n se descubri贸 que un valor de ganancia de mu de 0,01 funciona bien. La salida ahora deber铆a ser una muestra por s铆mbolo, es decir, nuestra salida son nuestros 鈥渟铆mbolos suaves鈥, con posible compensaci贸n de frecuencia incluida. La siguiente animaci贸n del gr谩fico de constelaci贸n se utiliza para verificar que estamos obteniendo s铆mbolos BPSK (con un desplazamiento de frecuencia que causa rotaci贸n):

Animation of BPSK rotating because fine frequency sync hasn't been performed yet

Si est谩 utilizando su propia se帽al de FM y no obtiene dos grupos distintos de muestras complejas en este punto, significa que la sincronizaci贸n del s铆mbolo anterior no logr贸 sincronizarse o que hay alg煤n problema con uno de los pasos anteriores. No es necesario animar la constelaci贸n, pero si graficarlas, aseg煤rate de evitar graficar todas las muestras, porque simplemente se ver谩 como un c铆rculo. Si traza s贸lo 100 o 200 muestras a la vez, tendr谩 una mejor idea de si est谩n en dos grupos o no, incluso si est谩n girando.

Sincronizaci贸n en Frecuencia Fina

# Fine freq sync
samples = x # for the sake of matching the sync chapter
N = len(samples)
phase = 0
freq = 0
# These next two params is what to adjust, to make the feedback loop faster or slower (which impacts stability)
alpha = 8.0
beta = 0.002
out = np.zeros(N, dtype=np.complex64)
freq_log = []
for i in range(N):
    out[i] = samples[i] * np.exp(-1j*phase) # adjust the input sample by the inverse of the estimated phase offset
    error = np.real(out[i]) * np.imag(out[i]) # This is the error formula for 2nd order Costas Loop (e.g. for BPSK)

    # Advance the loop (recalc phase and freq offset)
    freq += (beta * error)
    freq_log.append(freq * sample_rate / (2*np.pi)) # convert from angular velocity to Hz for logging
    phase += freq + (alpha * error)

    # Optional: Adjust phase so its always between 0 and 2pi, recall that phase wraps around every 2pi
    while phase >= 2*np.pi:
        phase -= 2*np.pi
    while phase < 0:
        phase += 2*np.pi
x = out

Tambi茅n copiaremos el c贸digo Python de sincronizaci贸n fina en frecuencia del cap铆tulo Sincronizaci贸n, que utiliza Costas Loop para eliminar cualquier desplazamiento de frecuencia residual, as铆 como alinear nuestro BPSK con el eje real (I), forzando Q sea lo m谩s cercano posible a cero. Cualquier cosa que quede en Q probablemente se deba al ruido en la se帽al, suponiendo que el bucle de Costas est茅 sintonizado correctamente. Solo por diversi贸n, veamos la misma animaci贸n que arriba excepto despu茅s de que se haya realizado la sincronizaci贸n de frecuencia (隆no m谩s giros!):

Animation of the frequency sync process using a Costas Loop

Adem谩s, podemos observar el error de frecuencia estimado a lo largo del tiempo para ver c贸mo funcionamiento de Costas Loop; observe c贸mo lo registramos en el c贸digo anterior. Parece que hubo alrededor de 13 Hz de compensaci贸n de frecuencia, ya sea debido a que el oscilador/LO del transmisor estaba apagado o al LO del receptor (muy probablemente el receptor). Si est谩 utilizando su propia se帽al de FM, es posible que necesite modificar alpha y beta hasta que la curva se vea similar; deber铆a lograr la sincronizaci贸n con bastante rapidez (por ejemplo, unos cientos de s铆mbolos) y mantenerla con m铆nima oscilaci贸n. El patr贸n que ve a continuaci贸n despu茅s de encontrar su estado estable es fluctuaci贸n de frecuencia, no oscilaci贸n.

The frequency sync process using a Costas Loop showing the estimated frequency offset over time

Demodulaci贸n BPSK

# Demod BPSK
bits = (np.real(x) > 0).astype(int) # 1's and 0's

Demodular el BPSK en este punto es muy f谩cil, recuerde que cada muestra representa un s铆mbolo suave, por lo que todo lo que tenemos que hacer es verificar si cada muestra est谩 por encima o por debajo de 0. El .astype(int) es as铆 podemos trabajar con una serie de enteros en lugar de una serie de booleanos. Quiz谩s te preguntes si por encima o por debajo de cero representa un 1 o un 0. Como ver谩s en el siguiente paso, 隆no importa!

Decodificaci贸n diferencial

# Differential decoding, so that it doesn't matter whether our BPSK was 180 degrees rotated without us realizing it
bits = (bits[1:] - bits[0:-1]) % 2
bits = bits.astype(np.uint8) # for decoder

La se帽al BPSK utiliz贸 codificaci贸n diferencial cuando se cre贸, lo que significa que cada 1 y 0 de los datos originales se transform贸 de manera que un cambio de 1 a 0 o de 0 a 1 se asign贸 a un 1, y ning煤n cambio se asign贸 a un 0. El gran beneficio de usar codificaci贸n diferencial es que no tienes que preocuparte por las rotaciones de 180 grados al recibir el BPSK, porque si consideramos que un 1 es mayor que cero o menor que cero ya no es un impacto, lo que importa es cambiando entre 1 y 0. Este concepto podr铆a ser m谩s f谩cil de entender si observa datos de ejemplo; a continuaci贸n se muestran los primeros 10 s铆mbolos antes y despu茅s de la decodificaci贸n diferencial:

[1 1 1 1 0 1 0 0 1 1] # before differential decoding
[- 0 0 0 1 1 1 0 1 0] # after differential decoding

Decodificaci贸n RDS

隆Finalmente tenemos nuestros fragmentos de informaci贸n y estamos listos para decodificar lo que significan! El enorme bloque de c贸digo que se proporciona a continuaci贸n es lo que usaremos para decodificar los 1 y 0 en grupos de bytes. Esta parte tendr铆a mucho m谩s sentido si primero cre谩ramos la parte del transmisor de RDS, pero por ahora solo sepa que en RDS, los bytes se agrupan en grupos de 12 bytes, donde los primeros 8 representan los datos y los 煤ltimos 4 act煤an como un palabra de sincronizaci贸n (llamadas 鈥減alabras desplazadas鈥). Los 煤ltimos 4 bytes no son necesarios para el siguiente paso (el analizador), por lo que no los incluimos en la salida. Este bloque de c贸digo toma los 1 y 0 creados anteriormente (en forma de una matriz 1D de uint8) y genera una lista de listas de bytes (una lista de 8 bytes donde esos 8 bytes est谩n en una lista). Esto lo hace conveniente para el siguiente paso, que recorrer谩 la lista de 8 bytes, un grupo de 8 a la vez.

La mayor parte del c贸digo de decodificaci贸n real a continuaci贸n gira en torno a la sincronizaci贸n (a nivel de bytes, no de s铆mbolos) y la verificaci贸n de errores. Funciona en bloques de 104 bits, cada bloque se recibe correctamente o con error (usando CRC para verificar), y cada 50 bloques verifica si m谩s de 35 de ellos se recibieron con error, en cuyo caso reinicia todo e intenta sincronizar nuevamente. El CRC se realiza mediante una verificaci贸n de 10 bits, con polinomio x^{10}+x^8+x^7+x^5+x^4+x^3+1; esto ocurre cuando reg se aplica xor con 0x5B9, que es el equivalente binario de ese polinomio. En Python, los operadores bit a bit para [y, o, no, xor] son & | ~ ^ respectivamente, exactamente igual que C++. Un desplazamiento de bit a la izquierda es x << y (igual que multiplicar x por 2**y), y un desplazamiento de bit a la derecha es x >> y (igual que dividir x por 2** y), tambi茅n como en C++.

Tenga en cuenta que no necesita revisar todo este c贸digo, ni nada de 茅l, especialmente si se est谩 concentrando en aprender el lado de la capa f铆sica (PHY) de DSP y SDR, ya que esto no representa la se帽al. Procesando. Este c贸digo es simplemente una implementaci贸n de un decodificador RDS y, esencialmente, nada de 茅l puede reutilizarse para otros protocolos, porque es muy espec铆fico de la forma en que funciona RDS. Si ya est谩 algo agotado con este cap铆tulo, si茅ntase libre de saltarse este enorme bloque de c贸digo que tiene un trabajo bastante simple pero lo hace de una manera compleja.

# Constants
syndrome = [383, 14, 303, 663, 748]
offset_pos = [0, 1, 2, 3, 2]
offset_word = [252, 408, 360, 436, 848]

# see Annex B, page 64 of the standard
def calc_syndrome(x, mlen):
    reg = 0
    plen = 10
    for ii in range(mlen, 0, -1):
        reg = (reg << 1) | ((x >> (ii-1)) & 0x01)
        if (reg & (1 << plen)):
            reg = reg ^ 0x5B9
    for ii in range(plen, 0, -1):
        reg = reg << 1
        if (reg & (1 << plen)):
            reg = reg ^ 0x5B9
    return reg & ((1 << plen) - 1) # select the bottom plen bits of reg

# Initialize all the working vars we'll need during the loop
synced = False
presync = False

wrong_blocks_counter = 0
blocks_counter = 0
group_good_blocks_counter = 0

reg = np.uint32(0) # was unsigned long in C++ (64 bits) but numpy doesn't support bitwise ops of uint64, I don't think it gets that high anyway
lastseen_offset_counter = 0
lastseen_offset = 0

# the synchronization process is described in Annex C, page 66 of the standard */
bytes_out = []
for i in range(len(bits)):
    # in C++ reg doesn't get init so it will be random at first, for ours its 0s
    # It was also an unsigned long but never seemed to get anywhere near the max value
    # bits are either 0 or 1
    reg = np.bitwise_or(np.left_shift(reg, 1), bits[i]) # reg contains the last 26 rds bits. these are both bitwise ops
    if not synced:
        reg_syndrome = calc_syndrome(reg, 26)
        for j in range(5):
            if reg_syndrome == syndrome[j]:
                if not presync:
                    lastseen_offset = j
                    lastseen_offset_counter = i
                    presync = True
                else:
                    if offset_pos[lastseen_offset] >= offset_pos[j]:
                        block_distance = offset_pos[j] + 4 - offset_pos[lastseen_offset]
                    else:
                        block_distance = offset_pos[j] - offset_pos[lastseen_offset]
                    if (block_distance*26) != (i - lastseen_offset_counter):
                        presync = False
                    else:
                        print('Sync State Detected')
                        wrong_blocks_counter = 0
                        blocks_counter = 0
                        block_bit_counter = 0
                        block_number = (j + 1) % 4
                        group_assembly_started = False
                        synced = True
            break # syndrome found, no more cycles

    else: # SYNCED
        # wait until 26 bits enter the buffer */
        if block_bit_counter < 25:
            block_bit_counter += 1
        else:
            good_block = False
            dataword = (reg >> 10) & 0xffff
            block_calculated_crc = calc_syndrome(dataword, 16)
            checkword = reg & 0x3ff
            if block_number == 2: # manage special case of C or C' offset word
                block_received_crc = checkword ^ offset_word[block_number]
                if (block_received_crc == block_calculated_crc):
                    good_block = True
                else:
                    block_received_crc = checkword ^ offset_word[4]
                    if (block_received_crc == block_calculated_crc):
                        good_block = True
                    else:
                        wrong_blocks_counter += 1
                        good_block = False
            else:
                block_received_crc = checkword ^ offset_word[block_number] # bitwise xor
                if block_received_crc == block_calculated_crc:
                    good_block = True
                else:
                    wrong_blocks_counter += 1
                    good_block = False

            # Done checking CRC
            if block_number == 0 and good_block:
                group_assembly_started = True
                group_good_blocks_counter = 1
                bytes = bytearray(8) # 8 bytes filled with 0s
            if group_assembly_started:
                if not good_block:
                    group_assembly_started = False
                else:
                    # raw data bytes, as received from RDS. 8 info bytes, followed by 4 RDS offset chars: ABCD/ABcD/EEEE (in US) which we leave out here
                    # RDS information words
                    # block_number is either 0,1,2,3 so this is how we fill out the 8 bytes
                    bytes[block_number*2] = (dataword >> 8) & 255
                    bytes[block_number*2+1] = dataword & 255
                    group_good_blocks_counter += 1
                    #print('group_good_blocks_counter:', group_good_blocks_counter)
                if group_good_blocks_counter == 5:
                    #print(bytes)
                    bytes_out.append(bytes) # list of len-8 lists of bytes
            block_bit_counter = 0
            block_number = (block_number + 1) % 4
            blocks_counter += 1
            if blocks_counter == 50:
                if wrong_blocks_counter > 35: # This many wrong blocks must mean we lost sync
                    print("Lost Sync (Got ", wrong_blocks_counter, " bad blocks on ", blocks_counter, " total)")
                    synced = False
                    presync = False
                else:
                    print("Still Sync-ed (Got ", wrong_blocks_counter, " bad blocks on ", blocks_counter, " total)")
                blocks_counter = 0
                wrong_blocks_counter = 0

A continuaci贸n se muestra un ejemplo de resultado de este paso de decodificaci贸n. Observe c贸mo en este ejemplo se sincroniz贸 con bastante rapidez pero luego pierde la sincronizaci贸n un par de veces por alg煤n motivo, aunque a煤n puede analizar todos los datos, como veremos. Si est谩 utilizando el archivo de muestra descargable de 1M, solo ver谩 las primeras l铆neas a continuaci贸n. El contenido real de estos bytes simplemente parece n煤meros/caracteres aleatorios dependiendo de c贸mo los muestre, pero en el siguiente paso los analizaremos para convertirlos en informaci贸n legible por humanos.

Sync State Detected
Still Sync-ed (Got  0  bad blocks on  50  total)
Still Sync-ed (Got  0  bad blocks on  50  total)
Still Sync-ed (Got  0  bad blocks on  50  total)
Still Sync-ed (Got  0  bad blocks on  50  total)
Still Sync-ed (Got  1  bad blocks on  50  total)
Still Sync-ed (Got  5  bad blocks on  50  total)
Still Sync-ed (Got  26  bad blocks on  50  total)
Lost Sync (Got  50  bad blocks on  50  total)
Sync State Detected
Still Sync-ed (Got  3  bad blocks on  50  total)
Still Sync-ed (Got  0  bad blocks on  50  total)
Still Sync-ed (Got  0  bad blocks on  50  total)
Still Sync-ed (Got  0  bad blocks on  50  total)
Still Sync-ed (Got  0  bad blocks on  50  total)
Still Sync-ed (Got  0  bad blocks on  50  total)
Still Sync-ed (Got  0  bad blocks on  50  total)
Still Sync-ed (Got  0  bad blocks on  50  total)
Still Sync-ed (Got  0  bad blocks on  50  total)
Still Sync-ed (Got  0  bad blocks on  50  total)
Still Sync-ed (Got  0  bad blocks on  50  total)
Still Sync-ed (Got  0  bad blocks on  50  total)
Still Sync-ed (Got  0  bad blocks on  50  total)
Still Sync-ed (Got  0  bad blocks on  50  total)
Still Sync-ed (Got  0  bad blocks on  50  total)
Still Sync-ed (Got  0  bad blocks on  50  total)
Still Sync-ed (Got  0  bad blocks on  50  total)
Still Sync-ed (Got  0  bad blocks on  50  total)
Still Sync-ed (Got  0  bad blocks on  50  total)
Still Sync-ed (Got  0  bad blocks on  50  total)
Still Sync-ed (Got  0  bad blocks on  50  total)
Still Sync-ed (Got  0  bad blocks on  50  total)
Still Sync-ed (Got  2  bad blocks on  50  total)
Still Sync-ed (Got  1  bad blocks on  50  total)
Still Sync-ed (Got  20  bad blocks on  50  total)
Lost Sync (Got  47  bad blocks on  50  total)
Sync State Detected
Still Sync-ed (Got  32  bad blocks on  50  total)

An谩lisis RDS

Ahora que tenemos bytes, en grupos de 8, podemos extraer los datos finales, es decir, el resultado final que sea comprensible para los humanos. Esto se conoce como an谩lisis de bytes y, al igual que el decodificador de la secci贸n anterior, es simplemente una implementaci贸n del protocolo RDS y en realidad no es tan importante entenderlo. Por suerte no es un mont贸n de c贸digo, si no incluyes las dos tablas definidas al principio, que son simplemente las tablas de b煤squeda para el tipo de canal FM y el 谩rea de cobertura.

Para aquellos que quieran aprender c贸mo funciona este c贸digo, les proporcionar茅 informaci贸n adicional. El protocolo utiliza este concepto de indicador A/B, lo que significa que algunos mensajes est谩n marcados como A y otros como B, y el an谩lisis cambia seg煤n cu谩l (si es A o B se almacena en el tercer bit del segundo byte). Tambi茅n utiliza diferentes tipos de 鈥済rupo鈥 que son an谩logos al tipo de mensaje, y en este c贸digo solo estamos analizando el tipo de mensaje 2, que es el tipo de mensaje que tiene el texto de radio, que es la parte interesante, es el texto que se desplaza por la pantalla de su autom贸vil. A煤n podremos analizar el tipo de canal y la regi贸n, ya que est谩n almacenados en cada mensaje. Por 煤ltimo, tenga en cuenta que radiotext es una cadena que se inicializa en todos los espacios, se completa lentamente a medida que se analizan los bytes y luego se restablece en todos los espacios si se recibe un conjunto espec铆fico de bytes. Si tiene curiosidad sobre qu茅 otros tipos de mensajes existen, la lista es: [鈥淏ASIC鈥, 鈥淧IN/SL鈥, 鈥淩T鈥, 鈥淎ID鈥, 鈥淐T鈥, 鈥淭DC鈥, 鈥淚H鈥, 鈥淩P鈥, 鈥 TMC鈥, 鈥淓WS鈥, 鈥淓ON鈥漖. El mensaje 鈥淩T鈥 es radiotexto que es el 煤nico que decodificamos. El bloque RDS GNU Radio tambi茅n decodifica 鈥淏ASIC鈥, pero para las estaciones que utilic茅 para probar no conten铆a mucha informaci贸n interesante y habr铆a agregado muchas l铆neas al c贸digo siguiente.

# Annex F of RBDS Standard Table F.1 (North America) and Table F.2 (Europe)
#              Europe                   North America
pty_table = [["Undefined",             "Undefined"],
             ["News",                  "News"],
             ["Current Affairs",       "Information"],
             ["Information",           "Sports"],
             ["Sport",                 "Talk"],
             ["Education",             "Rock"],
             ["Drama",                 "Classic Rock"],
             ["Culture",               "Adult Hits"],
             ["Science",               "Soft Rock"],
             ["Varied",                "Top 40"],
             ["Pop Music",             "Country"],
             ["Rock Music",            "Oldies"],
             ["Easy Listening",        "Soft"],
             ["Light Classical",       "Nostalgia"],
             ["Serious Classical",     "Jazz"],
             ["Other Music",           "Classical"],
             ["Weather",               "Rhythm & Blues"],
             ["Finance",               "Soft Rhythm & Blues"],
             ["Children鈥檚 Programmes", "Language"],
             ["Social Affairs",        "Religious Music"],
             ["Religion",              "Religious Talk"],
             ["Phone-In",              "Personality"],
             ["Travel",                "Public"],
             ["Leisure",               "College"],
             ["Jazz Music",            "Spanish Talk"],
             ["Country Music",         "Spanish Music"],
             ["National Music",        "Hip Hop"],
             ["Oldies Music",          "Unassigned"],
             ["Folk Music",            "Unassigned"],
             ["Documentary",           "Weather"],
             ["Alarm Test",            "Emergency Test"],
             ["Alarm",                 "Emergency"]]
pty_locale = 1 # set to 0 for Europe which will use first column instead

# page 72, Annex D, table D.2 in the standard
coverage_area_codes = ["Local",
                       "International",
                       "National",
                       "Supra-regional",
                       "Regional 1",
                       "Regional 2",
                       "Regional 3",
                       "Regional 4",
                       "Regional 5",
                       "Regional 6",
                       "Regional 7",
                       "Regional 8",
                       "Regional 9",
                       "Regional 10",
                       "Regional 11",
                       "Regional 12"]

radiotext_AB_flag = 0
radiotext = [' ']*65
first_time = True
for bytes in bytes_out:
    group_0 = bytes[1] | (bytes[0] << 8)
    group_1 = bytes[3] | (bytes[2] << 8)
    group_2 = bytes[5] | (bytes[4] << 8)
    group_3 = bytes[7] | (bytes[6] << 8)

    group_type = (group_1 >> 12) & 0xf # here is what each one means, e.g. RT is radiotext which is the only one we decode here: ["BASIC", "PIN/SL", "RT", "AID", "CT", "TDC", "IH", "RP", "TMC", "EWS", "___", "___", "___", "___", "EON", "___"]
    AB = (group_1 >> 11 ) & 0x1 # b if 1, a if 0

    #print("group_type:", group_type) # this is essentially message type, i only see type 0 and 2 in my recording
    #print("AB:", AB)

    program_identification = group_0     # "PI"

    program_type = (group_1 >> 5) & 0x1f # "PTY"
    pty = pty_table[program_type][pty_locale]

    pi_area_coverage = (program_identification >> 8) & 0xf
    coverage_area = coverage_area_codes[pi_area_coverage]

    pi_program_reference_number = program_identification & 0xff # just an int

    if first_time:
        print("PTY:", pty)
        print("program:", pi_program_reference_number)
        print("coverage_area:", coverage_area)
        first_time = False

    if group_type == 2:
        # when the A/B flag is toggled, flush your current radiotext
        if radiotext_AB_flag != ((group_1 >> 4) & 0x01):
            radiotext = [' ']*65
        radiotext_AB_flag = (group_1 >> 4) & 0x01
        text_segment_address_code = group_1 & 0x0f
        if AB:
            radiotext[text_segment_address_code * 2    ] = chr((group_3 >> 8) & 0xff)
            radiotext[text_segment_address_code * 2 + 1] = chr(group_3        & 0xff)
        else:
            radiotext[text_segment_address_code *4     ] = chr((group_2 >> 8) & 0xff)
            radiotext[text_segment_address_code * 4 + 1] = chr(group_2        & 0xff)
            radiotext[text_segment_address_code * 4 + 2] = chr((group_3 >> 8) & 0xff)
            radiotext[text_segment_address_code * 4 + 3] = chr(group_3        & 0xff)
        print(''.join(radiotext))
    else:
        pass
        #print("unsupported group_type:", group_type)

A continuaci贸n se muestra el resultado del paso de an谩lisis para una estaci贸n de FM de ejemplo. Observe c贸mo tiene que construir la cadena de radiotexto sobre m煤ltiples mensajes, y luego peri贸dicamente borra la cadena y comienza de nuevo. Si est谩 utilizando el archivo descargado de 1 millon de muestras, solo ver谩 las primeras l铆neas a continuaci贸n.

PTY: Top 40
program: 29
coverage_area: Regional 4
            ing.
            ing. Upb
            ing. Upbeat.
            ing. Upbeat. Rea

WAY-
WAY-FM U
WAY-FM Uplif
WAY-FM Uplifting
WAY-FM Uplifting. Up
WAY-FM Uplifting. Upbeat
WAY-FM Uplifting. Upbeat. Re

WayF
WayFM Up
WayFM Uplift
WayFM Uplifting.
WayFM Uplifting. Upb
WayFM Uplifting. Upbeat.
WayFM Uplifting. Upbeat. Rea

Resumen y c贸digo final

隆Lo hiciste! A continuaci贸n se muestra todo el c贸digo anterior, concatenado, deber铆a funcionar con la prueba de grabaci贸n de radio FM que puedes encontrar aqu铆, aunque deber铆a poder transmitir su propia se帽al siempre que se reciba con una SNR lo suficientemente alta, simplemente sintonice la frecuencia central de la estaci贸n y muestree a una velocidad de 250 kHz. Si descubre que tuvo que hacer ajustes para que funcione con su propia grabaci贸n o SDR en vivo, h谩game saber lo que tuvo que hacer; puede enviarlo como Pull Request (PR) de GitHub en la p谩gina de GitHub del libro de texto. Tambi茅n puede encontrar una versi贸n de este c贸digo con docenas de graficas/impresi贸n de depuraci贸n incluidos, que utilic茅 originalmente para hacer este cap铆tulo. aqui.

Final Code
import numpy as np
from scipy.signal import resample_poly, firwin, bilinear, lfilter
import matplotlib.pyplot as plt

# Read in signal
x = np.fromfile('/home/marc/Downloads/fm_rds_250k_from_sdrplay.iq', dtype=np.complex64)
sample_rate = 250e3
center_freq = 99.5e6

# Quadrature Demod
x = 0.5 * np.angle(x[0:-1] * np.conj(x[1:])) # see https://wiki.gnuradio.org/index.php/Quadrature_Demod

# Freq shift
N = len(x)
f_o = -57e3 # amount we need to shift by
t = np.arange(N)/sample_rate # time vector
x = x * np.exp(2j*np.pi*f_o*t) # down shift

# Low-Pass Filter
taps = firwin(numtaps=101, cutoff=7.5e3, fs=sample_rate)
x = np.convolve(x, taps, 'valid')

# Decimate by 10, now that we filtered and there wont be aliasing
x = x[::10]
sample_rate = 25e3

# Resample to 19kHz
x = resample_poly(x, 19, 25) # up, down
sample_rate = 19e3

# Symbol sync, using what we did in sync chapter
samples = x # for the sake of matching the sync chapter
samples_interpolated = resample_poly(samples, 32, 1) # we'll use 32 as the interpolation factor, arbitrarily chosen
sps = 16
mu = 0.01 # initial estimate of phase of sample
out = np.zeros(len(samples) + 10, dtype=np.complex64)
out_rail = np.zeros(len(samples) + 10, dtype=np.complex64) # stores values, each iteration we need the previous 2 values plus current value
i_in = 0 # input samples index
i_out = 2 # output index (let first two outputs be 0)
while i_out < len(samples) and i_in+32 < len(samples):
    out[i_out] = samples_interpolated[i_in*32 + int(mu*32)] # grab what we think is the "best" sample
    out_rail[i_out] = int(np.real(out[i_out]) > 0) + 1j*int(np.imag(out[i_out]) > 0)
    x = (out_rail[i_out] - out_rail[i_out-2]) * np.conj(out[i_out-1])
    y = (out[i_out] - out[i_out-2]) * np.conj(out_rail[i_out-1])
    mm_val = np.real(y - x)
    mu += sps + 0.01*mm_val
    i_in += int(np.floor(mu)) # round down to nearest int since we are using it as an index
    mu = mu - np.floor(mu) # remove the integer part of mu
    i_out += 1 # increment output index
x = out[2:i_out] # remove the first two, and anything after i_out (that was never filled out)

#new sample_rate should be 1187.5
sample_rate /= 16

# Fine freq sync
samples = x # for the sake of matching the sync chapter
N = len(samples)
phase = 0
freq = 0
# These next two params is what to adjust, to make the feedback loop faster or slower (which impacts stability)
alpha = 8.0
beta = 0.002
out = np.zeros(N, dtype=np.complex64)
freq_log = []
for i in range(N):
    out[i] = samples[i] * np.exp(-1j*phase) # adjust the input sample by the inverse of the estimated phase offset
    error = np.real(out[i]) * np.imag(out[i]) # This is the error formula for 2nd order Costas Loop (e.g. for BPSK)

    # Advance the loop (recalc phase and freq offset)
    freq += (beta * error)
    freq_log.append(freq * sample_rate / (2*np.pi)) # convert from angular velocity to Hz for logging
    phase += freq + (alpha * error)

    # Optional: Adjust phase so its always between 0 and 2pi, recall that phase wraps around every 2pi
    while phase >= 2*np.pi:
        phase -= 2*np.pi
    while phase < 0:
        phase += 2*np.pi
x = out

# Demod BPSK
bits = (np.real(x) > 0).astype(int) # 1's and 0's

# Differential decoding, so that it doesn't matter whether our BPSK was 180 degrees rotated without us realizing it
bits = (bits[1:] - bits[0:-1]) % 2
bits = bits.astype(np.uint8) # for decoder

###########
# DECODER #
###########

# Constants
syndrome = [383, 14, 303, 663, 748]
offset_pos = [0, 1, 2, 3, 2]
offset_word = [252, 408, 360, 436, 848]

# see Annex B, page 64 of the standard
def calc_syndrome(x, mlen):
    reg = 0
    plen = 10
    for ii in range(mlen, 0, -1):
        reg = (reg << 1) | ((x >> (ii-1)) & 0x01)
        if (reg & (1 << plen)):
            reg = reg ^ 0x5B9
    for ii in range(plen, 0, -1):
        reg = reg << 1
        if (reg & (1 << plen)):
            reg = reg ^ 0x5B9
    return reg & ((1 << plen) - 1) # select the bottom plen bits of reg

# Initialize all the working vars we'll need during the loop
synced = False
presync = False

wrong_blocks_counter = 0
blocks_counter = 0
group_good_blocks_counter = 0

reg = np.uint32(0) # was unsigned long in C++ (64 bits) but numpy doesn't support bitwise ops of uint64, I don't think it gets that high anyway
lastseen_offset_counter = 0
lastseen_offset = 0

# the synchronization process is described in Annex C, page 66 of the standard */
bytes_out = []
for i in range(len(bits)):
    # in C++ reg doesn't get init so it will be random at first, for ours its 0s
    # It was also an unsigned long but never seemed to get anywhere near the max value
    # bits are either 0 or 1
    reg = np.bitwise_or(np.left_shift(reg, 1), bits[i]) # reg contains the last 26 rds bits. these are both bitwise ops
    if not synced:
        reg_syndrome = calc_syndrome(reg, 26)
        for j in range(5):
            if reg_syndrome == syndrome[j]:
                if not presync:
                    lastseen_offset = j
                    lastseen_offset_counter = i
                    presync = True
                else:
                    if offset_pos[lastseen_offset] >= offset_pos[j]:
                        block_distance = offset_pos[j] + 4 - offset_pos[lastseen_offset]
                    else:
                        block_distance = offset_pos[j] - offset_pos[lastseen_offset]
                    if (block_distance*26) != (i - lastseen_offset_counter):
                        presync = False
                    else:
                        print('Sync State Detected')
                        wrong_blocks_counter = 0
                        blocks_counter = 0
                        block_bit_counter = 0
                        block_number = (j + 1) % 4
                        group_assembly_started = False
                        synced = True
            break # syndrome found, no more cycles

    else: # SYNCED
        # wait until 26 bits enter the buffer */
        if block_bit_counter < 25:
            block_bit_counter += 1
        else:
            good_block = False
            dataword = (reg >> 10) & 0xffff
            block_calculated_crc = calc_syndrome(dataword, 16)
            checkword = reg & 0x3ff
            if block_number == 2: # manage special case of C or C' offset word
                block_received_crc = checkword ^ offset_word[block_number]
                if (block_received_crc == block_calculated_crc):
                    good_block = True
                else:
                    block_received_crc = checkword ^ offset_word[4]
                    if (block_received_crc == block_calculated_crc):
                        good_block = True
                    else:
                        wrong_blocks_counter += 1
                        good_block = False
            else:
                block_received_crc = checkword ^ offset_word[block_number] # bitwise xor
                if block_received_crc == block_calculated_crc:
                    good_block = True
                else:
                    wrong_blocks_counter += 1
                    good_block = False

            # Done checking CRC
            if block_number == 0 and good_block:
                group_assembly_started = True
                group_good_blocks_counter = 1
                bytes = bytearray(8) # 8 bytes filled with 0s
            if group_assembly_started:
                if not good_block:
                    group_assembly_started = False
                else:
                    # raw data bytes, as received from RDS. 8 info bytes, followed by 4 RDS offset chars: ABCD/ABcD/EEEE (in US) which we leave out here
                    # RDS information words
                    # block_number is either 0,1,2,3 so this is how we fill out the 8 bytes
                    bytes[block_number*2] = (dataword >> 8) & 255
                    bytes[block_number*2+1] = dataword & 255
                    group_good_blocks_counter += 1
                    #print('group_good_blocks_counter:', group_good_blocks_counter)
                if group_good_blocks_counter == 5:
                    #print(bytes)
                    bytes_out.append(bytes) # list of len-8 lists of bytes
            block_bit_counter = 0
            block_number = (block_number + 1) % 4
            blocks_counter += 1
            if blocks_counter == 50:
                if wrong_blocks_counter > 35: # This many wrong blocks must mean we lost sync
                    print("Lost Sync (Got ", wrong_blocks_counter, " bad blocks on ", blocks_counter, " total)")
                    synced = False
                    presync = False
                else:
                    print("Still Sync-ed (Got ", wrong_blocks_counter, " bad blocks on ", blocks_counter, " total)")
                blocks_counter = 0
                wrong_blocks_counter = 0

###########
# PARSER  #
###########

# Annex F of RBDS Standard Table F.1 (North America) and Table F.2 (Europe)
#              Europe                   North America
pty_table = [["Undefined",             "Undefined"],
             ["News",                  "News"],
             ["Current Affairs",       "Information"],
             ["Information",           "Sports"],
             ["Sport",                 "Talk"],
             ["Education",             "Rock"],
             ["Drama",                 "Classic Rock"],
             ["Culture",               "Adult Hits"],
             ["Science",               "Soft Rock"],
             ["Varied",                "Top 40"],
             ["Pop Music",             "Country"],
             ["Rock Music",            "Oldies"],
             ["Easy Listening",        "Soft"],
             ["Light Classical",       "Nostalgia"],
             ["Serious Classical",     "Jazz"],
             ["Other Music",           "Classical"],
             ["Weather",               "Rhythm & Blues"],
             ["Finance",               "Soft Rhythm & Blues"],
             ["Children鈥檚 Programmes", "Language"],
             ["Social Affairs",        "Religious Music"],
             ["Religion",              "Religious Talk"],
             ["Phone-In",              "Personality"],
             ["Travel",                "Public"],
             ["Leisure",               "College"],
             ["Jazz Music",            "Spanish Talk"],
             ["Country Music",         "Spanish Music"],
             ["National Music",        "Hip Hop"],
             ["Oldies Music",          "Unassigned"],
             ["Folk Music",            "Unassigned"],
             ["Documentary",           "Weather"],
             ["Alarm Test",            "Emergency Test"],
             ["Alarm",                 "Emergency"]]
pty_locale = 1 # set to 0 for Europe which will use first column instead

# page 72, Annex D, table D.2 in the standard
coverage_area_codes = ["Local",
                       "International",
                       "National",
                       "Supra-regional",
                       "Regional 1",
                       "Regional 2",
                       "Regional 3",
                       "Regional 4",
                       "Regional 5",
                       "Regional 6",
                       "Regional 7",
                       "Regional 8",
                       "Regional 9",
                       "Regional 10",
                       "Regional 11",
                       "Regional 12"]

radiotext_AB_flag = 0
radiotext = [' ']*65
first_time = True
for bytes in bytes_out:
    group_0 = bytes[1] | (bytes[0] << 8)
    group_1 = bytes[3] | (bytes[2] << 8)
    group_2 = bytes[5] | (bytes[4] << 8)
    group_3 = bytes[7] | (bytes[6] << 8)

    group_type = (group_1 >> 12) & 0xf # here is what each one means, e.g. RT is radiotext which is the only one we decode here: ["BASIC", "PIN/SL", "RT", "AID", "CT", "TDC", "IH", "RP", "TMC", "EWS", "___", "___", "___", "___", "EON", "___"]
    AB = (group_1 >> 11 ) & 0x1 # b if 1, a if 0

    #print("group_type:", group_type) # this is essentially message type, i only see type 0 and 2 in my recording
    #print("AB:", AB)

    program_identification = group_0     # "PI"

    program_type = (group_1 >> 5) & 0x1f # "PTY"
    pty = pty_table[program_type][pty_locale]

    pi_area_coverage = (program_identification >> 8) & 0xf
    coverage_area = coverage_area_codes[pi_area_coverage]

    pi_program_reference_number = program_identification & 0xff # just an int

    if first_time:
        print("PTY:", pty)
        print("program:", pi_program_reference_number)
        print("coverage_area:", coverage_area)
        first_time = False

    if group_type == 2:
        # when the A/B flag is toggled, flush your current radiotext
        if radiotext_AB_flag != ((group_1 >> 4) & 0x01):
            radiotext = [' ']*65
        radiotext_AB_flag = (group_1 >> 4) & 0x01
        text_segment_address_code = group_1 & 0x0f
        if AB:
            radiotext[text_segment_address_code * 2    ] = chr((group_3 >> 8) & 0xff)
            radiotext[text_segment_address_code * 2 + 1] = chr(group_3        & 0xff)
        else:
            radiotext[text_segment_address_code *4     ] = chr((group_2 >> 8) & 0xff)
            radiotext[text_segment_address_code * 4 + 1] = chr(group_2        & 0xff)
            radiotext[text_segment_address_code * 4 + 2] = chr((group_3 >> 8) & 0xff)
            radiotext[text_segment_address_code * 4 + 3] = chr(group_3        & 0xff)
        print(''.join(radiotext))
    else:
        pass
        #print("unsupported group_type:", group_type)

Una vez m谩s, el ejemplo de grabaci贸n de FM que funciona con este c贸digo lo puede encontra aqui.

Para aquellos interesados en demodular la se帽al de audio real, simplemente agregue las siguientes l铆neas justo despu茅s de la secci贸n 鈥淎dquirir una se帽al鈥 (agradecimiento especial a Joel Cordeiro por el c贸digo):

# Add the following code right after the "Acquiring a Signal" section

from scipy.io import wavfile

# Demodulation
x = np.diff(np.unwrap(np.angle(x)))

# De-emphasis filter, H(s) = 1/(RC*s + 1), implemented as IIR via bilinear transform
bz, az = bilinear(1, [75e-6, 1], fs=sample_rate)
x = lfilter(bz, az, x)

# decimate by 6 to get mono audio
x = x[::6]
sample_rate_audio = sample_rate/6

# normalize volume so its between -1 and +1
x /= np.max(np.abs(x))

# some machines want int16s
x *= 32767
x = x.astype(np.int16)

# Save to wav file, you can open this in Audacity for example
wavfile.write('fm.wav', int(sample_rate_audio), x)

La parte m谩s complicada es el filtro de de-emphasis, sobre el cual puedes aprender aqu铆, aunque en realidad es un paso opcional si est谩s de acuerdo con el audio que tiene un equilibrio deficiente de graves y agudos. Para aquellos curiosos, aqu铆 est谩 cu谩l es la respuesta de frecuencia del IIR parece que el filtro de de-emphasis no filtra completamente ninguna frecuencia, es m谩s bien un filtro 鈥渇ormador鈥.

../_images/fm_demph_filter_freq_response.svg

Agradecimientos GNU Radio

La mayor铆a de los pasos anteriores utilizados para recibir RDS se adaptaron de la implementaci贸n de RDS de GNU Radio, que se encuentra en el m贸dulo fuera del 谩rbol de GNU Radio llamado gr-rds, creado originalmente por Dimitrios Symeonidis y mantenido por Bastian Bloessl, y me gustar铆a reconocer el trabajo de estos autores. Para crear este cap铆tulo, comenc茅 a usar gr-rds en GNU Radio, con una grabaci贸n de FM funcional, y poco a poco convert铆 cada uno de los bloques (incluidos muchos bloques integrados) a Python. Tom贸 bastante tiempo, hay algunos matices en los bloques integrados que son f谩ciles de pasar por alto, y pasar del procesamiento del formato de la se帽al (es decir, usar una funci贸n de trabajo que procesa unos pocos miles de muestras a la vez en un manera estable) a un bloque de Python no siempre es sencillo. GNU Radio es una herramienta incre铆ble para este tipo de creaci贸n de prototipos y nunca habr铆a podido crear todo este c贸digo Python funcional sin ella.

Lectura Futuras.

  1. https://en.wikipedia.org/wiki/Radio_Data_System

  2. https://www.sigidwiki.com/wiki/Radio_Data_System_(RDS)

  3. https://github.com/bastibl/gr-rds

  4. https://www.gnuradio.org/