La señal horaria TDF en 162 kHz (II)

En la anterior entrada hablaba de las señales horarias que se encuentran en la radio y explicaba con mayor detalle la modulación de la emisora de onda larga TDF en 162 kHz que transmite desde la localidad francesa de Allouis. En este artículo muestro un script de octave que decodifica los datos. Se repasan conceptos de comunicaciones digitales y teoría de señal como son la relación señal/ruido, la sincronización de símbolo o la detección de errores. Al final enumero posibles mejoras futuras y aplicaciones del script.

En España, con una antena dipolo típica para HF y un transceptor para estas bandas, es posible recibir la señal suficientemente bien, aunque débil a causa de la gran desadaptación de impedancias y de que los equipos no están optimizados para frecuencias tan bajas y suelen sufrir saturación de señales más fuertes de onda media. Esto se puede mejorar con un filtro paso banda como éste, diseñado y construido con éxito por Ramiro EA4NZ. La sección paso bajo formada por L4-C8 y el atenuador R1-R2-R3 fueron ideados por Carlos EB4FBZ:

Esto mejora la recepción de la radiodifusión en onda larga, radiofaros y señales de aficionados en 137 o 475 kHz. Si construimos nuestro propio acoplador de antena se puede incorporar ese circuito como opción seleccionable.

Volviendo al procesamiento de la señal, el punto de partida es recibirla en USB como tono de 1 kHz y digitalizarla con una frecuencia de muestreo de 8 kHz. Esto simplifica mucho el proceso puesto que desfasar 90º es tan simple como retrasar dos muestras:

Es conveniente que el receptor aplique el filtro más estrecho disponible. Por ejemplo podemos ajustar el equipo en modo CW USB con un filtro de 300 Hz de anchura. El modo CW suele desplazar la frecuencia de la pantalla a un tono audible de 700 Hz por lo que deberíamos sintonizar el equipo en 161,7 kHz para que quede en 1 kHz, y desplazar análogamente el centro del filtro de CW. En cualquier caso es necesario un filtrado posterior en el ordenador, como hace el script propuesto.

Como decíamos, al multiplicar la señal por una copia de sí misma desplazada obtenemos dos componentes de espectrales: una de frecuencia doble (2 kHz) y otra en banda base, de amplitud proporcional al coseno del desfase. Si el desfase es de dos muestras, esto son 90º en un tono de 1 kHz muestreado a 8 kHz, y por tanto la componente DC será nula. Ahora bien, cuando el emisor varía la fase para enviar información, la frecuencia ya no es exactamente 1 kHz y la componente de baja frecuencia representa la demodulación de la señal. Como indica la entrada en la wiki, hay modulación todo el tiempo excepto durante el segundo 59 de cada minuto. Esto se aprecia claramente en la siguiente imagen:

El eje horizontal indica el número de muestras. La frecuencia de muestreo se ha reducido de 8000 a 500 Hz. El círculo marca la transmisión del primer pulso en el segundo 0. La señal primero baja, luego sube marcando un pulso más ancho y luego baja de nuevo. Esto se corresponde con la modulación transmitida, si bien la gráfica superior representa frecuencia (derivada de la fase respecto del tiempo) y la inferior representa fase. Hay, además, una inversión. La rampa creciente de fase (tiempo de 950 a 975 ms) se corresponde con una frecuencia más alta mientras que la señal demodulada más arriba se hace negativa en ese tiempo. Eso es una simple inversión matemática consecuencia de que la demodulación resulta proporcional a cos(π/2+θ)=sen(-θ)=-sen(θ)≈-θ y por tanto negativa; se deja así por comodidad. Recordemos que el incremento de fase θ se corresponde con la rampa lineal de 40 radianes por segundo observada entre tiempos que distan 2 muestras con una frecuencia de muestreo de 8 kHz, por lo tanto: 40*2/8000=0,01 rad.

La señal demodulada tiene ruido de frecuencias altas como se hace evidente en la gráfica y también se observa en esta otra del espectro:

Intuitivamente podemos pensar en aplicar un filtro pero no es simple la elección. Un filtro podría ensanchar la forma temporal del pulso de modo que se extendiese hasta el instante en el que pretendemos observar el siguiente pulso: esto se llama interferencia entre símbolos (ISI). Sería estúpido distorsionar la señal y perder capacidad de demodularla correctamente porque hacemos que se pise a sí misma. Está demostrado que el mejor filtro en recepción es el mismo que se aplicase en transmisión, esto es, uno que responda con exactamente la forma de onda del pulso transmitido. Consideramos el canal de transmisión hasta la propia demodulación de fase y de nuestra señal ruidosa obtenemos el pulso patrón promediando todos los recibidos durante un minuto, de lo que resulta esto:

Ese pulso debería representar la excursión de frecuencia, o rampa lineal de fase, si bien tiene algunas diferencias con la forma esperada. En lugar de ver escalones planos de 25, 50 y 25 ms observamos curvas más suavizadas y la duración total se alarga hasta 114 ms. Eso puede deberse a cierto filtrado del equipo de radio y otros procesos. Pero sin embargo la zona central tiene un acusado hundimiento que no sabía explicar. Se puede obtener una demodulación precisa de la fase, que no de la frecuencia, mediante un proceso más general pero también más delicado y más intensivo en uso de cpu, que consiste en multiplicar la secuencia por un seno y un coseno, quedarnos con la componente DC y calcular el arco tangente del cociente:

corre_I=graba'.*seno;
corre_Q=graba'.*coseno;
% Filtro Butterworth 100 Hz
Wc=200/FS;
[coef_b,coef_a]=butter(4,Wc);
corre_filt_I=filter(coef_b,coef_a,corre_I);
corre_filt_Q=filter(coef_b,coef_a,corre_Q);
fase=atan(corre_filt_I./corre_filt_Q);

Representamos la señal así demodulada:

El pulso dura 100 ms y alcanza excursiones de fase de +1 y -1 radianes. Se parece bastante al pulso ideal publicado en la wikipedia salvo porque en el instante de tiempo 0 parece que la rampa se detiene y se aplana brevemente. Eso explicaría el hundimiento en el pulso obtenido como demodulación de frecuencia (que se corresponde con la derivada del pulso obtenido como demodulación de fase).

Por lo tanto la forma de onda del pulso patrón se puede tomar por buena y se puede emplear para filtrar la señal original. Esto hace que el resultado de ese filtrado sea, además, la correlación de la señal recibida con la esperada, y los valores pico se correspondan con los instantes óptimos de muestreo y decisión de los símbolos transmitidos. Filtramos no con el pulso patrón sino con la inversión temporal para que la correlación sea máxima. Daría lo mismo si fuese simétrico. Es decir, tras filtrar la señal, sólo nos queda sincronizarnos adecuadamente y limitarnos a observar los valores de la secuencia en los instantes temporales adecuados, para saber si se transmitió un bit 0 o 1. Tal como dice la wikipedia, y como se puede observar en la siguiente imagen, el bit 0 se señaliza con un único pulso mientras que el bit 1 contiene dos pulsos iguales seguidos.

Para localizar el instante preciso en el que se debe observar el valor de los bits transmitidos parece razonable empezar localizando el segundo 59, caracterizado por no tener modulación. Esto se puede conseguir con la función «xcorr», que muestra la correlación cruzada entre dos señales. Una de ellas es la señal de trabajo con la fase demodulada, pero rectificada para observar sólo el valor absoluto, y otra una simple serie de 1 segundo (500 muestras) de duración con valor constante 1. Cuando hay modulación de fase la señal tiene variaciones y por tanto una envolvente mucho más destacada que durante el segundo 59, donde sólo queda algo de ruido. Observamos claramente unos picos correspondientes al segundo 59:

Las marcas de los minutos se repiten, evidentemente, cada minuto. Para afinar la sincronización, mejorar la inmunidad al ruido y aprovechar toda la información disponible, podemos tomar la señal de la gráfica anterior en secuencias de 1 minuto de duración y sumarlas superponíendolas entre sí. Tenemos dos minutos enteros y parte de otro, por lo que ese último fragmento se rellena con el valor promedio. Esta operación realza el pico. La localización del inicio del minuto es simplemente hallar el instante temporal con el valor mínimo:

Una vez hecha esa sincronización gruesa, falta localizar la posición precisa del pulso inicial. Sabemos aproximadamente dónde debe estar, y nos dicen que los trece primeros segundos (0 al 12) tienen sólo el primer pulso, puesto que transmiten siempre un bit 0. Por lo tanto observamos unas ventanas que deben englobar el primer pulso de cada uno de los trece segundos, y las sumamos entre sí. El instante óptimo de muestreo es aquel en el que se obtenga el valor máximo:

Una vez estamos perfectamente sincronizados, podemos observar todos los segundos en los instantes del pulso de bit 0 y bit 1. Si el bit es 1 se transmite un segundo pulso 100 ms después del primero. Como no conocemos a priori la amplitud que debe alcanzar, he optado por comparar la amplitud del hipotético segundo pulso con la del primero: si es mayor que la mitad interpretamos que hay un pulso, y por tanto el bit es 1, si es menor interpretamos que es sólo ruido, y por tanto el bit es 0.

Finalmente se obtienen los bits de los 59 segundos y el script devuelve la información en un formato legible:

Procesamiento de la señal horaria TDF de Allouis
Número de muestras: 1100729. Frecuencia de muestreo: 8000 Hz
Duración de la grabación 137 s
SNR inicial 17
La grabación contiene 2 minutos íntegros
Procesamiento del minuto 1
Instante de muestreo (-125..+125):-4
SNR: 6
Bits iniciales segundos 0 a 12: 0 0 0 1 1 1 0 0 0 0 0 0 0
Festivo mañana: 0
Festivo hoy: 0
Operación anormal: 0
Horario de verano: 0
Horario de invierno: 1
Anuncio leap: 0
Inicio de hora: 1
Paridad par minutos: 0
Paridad par horas: 0
Hora: 19:44
Paridad fecha: 0
Día del mes: 10
Día de la semana: 5
Mes: 2
Año: 17

A pesar de lo publicado en la wikipedia, no todos los bits de los segundos 0 al 12 valen 0, sino que algunos valen 1. En la DCF lleva info meteorológica o protección civil.

El script realiza una serie de verificaciones sobre la señal para detectar si los datos de partida no son válido y reportar el error adecuado:

• Que la grabación dure al menos 1 minuto
• Si dura más de 10 minutos trunca la duración y se queda sólo con los 10 primeros
• Que la frecuencia de muestreo sea de 8 kHz
• Que haya un tono en 1 kHz
• Que se pueda sincronizar el inicio del minuto
• Que la grabación contenga una secuencia íntegra de un minuto
• Que se pueda sincronizar el símbolo

Además, la información original va acompañada de bits de paridad, por lo que es posible saber si ha habido un error en algún bit.

El script se puede descargar aquí. También hay que copiar el fichero que contiene la secuencia patrón con la que correlar los pulsos (aquí). Hay que ejecutar octave y, desde su línea de comandos, lanzar el script con «tdf1». Por si sirve de ayuda también comparto el fichero con la grabación. Para usar otro fichero diferente se debe editar la tercera línea del script para indicar el nombre del fichero wav que queramos utilizar. Debe grabarse con un receptor USB centrado en 161 kHz para que el tono suene a 1 kHz, y digitalizarse a 8 kHz.

Es posible demodular la señal en tiempo real con el programa «Clock» de F6CTE, autor también del Multipsk. Estos programas funcionan en Windows.

Comparto el script para su descarga a pesar de que todavía cabrían una serie de mejoras:

• Asegurarse de que 50 muestras después es el instante óptimo de muestreo del segundo pulso.
• Mostrar diagrama de ojo del pulso de bit 1 normalizado al valor del pulso de bit 0 a modo de AGC
• En caso de lsb multiplicar correlación por -1 o desplazar 2 muestras en sentido contrario. Idealmente, autodetectar usb/lsb.

Ensayos posteriores:

• Hasta qué valores de calidad medidos es buena la decodificación.
• Monitorizar propagación a lo largo de un día entero. Programar grabaciones a cada media hora. Apuntar SNR y graficar. Calcular la elevación del sol sobre el punto medio del camino.
• Como mejora del programa, demodular con pll en vez de con la propia señal desplazada, para reducir ruido.
• Investigar bits de datos.
• Comparar frecuencia de muestreo de la tarjeta de sonido correlando pulsos de segundo 0 de minutos consecutivos antes de submuestrear. Son 480.000 muestras o sea detecta hasta 2 ppm de error.
• Si la correlación con desplazamiento de 2 muestras tiene valor medio es que el receptor no entrega exactamente 1 kHz. Serviría para detectar un desajuste en su referencia de reloj.
• Señal Droitwich de 198 kHz.