Modem QAM: módulo RX
Modem QAM: módulo RX
Esta é a segunda parte do artigo que implementa um modem QAM em linguagem Python. Nesta parte revisaremos a implementação do receptor (RX).
Decodificar um sinal QAM de forma confiável é muito mais complicado que gerá-lo. Por conta disso, o RX tem muito mais código que o módulo TX, e devo explicar o código do RX parte por parte.
Primeiro, algum código burocrático:
#!/usr/bin/env python
import wave, struct, math, numpy, random
SAMPLE_RATE = 44100.0 # Hz
CARRIER = 1800.0 # Hz
BAUD_RATE = 360.0 # Hz
PHASES = 8
AMPLITUDES = 4
BITS_PER_SYMBOL = 5
# FIR lowpass filter. This time we calculated the coefficients
# directly, without resorting to FFT.
lowpass = []
CUTOFF = BAUD_RATE * 4
w = int(SAMPLE_RATE / CUTOFF)
R = (SAMPLE_RATE / 2.0 / CUTOFF)
for x in range(-w, w):
if x == 0:
y = 1
else:
y = math.sin(x * math.pi / R) / (x * math.pi / R)
lowpass.append(y / R)
# Modem constellation (the same as TX side)
constellation = {}
invconstellation = {}
for p in range(0, PHASES):
invconstellation[PHASES - p - 1] = {}
for a in range(0, AMPLITUDES):
s = p * AMPLITUDES + a
constellation[s] = (PHASES - p - 1, a + 1)
invconstellation[PHASES - p - 1][a + 1] = s
qam = wave.open("qam.wav", "r")
n = qam.getnframes()
qam = struct.unpack('%dh' % n, qam.readframes(n))
qam = [ sample / 32768.0 for sample in qam ]
Parte do código foi emprestado do TX: parâmetros de modulação, gerador de constelação, etc. Estas coisas deveriam estar num módulo comum separado, mas (à época em que isto foi escrito) estava com preguiça de fazer esse refatoramento. O loop de constelação também gera um dicionário reverso, já que estamos fazendo o contrário do TX (dada amplitude e fase, me diga os bits).
Um item novo é o gerador de filtro FIR. Poderia ter usado FFT para gerar os coeficientes, mas aqui preferi usar uma fórmula pronta. Tenho ciência que este filtro FIR não é o melhor, mas serviu bem o suficiente para esta aplicação.
# Demodulate in "real" and "imaginary" parts. The "real" part
# is the correlation with carrier. The "imaginary" part is the
# correlation with carrier offset 90 degrees. Having both values
# allows us to find QAM signal's phase.
real = []
imag = []
amplitude = []
# This proofs that we don't need to replicate the original carrier
# (at least, not in the same phase). Sometimes this random carrier
# phase makes the receiver to interpret training sequence as a
# character, but nothing serious.
t = int(random.random() * CARRIER)
# t = -1
# In the other hand, frequency must match carrier exactly; if we
# need to allow for frequency deviations, we must build a digital
# equivalent to PLL / VCO oscilator.
for sample in qam:
t += 1
angle = CARRIER * t * 2 * math.pi / SAMPLE_RATE
real.append(sample * math.cos(angle))
imag.append(sample * -math.sin(angle))
amplitude.append(abs(sample))
del qam
O bloco acima constitui a parte "analógica" do decodificador, que é muito simples de implementar, porque a trigonometria faz sozinha toda a mágica.
Com a modulação AM analógica nós aprendemos que um sinal demodulado tem dois componentes: a informação original que nós queremos, e uma remodulação em alta freqüência (em torno do dobro da portadora) que precisamos descartar.
A decodificação QAM produz dois resultados brutos: as partes "real" e "imaginária". Ambas são contaminadas pelo ruído da remodulação. Se você analisar a forma de onda delas, verá que não fazem sentido. Ambas precisam de filtragem.
Um pouco de matemática ajuda a reforçar esse ponto, que se aplica tanto a AM quanto a QAM. A informação original f(t) é modulada pela função cos(p.t), e o sinal transmitido é igual a f(t).cos(p.t). A demodulação converte isso para f(t).cos(p.t).cos(p.t), que é equivalente a f(t)+[f(t).cos(2.p.t)]. Esta última fórmula ilustra que o sinal demodulado é soma de duas partes: o sinal f(t) desejado, e o indesejável sinal remodulado por cos(2.p.t).
Já que fizemos uma espécie de modulação AM no transmissor, multiplicando um sinal digital por uma portadora cossenóide, fazemos o mesmo do lado RX: multiplicar de novo pela mesma portadora para obter o sinal digital original. Não sabemos exatamente a fase da portadora QAM, nem precisamos saber. Apenas multiplicamos por duas portadoras locais em quadratura (math.cos() e math.sin()) idênticas exceto pela defasagem de 90 graus.
Alguns leitores ficaram um pouco confusos nessa parte: como podemos extrair dois valores (real e imaginário) de um único sinal? Ou, como podemos medir a fase de um sinal apenas multiplicando-o por cos() and -sin()?
Esta é realmente uma diferença fundamental entre AM e QAM. Em AM, a portadora nunca muda de fase; apenas a variação de amplitude interessa; e apenas um sinal pode ser transmitido. Em QAM, amplitude e fase são ambas manipuladas, de modo que dois sinais podem ser transmitidos ao mesmo tempo. Métodos de modulação como AM e PSK (codificação por mudanças de fase) são simplificações de QAM onde fase ou amplitude nunca mudam, respectivamente.
Tanto modulação quanto demodulação usam duas portadoras ortogonais de uma mesma freqüência: cos() e -sin(). No modulador, manipulamos a amplitude de cada portadora para gerar um sinal QAM de qualquer fase (e amplitude) desejada. No demodulador, a trigonometria nos permite reutilizar o mesmo truque.
Mais um pouco de matemática pode convencer que o truque funciona na demodulação. Suponha que a informação original seja: real=0, imaginário=-1. Isto produz um sinal QAM igual a sin(p.t). No lado RX, este sinal é multiplicado pelas portadoras em quadratura:
QAM = sin(p.t) real = sin(p.t) * cos(p.t) = 0 + sin(2.p.t) imag = sin(p.t) * -sin(p.t) = -1 + cos(2.p.t) Depois de remover por filtragem a remodulação em 2pt: real = 0 imag = -1
Há muitas maneiras de entender como a quadradura QAM realmente funciona. Ela pode ser vista como dois sinais AM simultâneos defasados de 90 graus, ou pode ser encarada como uma mídia capaz de transmitir números complexos. Ou, como eu disse antes, notar que duas grandezas físicas estão sendo manipuladas (amplitude e fase) e portanto podemos transmitir um sinal através de cada grandeza.
É claro, no pequeno exercício matemático acima, nós assumimos que o modulador e o demodulador tenham suas portadoras perfeitamente em fase, o que é impossível na prática. Por conta disso, o receptor QAM procura por mudanças na fase, não pela fase absoluta, que será diferente em cada lado.
Note esta linha de código, retirada do trecho de código exibido antes:
t = int(random.random() * CARRIER)
A variável t ganha um valor aleatório inicial, o que significa que a fase da portadora também é aleatória. Isto assegura que a portadora RX quase nunca estará em fase com TX, e mesmo assim RX consegue receber dados, e fica provado que usar fase relativa em vez de absoluta pode funcionar.
A saída do demodulador é nada mais que a posição cartesiana X,Y do símbolo na constelação. As posições computadas estarão rotacionadas, porque a portadora RX não está em fase com TX, mas como nós integramos as alterações de fase no bloco TX, só precisamos considerar a diferença de fase de cada símbolo em relação ao anterior.
Finalmente, nós "demodulamos" a amplitude num canal separado, usando abs(), por pura preguiça, já que poderíamos calculá-la a partir das partes real e imaginária.
Tudo parece fácil até agora porque não tivemos de lidar com distorções de freqüência da portadora. Na prática, a portadora RX teria de "afinar-se" com o sinal recebido, porque a demodulação AM, AM-SSB e QAM dependem de uma portadora com exatamente a mesma freqüência da portadora TX. E todos os meios de transmissão do mundo real distorcem freqüência em algum grau. (Nossos arquivos WAV não distorcem, o que simplifica nossos exemplos.)
# Pass all signals through lowpass filter
real = numpy.convolve(real, lowpass)
imag = numpy.convolve(imag, lowpass)
amplitude = numpy.convolve(amplitude, lowpass)
# After lowpass filter, all three lists show something like
# a square wave, pretty much like the original bitstream.
# If you want to see the result of early detection phase,
# uncomment this:
# f = wave.open("teste.wav", "w")
# f.setnchannels(1)
# f.setsampwidth(2)
# f.setframerate(44100)
# bla = [ int(sample * 32767) for sample in imag ]
# f.writeframes(struct.pack('%dh' % len(bla), *bla))
Como dito antes, a demodulação AM e QAM deixam um resíduo de alta freqüência, então passamos nosso filtro FIR em todos os sinais (real, imaginário e amplitude). Depois dessa filtragem, os sinais são ondas quadradas praticamente perfeitas:
O sinal filtrado ainda tem "rebarbas" porque, quando o símbolo muda, o sinal QAM sofre uma mudança abrupta de fase e amplitude. Essa mudança brusca injeta ruído em todas as freqüências (assim como um raio gera todo tipo de ruído, desde o trovão audível até interferência eletromagnética). O filtro passa-baixas remove a maior parte, mas não tudo, desse ruído.
Esta é a razão pelo qual o filtro passa-baixas precisa ser seletivo em relação a freqüências acima da taxa de baud. Num primeiro olhar, um filtro mais simples e mais rápido que remova o ruído de demodulação parece bom o suficiente, mas ele deixaria muito ruído inter-símbolos passar batido.
As coisas talvez sejam diferentes em QAM analógico. Sim, além do QAM digital existe o QAM analógico — ele é utilizado para transmitir informação de cor na TV analógica (que aliás está em vias de extinção). Em QAM analógico, a informação não possui um "baud rate" definido e não há saltos abruptos, então a filtragem pós-modulação pode ser mais simples.
Também é importante mencionar a taxa de amostragem neste ponto. A taxa de amostragem precisa ser alta o suficiente para conter o ruído de demodulação AM. Por exemplo, se a portadora é 1800Hz, o ruído de demodulação ficará centrado em 3600Hz. Se a taxa de baud for de 600Hz, o sinal QAM modulado cobre a banda 1200-2400Hz, e o ruído de demodulação fica na banda 3000-4200Hz. A taxa de amostragem precisa ser de, no mínimo, 4200x2=8400Hz.
O senso comum nos diz que, se utilizarmos uma taxa de amostragem apenas um pouco superior à portadora, o ruído de demodulação acima do limite de Nyquist é simplesmente deixado de fora. Infelizmente não é assim que a coisa funciona: sinal (ou ruído) acima do limite de Nyquist sofrerá aliasing ou "dobramento". Por exemplo, amostrar um sinal de 3500Hz com uma taxa de 3600Hz (cujo limite de Nyquist é 1800Hz) fará aparecer um sinal de 100Hz nas amostras.
Isto significa que, se você tentar usar uma taxa de amostragem de 4000Hz para demodular um sinal QAM cuja portadora é 1800Hz, o ruído de demodulação centrado em 3600Hz será "dobrado" para a faixa dos 400Hz (4000-1800x2), misturando-se com a mensagem e provavelmente destruindo-a.
Então, para fechar o assunto, é necessário usar uma taxa de amostragem relativamente alta, e é preciso empregar alguma espécie de filtro passa-baixas para condicionar a mensagem recebida. Não há atalhos.
De volta para o gráfico dos sinais demodulados, depois de filtrar:
Temos três ondas quadradas semelhantes à figura acima: sinal real, sinal imaginário e amplitude. Agora a missão é extrair os símbolos dessa trinca. Primeiro, precisamos calcular a fase dos símbolos:
# Detect phase based on real and imaginary values
phase = []
wrap = 2.0 * math.pi * (PHASES - 0.5) / PHASES
for t in range(0, len(real)):
# This converts (real, imag) to an angle
angle = math.atan2(imag[t], real[t])
if angle < 0:
angle += 2 * math.pi
# Move angle near to 2 pi to 0
if angle > wrap:
angle = 0.0
phase.append(angle)
Uma vez que nossa constelação não emprega números complexos, as coordenadas X,Y encontradas pelo demodulador não são diretamente utilizáveis na procura pelo símbolo, então precisamos realmente calcular o ângulo. O coração do cálculo é a função atan2().
# Find maximum amplitude based on training whistle (1s), and
# measure how much our local carrier is out-of-phase in
# relationship to QAM signal
max_amplitude = 0.0
carrier_real = 0.0
carrier_imag = 0.0
for t in range(int(0.1 * SAMPLE_RATE), int(0.9 * SAMPLE_RATE)):
max_amplitude += amplitude[t]
carrier_real += real[t]
carrier_imag += imag[t]
max_amplitude /= int(0.8 * SAMPLE_RATE)
carrier_real /= int(0.8 * SAMPLE_RATE)
carrier_imag /= int(0.8 * SAMPLE_RATE)
skew = math.atan2(carrier_imag, carrier_real)
print "Carrier phase difference: %d degrees" % (skew * 180 / math.pi)
del imag
del real
No código acima nós analisamos o assobio de "treinamento" enviado pelo TX, para determinar a faixa dinâmica de amplitude, assim como a diferença de fase entre nossa portadora e a portadora TX. Fazemos isso de forma muito preguiçosa aqui, esta implementação não seria boa o suficiente para um modem comercial.
Qualquer meio de transmissão real, seja por fio, fibra ótica ou rádio, distorce amplitude, fase e freqüência. Assim, não seria suficiente estimar a amplitude apenas no início de uma sessão. Modems comerciais usam PLL para gerar a portadora local e mantê-la ajustada de acordo com o sinal recebido.
Mas, para nossa sorte, arquivos WAV não nos trazem esse tipo de problema :)
# Normalize/quantify amplitude and phase to constellation steps
qsymbols = []
for t in range(0, len(phase)):
p = phase[t]
a = amplitude[t]
a /= max_amplitude
a = int(a * AMPLITUDES + 0.49) # adding 0.49 is a rounding trick
# Compensate for local carrier phase difference
# This ensures the best phase quantization (avoiding that
# quantization edge is too near a constellation point)
p += 2 * math.pi - skew
p /= 2 * math.pi
p = int(p * PHASES + 0.49) % PHASES
qsymbols.append((p, a))
del phase
del amplitude
No código acima, nós pegamos os valores "analógicos" de fase e amplitude, e quantizamos os mesmos dentro da constelação. Neste ponto, também corrigimos a fase, para que os valores quantificados estejam alinhados com o transmissor e com a constelação.
O problema agora é: temos milhares de amostras que contém apenas um punhado de símbolos. Precisamos detectar quantos símbolos existem, e onde eles estão.
# Try to detect edges between symbols
edge = []
settle_time = int(SAMPLE_RATE / BAUD_RATE / 8)
settle_timer = -1
in_symbol = 0
first_symbol = 0
last_a = last_p = 0
for t in range(0, len(qsymbols)):
s = 0
p, a = qsymbols[t]
if in_symbol > 0:
# we presume we are flying over a valid symbol
in_symbol -= 1
elif a != last_a or p != last_p:
# change detected, look for stabilization in future
settle_timer = settle_time
elif settle_timer > 0:
# one sample closer a good symbol
settle_timer -= 1
elif settle_timer == 0:
# signal settled for (settle_time) samples;
# take the current signal as a valid symbol
settle_timer = -1
in_symbol = int(0.85 * SAMPLE_RATE / BAUD_RATE)
if t > 0.5 * SAMPLE_RATE:
# make sure we are not at the beginning
s = 1
if not first_symbol:
first_symbol = t
edge.append(s)
last_p = p
last_a = a
O trecho de código acima conta com o fato que existe um "degrau" abrupto entre um símbolo e outro, na maioria dos casos. Nós anotamos os pontos onde esse degrau existe, para referência futura.
Esse passo não encontra todos os símbolos, porque a mensagem pode ter símbolos repetidos adjacentes, entre os quais não há "degrau". Isto é resultado direto da má escolha da constelação: teria ajudado evitar pontos com 0º de fase, e todos os símbolos seriam marcados por degraus.
# Find symbols based on edges and known baud rate
symbols = []
expected_next_symbol = SAMPLE_RATE / BAUD_RATE
lost = 0
for t in range(first_symbol - int(SAMPLE_RATE * 0.1), len(qsymbols)):
p, a = qsymbols[t]
s = edge[t]
if s:
# amplitude/phase edge, strong hint for symbol
lost = 0
symbols.append([p, a])
else:
lost += 1
if lost > expected_next_symbol * 1.5:
# no edge yet, take this as a non-transition signal
lost -= expected_next_symbol
symbols.append([p, a])
del qsymbols
del edge
No código acima, terminamos de detectar os símbolos em que estamos interessados, e podemos descartar uma parte das informações coletadas.
A estratégia do detector de símbolos é muito smples. Cada "degrau" detectado com certeza precede um novo símbolo. Por outro lado, se nenhum "degrau" é encontrado num certo intervalo (baud rate mais 50% de tolerância), assumimos que temos um símbolo repetido.
Mesmo que as taxas de baud TX e RX fossem um pouco diferentes, a estratégia acima seria robusta, poirque a maior parte dos símbolos é precedida por "degraus" e o algoritmo não trabalha às cegas por mais que 2 ou 3 símbolos, assumindo que todos os símbolos são igualmente prováveis.
print "%d symbols found" % len(symbols)
# Differentiate phase (because we had integrated it at transmitter)
last_phase = 0
for t in range(0, len(symbols)):
p = symbols[t][0]
symbols[t][0] = (p - last_phase + PHASES) % PHASES
last_phase = p
Já que tínhamos integrado (somado) as mudanças de fase, agora nós fazemos a operação oposta, que é diferenciar.
# Translate constellation points into bit sequences
bitgroups = []
for t in range(0, len(symbols)):
p, a = symbols[t]
try:
bitgroups.append(invconstellation[p][a])
except KeyError:
print "Bad symbol:", p, a, ", ignored"
# Translate bit groups into individual bits
bits = []
for bitgroup in bitgroups:
for bit in range(0, BITS_PER_SYMBOL):
bits.append((bitgroup >> (BITS_PER_SYMBOL - bit - 1)) % 2)
Finalmente, os símbolos compostos por amplitude e fase podem ser transformados em seqüências de bits. A demodulação propriamente dita está feita. Agora precisamos recuperar os bytes a partir dos bits, semelhantemente a um receptor RS-232.
# Find bytes based on start and stop bits. This will make you feel
# like a poor RS-232 port :)
t = 0
state = 0
byte = []
bytes = []
while t < len(bits):
bit = bits[t]
if state == 0:
# did not find start bit
if not bit:
# just found start bit
state = 1
elif state == 1:
# found start bit, accumulating a byte
if len(byte) < 8:
byte.append(bit)
else:
# byte complete, test stop bit
if bit:
# stop bit is there
bytes.append(byte)
else:
# oops, stop bit not there, we were cheated!
# backtrack to the point fake start bit was 'found'
t -= 8
byte = []
state = 0
t += 1
# Make a string from the bytes
msg = ""
print "Bytes:",
for byte in bytes:
value = 0
for bit in range(0, 8):
value += byte[bit] << (8 - bit - 1)
msg += chr(value)
print value,
print
print msg
O código procura por start bits (0's) e stop bits (1's). Encontrada uma seqüência de bits que satisfaça essas condições, um byte é recriado a partir dela. O algoritmo pode ser enganado por ruídos, e acabar achando um ou outro byte onde não há nenhum. Mas a tendência é ressincronizar rapidamente e funcionar bem.
Outros artigos dessa mesma série discutem e implementam diversas melhorias no modem QAM Python, como o uso de PLL para sincronizar a portadora RX, o uso de randomizador para garantir que todos os símbolos sejam igualmente prováveis, a escolha de uma constelação melhor, e também o uso de técnicas mais avançadas de demodulação, mais apropriadas a um modem implementado digitalmente e menos imitativas do demodulador AM analógico.