Recebendo dados de baloes meteorologicos da Aeronautica

* This post is in portuguese

Hoje finalmente consegui rastrear um dos balões meteorológicos que a aeronáutica lança duas vezes por dia aqui em Porto Alegre / RS. A aeronáutica utiliza as sondas da Vaisala (uma empresa finlandesa) modelo RS-92SGP para realizar as medições de umidade, temperatura e pressão. Estes dados são geralmente utilizados para as previsões de tempo da região; existe um datasheet com mais dados sobre o equipamento que eles utilizam, neste datasheet tem informações importantes como por exemplo a frequência em que o aparelho envia os dados de telemetria. Aqui em Porto Alegre / RS a aeronáutica está utilizando a faixa de operação em 402.700Mhz, que também é coberta pelos dongles USB RTLSDR como o que eu utilizo.

O equipamento da Vaisala é um equipamento que utiliza 60mW de potência na transmissão (eu já vi balões transmitindo até 600km nessa frequência com apenas 10mW e com uma antena decente é claro) e utiliza modulação GFSK. Para decodificar o protocolo e GFSK podemos utilizar  SDR# juntamente com o Virtual Cable (ou algo semelhante para redirecionar os dados do SDR# para o SondeMonitor que é o software que irá fazer a decodificação dos dados (infelizmente o software é pag e só roda apenas em Windows, mas ao menos vem com alguns dias de trial).

No meu setup eu estou utilizando um dongle RTLSDR R820T juntamente com um Low Noise Amplifier (LNA4ALL na foto abaixo) e uma antena de 5dB omnidirecional:

LNA e regulador de tensão

LNA e regulador de tensão

Abaixo segue a foto do equipamento lançado como payload do balão meteorológico:

Sonda RS92-SGP

Sonda RS92-SGP

Estes balões geralmente atingem uma altitude de uns 20km a 35km, mas isto depende de vários fatores como por exemplo os ventos, a quantidade de gás que foi utilizada no balão, a espessura do latex do balão e outros fatores. Quando o balão estoura este fenômeno é geralmente chamado de “burst” e após este estouro o balão acaba caindo por terra (ele tem uma bateria que não agride a natureza).

Seguem abaixo os screenshots do recebimento dos dados, neste momento eu ainda não havia conseguido receber toda calibração do aparelho:

SDR# recebendo o sinal da sonda

SDR# recebendo o sinal da sonda

Screenshot de alguém mais fazendo o tracking da sonda e jogando para o APRS aqui de Porto Alegre:

Localização do balão no APRS

Localização do balão no APRS

Imagem do rastremento do balão no APRS:

Rastreamento da sonda

O próximo passo agora é conseguir uma antena direcional para melhorar a recepção =)

Para quem tiver interesse em receber os dados, os balões são lançados diariamente as 00:00 UTC e às 12:00 UTC.

- Christian S. Perone

Simple and effective coin segmentation using Python and OpenCV

The new generation of OpenCV bindings for Python are getting better and better with the hard work of the community. The new bindings, called “cv2″ are the replacement of the old “cv” bindings; in this new generation of bindings, almost all operations returns now native Python objects or Numpy objects, which is pretty nice since it simplified a lot and also improved performance on some areas due to the fact that you can now also use the optimized operations from Numpy and also enabled the integration with other frameworks like the scikit-image which also uses Numpy arrays for image representation.

In this example, I’ll show how to segment coins present in images or even real-time video capture with a simple approach using thresholding, morphological operators and contour approximation. This approach is a lot simpler than the approach using Otsu’s thresholding and Watershed segmentation here in OpenCV Python tutorials, which I highly recommend you to read due to its robustness. Unfortunatelly, the approach using Otsu’s thresholding is highly dependent on a illumination normalization. One could extract small patches of the image to implement something similar to an adaptive Otsu’s binarization (like the one implemented in Letptonica – the framework used by Tesseract OCR) to overcome this problem, but let’s see another approach. For reference, see the output of the Otsu’s thresholding using an image taken with my webcam with a non-normalized illumination:

 

Original image vs Otsu binarization

Original image vs Otsu binarization

1. Setting the Video Capture configuration

The first step to create a real-time Video Capture using the Python bindings is to instantiate the VideoCapture class, set the properties and then start reading frames from the camera:

1
2
3
4
5
6
import numpy as np
import cv2

cap = cv2.VideoCapture(0)
cap.set(cv2.cv.CV_CAP_PROP_FRAME_WIDTH, 1280)
cap.set(cv2.cv.CV_CAP_PROP_FRAME_HEIGHT, 720)

In newer versions (unreleased yet), the constants for CV_CAP_PROP_FRAME_WIDTH are now in the cv2 module, for now let’s just use the cv2.cv module.

2. Reading image frames

The next step is to use the VideoCapture object to read the frames and then convert them to gray color (we are not going to use color information to segment the coins):

1
2
3
4
while True:
    ret, frame = cap.read()
    roi = frame[0:500, 0:500]
    gray = cv2.cvtColor(roi, cv2.COLOR_BGR2GRAY)

Note that here I’m extracting a smal portion of the complete image (where the coins are located), but you don’t have to do that if you have only coins on your image. At this moment, we have the following gray image:

The original Gray image captured.

The original Gray image captured.

3. Applying adaptive thresholding

In this step we will apply the Adaptive Thresholding after applying a Gaussian Blur kernel to eliminate the noise that we have in the image:

1
2
3
gray_blur = cv2.GaussianBlur(gray, (15, 15), 0)
thresh = cv2.adaptiveThreshold(gray_blur, 255, cv2.ADAPTIVE_THRESH_GAUSSIAN_C,
    cv2.THRESH_BINARY_INV, 11, 1)

See the effect of the Gaussian Kernel in the image:

The original gray image and the image after applying the Gaussian Kernel.

The original gray image and the image after applying the Gaussian Kernel.

And now the effect of the Adaptive Thresholding with the blurry image:

The effect of the adaptive thresholding into the blurry image

Note that at that moment we already have the coins segmented except for the small noisy inside the center of the coins and also in some places around them.

4. Morphology

The Morphological Operators are used to dilate, erode and other operations on the pixels of the image. Here, due to the fact that sometimes the camera can present some artifacts we will use the Morphological Operation of Closing to make sure that the borders of the coins are always close, otherwise we may found a coin with a semi circle or something like that. To understand the effect of the Closing operation (which is the operation of erosion of the pixels already dilated) see the image below:

Morphological Closing

You can see that after some iterations of the operation, the circles starts to become filled. To use the Closing operation, we’ll use the morphologyEx function from the OpenCV Python bindings:

1
2
kernel = np.ones((3, 3), np.uint8)
closing = cv2.morphologyEx(thresh, cv2.MORPH_CLOSE, kernel, iterations=4)

See now the effect of the Closing operation on our coins:

Closing Operation in the coins

The operations of Morphological Operators are very simple, the main principle is the application of a element (in our case we have a block element of 3×3) into the pixels of the image. If you want to understand it, please see this animation explaning the operation of Erosion.

5. Contour detection and filtering

After applying the morphological operators, all we have to do is to find the contour of each coin and then filter the contours having an area smaller or larger than a coin area. You can imagine the procedure of finding contours in OpenCV as the operation of finding connected components and their boundaries. To do that, we’ll use the OpenCV findContours function.

1
2
3
cont_img = closing.copy()
contours, hierarchy = cv2.findContours(cont_img, cv2.RETR_EXTERNAL,
    cv2.CHAIN_APPROX_SIMPLE)

Note that we made a copy of the closing image because the function findContours will change the image passed as the first parameter, we’re also using the RETR_EXTERNAL flag, which means that the contours returned are only the extreme outer contours. The parameter CHAIN_APPROX_SIMPLE will also return a compact representation of the contour, for more information see here.

After finding the contours, we need to iterate into each one and check the area of them to filter the contours containing an area greater or smaller than the area of a coin. We also need to fit an ellipse to the contour found. We could have done this using the minimum enclosing circle, but since my cameara isn’t perfectly above the coins, the coins appear with a small inclination describing an ellipse.

1
2
3
4
5
6
7
8
9
10
for cnt in contours:
    area = cv2.contourArea(cnt)
    if area < 2000 or area > 4000:
        continue

    if len(cnt) < 5:
        continue

    ellipse = cv2.fitEllipse(cnt)
    cv2.ellipse(roi, ellipse, (0,255,0), 2)

Note that in the code above we are iterating on each contour, filtering coins with area smaller than 2000 or greater than 4000 (these are hardcoded values I found for the brazilian coins at this distance from the camera), later we check for the number of points of the contour because the function fitEllipse needs a number of points greater or equal than 5 and finally we use the ellipse function to draw the ellipse in green over the original image.

To show the final image with the contours we just use the imshow function to show a new window with the image:

1
cv2.imshow('final result', roi)

And finally, this is the result in the end of all steps described above:

The final image with the contours detected

The final image with the contours detected

The complete source-code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
import numpy as np
import cv2

def run_main():
    cap = cv2.VideoCapture(0)
    cap.set(cv2.cv.CV_CAP_PROP_FRAME_WIDTH, 1280)
    cap.set(cv2.cv.CV_CAP_PROP_FRAME_HEIGHT, 720)

    while(True):
        ret, frame = cap.read()
        roi = frame[0:500, 0:500]
        gray = cv2.cvtColor(roi, cv2.COLOR_BGR2GRAY)

        gray_blur = cv2.GaussianBlur(gray, (15, 15), 0)
        thresh = cv2.adaptiveThreshold(gray_blur, 255, cv2.ADAPTIVE_THRESH_GAUSSIAN_C,
            cv2.THRESH_BINARY_INV, 11, 1)

        kernel = np.ones((3, 3), np.uint8)
        closing = cv2.morphologyEx(thresh, cv2.MORPH_CLOSE,
            kernel, iterations=4)

        cont_img = closing.copy()
        contours, hierarchy = cv2.findContours(cont_img, cv2.RETR_EXTERNAL,
            cv2.CHAIN_APPROX_SIMPLE)
       
        for cnt in contours:
            area = cv2.contourArea(cnt)
            if area < 2000 or area > 4000:
                continue

            if len(cnt) < 5:
                continue

            ellipse = cv2.fitEllipse(cnt)
            cv2.ellipse(roi, ellipse, (0,255,0), 2)

        cv2.imshow("Morphological Closing", closing)
        cv2.imshow("Adaptive Thresholding", thresh)
        cv2.imshow('Contours', roi)

        if cv2.waitKey(1) & 0xFF == ord('q'):
            break

    cap.release()
    cv2.destroyAllWindows()

if __name__ == "__main__":
    run_main()

Despesas de Custeio e Lei de Benford

* This post is in Portuguese.

Há poucos dias, a prefeitura de Porto Alegre liberou os datasets com os dados de despesas de custeio de vários órgãos municipais (Secretaria Municipal de Saúde, Secretaria Municipal de Cultura, Gabinete do Prefeito, etc.).  O plot abaixo mostra a quantidade de empenhos para cada órgão municipal:

Plot - Qtd Empenhos vs Órgãos

Plot – Qtd Empenhos vs Órgãos

Uma das maneiras utilizadas geralmente para verificar fraudes é o uso da Lei de Benford [1] [2] [3], que fala sobre a distribuição das frequências de dígitos em vários datasets da vida real, incluindo valores de ações, número de populações, tamanhos de rios, etc.

Ao correlacionar a distribuição de números dos primeiros digitos dos valores de empenhos dos dados de Despesas de Custeio do 2º bimestre de 2014 com a distribuição da Lei de Benford,  a correlação ficou muito clara:

 

Lei de Benford vs Despesas de Custeio (Empenho)

Lei de Benford vs Despesas de Custeio (Empenho)

Segue aí mais um exemplo de correlação da Lei de Benford. Um sistema legal para ser construído seria um monitor de despesas que verificasse a correlação da Lei de Benford automaticamente e alertasse a cada anomalia encontrada.

Universality, primes and space communication

So, in mathematics we have the concept of universality in which we have laws like the law of large numbers, the Benford’s law (that I cited a lot in previous posts), the central limit theorem and many other laws that acts like laws of physics for the world of mathematics. These laws are not our inventions, I mean, the concepts are our inventions but the laws per se are universal, they are true no matter where you are on the earth or if you live far away on the universe. And that is why Frank Drake, one of the founders of SETI and also one of the pioneers in search for extraterrestrial intelligence came with this brilliant idea of using prime numbers (another example of universality) to communicate with distant worlds. The idea that Frank Drake had was the use of prime numbers to hide (not actually hide, but to make self evident, you’ll understand later) the dimension of a transmitted image in the image size itself.

So, imagine you are receiving a message that is a sequence of dashes and dots like “—.-.—.-.——–…-.—” that repeats after a short pause and then again and again. Let’s suppose that this message has the size of 1679 symbols. So you begin analyzing the number, which is in fact a semiprime number (the same used in cryptography, a number that is a product of two prime numbers) that can be factored in prime factors as 23*73=1679, and this is the only way to factor it in prime factors (actually all numbers have only a single set of prime factors that are unique, see Fundamental theorem of arithmetic). So, since there are only two prime factors, you will try to reshape the signal in a 2D image and this image can have the dimension of 23×73 or 73×23, when you arrange the image in one of these dimensions you’ll see that the image makes sense and the other will be just a random and strange sequence. By using prime numbers (or semiprimes) you just used the total image size to define the only two possible ways of arranging the image dimension.

Arecibo Radio Telescope

This idea was actually used in reality in 1974 by the Arecibo radio telescope when a message was broadcast in frequency modulation (FM) aiming the M13 globular star cluster at 25.000 light-years away:

M13 Globular Star Cluster

This message had the size (surprise) of 1679 binary digits and carried a lot of information of your world like: a graphical representation of an human, numbers from 1 to 10, a graphical representation of the Arecibo radio telescope, etc.

The message decoded as 23 rows and 73 columns is this:

Arecibo Message Shifted (Source: Wikipedia)

As you can see, the message looks a lot nonsensical, but when it is decoded as an image with 73 rows and 23 columns, it will show its real significance:

Arecibo Message with the correct dimension (Source: Wikipedia)

Amazing, don’t you think ? I hope you liked it !

- Christian S. Perone

 

The beauty of Bitcoin P2P network

So, in the last days I just released Protocoin, a framework in pure Python with a Bitcoin P2P network implementation. While I’m in process of development of the v.0.2 of the framework (with new and nice features like Bitcoin keys management – you can see some preview here) I would like to show a real-time visualization I’ve made with Protocoin and Ubigraph of a node connecting to a seed node and then issuing GetAddr message for each node and connecting on the received nodes in a breadth-first search fashion. I’ll release the code used to create this visualization in the next release of Protocoin as soon as possible. I hope you enjoy it !

Color legend

Yellow = Connecting
Green = Connected
Blue = Disconnected after connection

Video

Protocoin – a pure Python Bitcoin protocol implementation

Just release the first version of Protocoin, a pure Python Bitcoin protocol implementation, for more information see the documentation or the project in Github.

If you want to donate for the project, my Bitcoin address is: 1Q6JZEE5turJXaTn1fburWkqjhC4oMJ4yV.

I hope you like it !

- Christian S. Perone

Mapa de calor dos dados de acidentes de transito do DataPoa

Esta semana será disponibilzada a nova versão do Django GIS Brasil, segue abaixo um exemplo de mapa criado usando os novos dados do Django GIS Brasil importados do DataPoa.

O exemplo abaixo é um mapa de calor utilizando os dados de acidentes de trânsito em Porto Alegre /RS durante os anos de 2000 até 2012. Os eixos (ruas, avenidas, etc.) também estarão presentes no Django GIS Brasil.

 

Heatmap POA

Mapa de Acidentes de Trânsito em PoA/RS. Clique para ampliar.

Book Suggestion: Codex Seraphinianus

Today the Codex Seraphinianus just arrived (after months waiting in the pre-order state). I bought it from Amazon and I really recommend this edition for those who are interested because this is a very large edition with high quality textured paper and beautiful printing style. The book has also in the end a pocket with a small brochure called “Decodex” with a letter from Luigi Serafini.

The book is a very impressive creation by Luigi Serafini (or by the cat) dating from 1981 and presenting an impossible world that will cause to you the most strange feelings. See the photo of the cover and some pages below.

- Christian S. Perone

Codex Page

Codex Title

Machine Learning :: Cosine Similarity for Vector Space Models (Part III)

* It has been a long time since I wrote the TF-IDF tutorial (Part I and Part II) and as I promissed, here is the continuation of the tutorial. Unfortunately I had no time to fix the previous tutorials for the newer versions of the scikit-learn (sklearn) package nor to answer all the questions, but I hope to do that in a close future.

So, on the previous tutorials we learned how a document can be modeled in the Vector Space, how the TF-IDF transformation works and how the TF-IDF is calculated, now what we are going to learn is how to use a well-known similarity measure (Cosine Similarity) to calculate the similarity between different documents.

The Dot Product

Let’s begin with the definition of the dot product for two vectors: \vec{a} = (a_1, a_2, a_3, \ldots) and \vec{b} = (b_1, b_2, b_3, \ldots), where a_n and b_n are the components of the vector (features of the document, or TF-IDF values for each word of the document in our example) and the \mathit{n} is the dimension of the vectors:

  \vec{a} \cdot \vec{b} = \sum_{i=1}^n a_ib_i = a_1b_1 + a_2b_2 + \cdots + a_nb_n

As you can see, the definition of the dot product is a simple multiplication of each component from the both vectors added together. See an example of a dot product for two vectors with 2 dimensions each (2D):

  \vec{a} = (0, 3) \\   \vec{b} = (4, 0) \\   \vec{a} \cdot \vec{b} = 0*4 + 3*0 = 0

The first thing you probably noticed is that the result of a dot product between two vectors isn’t another vector but a single value, a scalar.

This is all very simple and easy to understand, but what is a dot product ? What is the intuitive idea behind it ? What does it mean to have a dot product of zero ? To understand it, we need to understand what is the geometric definition of the dot product:

  \vec{a} \cdot \vec{b} = \|\vec{a}\|\|\vec{b}\|\cos{\theta}

Rearranging the equation to understand it better using the commutative property, we have:

   \vec{a} \cdot \vec{b} = \|\vec{b}\|\|\vec{a}\|\cos{\theta}

So, what is the term \displaystyle \|\vec{a}\|\cos{\theta} ? This term is the projection of the vector \vec{a} into the vector \vec{b} as shown on the image below:

The projection of the vector A into the vector B. By Wikipedia.

Now, what happens when the vector \vec{a} is orthogonal (with an angle of 90 degrees) to the vector \vec{b} like on the image below ?

Two orthogonal vectors (with 90 degrees angle).

There will be no adjacent side on the triangle, it will be equivalent to zero, the term \displaystyle \|\vec{a}\|\cos{\theta} will be zero and the resulting multiplication with the magnitude of the vector \vec{b} will also be zero. Now you know that, when the dot product between two different vectors is zero, they are orthogonal to each other (they have an angle of 90 degrees), this is a very neat way to check the orthogonality of different vectors. It is also important to note that we are using 2D examples, but the most amazing fact about it is that we can also calculate angles and similarity between vectors in higher dimensional spaces, and that is why math let us see far than the obvious even when we can’t visualize or imagine what is the angle between two vectors with twelve dimensions for instance.

The Cosine Similarity

The cosine similarity between two vectors (or two documents on the Vector Space) is a measure that calculates the cosine of the angle between them. This metric is a measurement of orientation and not magnitude, it can be seen as a comparison between documents on a normalized space because we’re not taking into the consideration only the magnitude of each word count (tf-idf) of each document, but the angle between the documents. What we have to do to build the cosine similarity equation is to solve the equation of the dot product for the \cos{\theta}:

   \displaystyle   \vec{a} \cdot \vec{b} = \|\vec{a}\|\|\vec{b}\|\cos{\theta} \\   \cos{\theta} = \frac{\vec{a} \cdot \vec{b}}{\|\vec{a}\|\|\vec{b}\|}

And that is it, this is the cosine similarity formula. Cosine Similarity will generate a metric that says how related are two documents by looking at the angle instead of magnitude, like in the examples below:

The Cosine Similarity values for different documents, 1 (same direction), 0 (90 deg.), -1 (opposite directions).

Note that even if we had a vector pointing to a point far from another vector, they still could have an small angle and that is the central point on the use of Cosine Similarity, the measurement tends to ignore the higher term count on documents. Suppose we have a document with the word “sky” appearing 200 times and another document with the word “sky” appearing 50, the Euclidean distance between them will be higher but the angle will still be small because they are pointing to the same direction, which is what matters when we are comparing documents.

Now that we have a Vector Space Model of documents (like on the image below) modeled as vectors (with TF-IDF counts) and also have a formula to calculate the similarity between different documents in this space, let’s see now how we do it in practice using scikit-learn (sklearn).

Vector Space Model

Practice Using Scikit-learn (sklearn)

* In this tutorial I’m using the Python 2.7.5 and Scikit-learn 0.14.1.

The first thing we need to do is to define our set of example documents:

1
2
3
4
5
6
documents = (
"The sky is blue",
"The sun is bright",
"The sun in the sky is bright",
"We can see the shining sun, the bright sun"
)

And then we instantiate the Sklearn TF-IDF Vectorizer and transform our documents into the TF-IDF matrix:

1
2
3
4
5
from sklearn.feature_extraction.text import TfidfVectorizer
tfidf_vectorizer = TfidfVectorizer()
tfidf_matrix = tfidf_vectorizer.fit_transform(documents)
print tfidf_matrix.shape
(4, 11)

Now we have the TF-IDF matrix (tfidf_matrix) for each document (the number of rows of the matrix) with 11 tf-idf terms (the number of columns from the matrix), we can calculate the Cosine Similarity between the first document (“The sky is blue”) with each of the other documents of the set:

1
2
3
from sklearn.metrics.pairwise import cosine_similarity
cosine_similarity(tfidf_matrix[0:1], tfidf_matrix)
array([[ 1.        ,  0.36651513,  0.52305744,  0.13448867]])

The tfidf_matrix[0:1] is the Scipy operation to get the first row of the sparse matrix and the resulting array is the Cosine Similarity between the first document with all documents in the set. Note that the first value of the array is 1.0 because it is the Cosine Similarity between the first document with itself. Also note that due to the presence of similar words on the third document (“The sun in the sky is bright”), it achieved a better score.

If you want, you can also solve the Cosine Similarity for the angle between vectors:

   \cos{\theta} = \frac{\vec{a} \cdot \vec{b}}{\|\vec{a}\|\|\vec{b}\|}

We only need to isolate the angle (\theta) and move the \cos to the right hand of the equation:

   \theta = \arccos{\frac{\vec{a} \cdot \vec{b}}{\|\vec{a}\|\|\vec{b}\|}}

The \arccos is the same as the inverse of the cosine (\cos^-1).

 Lets for instance, check the angle between the first and third documents:
1
2
3
4
5
6
import math
# This was already calculated on the previous step, so we just use the value
cos_sim = 0.52305744
angle_in_radians = math.acos(cos_sim)
print math.degrees(angle_in_radians)
58.462437107432784
 And that angle of ~58.5 is the angle between the first and the third document of our document set.
That is it, I hope you liked this third tutorial !

Related Material

A video about Dot Product on The Khan Academy

Wikipedia: Dot Product

Wikipedia: Cosine Similarity

Scikit-learn (sklearn) – The de facto Machine Learning package for Python

Rastreamento em tempo real de avioes em Porto Alegre utilizando Raspberry Pi + Radio UHF (SDR RTL2832U)

* This post is in Portuguese.

SDR (Software-Defined Radio)

SDR é uma área de radiocomunicação baseada em uma ideia muito simples: implementar em software o que antes era implementado em hardware (ex: mod/demod, filtros, etc). O fato do SDR estar se tornando uma tendência hoje se dá principalmente pelo baixo custo de alguns receptores como por exemplo o RTL2832U (mais sobre ele depois) que hoje você pode encontrar facilmente por uns U$ 20.00, preço este muito barato se você levar em consideração o intervalo de cobertura das frequências de 22Mhz até 2200Mhz (dependendo do tuner). Neste intervalo dá pra se ter uma ampla cobertura de sinais de rádio AM/FM, rádio da polícia, TV, GSM, GPS, ADSB (este post é sobre o ADSB), comunicação marítma, LTE (ainda em definição no Brasil) e outras aplicações que estão inseridas neste intervalo de frequência.

Um sistema SDR é composto geralmente por: um receptor ligado ao computador ou qualquer outro dispositivo embarcado com um poder de processamento razoável (por meio de placa de captura de áudio, USB, etc) e um software que irá fazer o tratamento do sinal recebido. Neste post eu vou utilizar o Raspberry Pi como dispositivo embarcado para capturar e demodular os dados enviados pelos transponders presentes em aeronaves comerciais e domésticas.

ADS-B (Automatic dependent surveillance-broadcast)

Grande parte dos aviões modernos estão sendo equipados com um dispositivo que tem o objetivo de substituir os radares que existem hoje. Até o ano de 2020 todos os aviões que entrarem no espaço aéreo estadunidense deverão ter como item obrigatório um dispositivo compatível com ADS-B. O dispositivo ADS-B faz com que as aeronaves sejam visíveis aos radares em terra e também para outros aviões através do broadcast de mensagens com sua altura, velocidade, posição e muitas outras informações relevantes. A transmissão destas mensagens se dá através da frequência 1090Mhz, uma frequência dentro do intervalo de captura de maioria dos receptores SDR usando o chipset RTL2832U. A ideia deste post é usar o Raspberry Pi para receber este broadcast enviado diretamente pelo transponder dos aviões, demodular os frames e então utilizar um software para decodificar/interpretar os frames e plotar em um mapa a posição atual dos aviões em Porto Alegre / RS.

Realtek RTL2832U


Alguns receptores de TV digital USB utilizam o chipset da Realtek RTL2832U, como este meu acima. Em 2012 foi descoberto que este chipset permitia o envio de dados brutos do receptor para o host, permitindo assim seu uso para SDR. Existem alguns projetos com drivers para se comunicar com o dispositivo e receber estes dados (Linux e Windows), entre os quais o mais utilizado em Linux é o projeto rtl-sdr, que agrega além do driver alguns utilitários de linha de comando como por exemplo o “rtl_adsb” que utilizarei no Raspberry Pi para demodular o sinal ADS-B enviado pelas aeronaves; outro componente importante nos receptores USB é o tuner, que é responsável pelo ajuste da frequencia do rádio, no caso do meu dongle USB ele é o R820T que tem uma ótima sensibilidade mas tem um intervalo menor de cobertura do espectro quando comparado ao E4000 (da Elonics), veja a tabela no site do projeto para saber o intervalo de cobertura de cada tuner e de outros hardwares suportados pelo rtl-sdr.

Compilando o rtl-sdr no Raspberry Pi

Estou utilizando o Raspberry Pi como host do dongle USB porque ele é um aparelho barato e o consumo de energia é muito baixo, ou seja, você pode deixá-lo ligado capturando o sinal ADS-B por quanto tempo quiser sem se preocupar com um gasto maior que 5 watts no pior caso.

Baixando repositório

1
git clone git://git.osmocom.org/rtl-sdr.git

Dependências: libusb-1.0-0-dev, cmake, compilador (gcc)

 Criando arquivos para compilação e compilando:
1
2
3
4
5
6
7
cd rtl-sdr/
mkdir build
cd build
cmake ../
make
sudo make install
sudo ldconfig
E pronto, já estamos com o rtl-sdr instalado (driver a aplicativos extras).

Recebendo sinal ADS-B

Antes de mais nada, um pequeno teste para verificar se o rtl-sdr está encontrando corretamente o dongle USB:

1
2
3
4
5
6
7
8
9
10
11
pi@raspberrypi ~ $ rtl_test -t
Found 1 device(s):
 0: ezcap USB 2.0 DVB-T/DAB/FM dongle</pre>
<pre>Using device 0: ezcap USB 2.0 DVB-T/DAB/FM dongle
Found Rafael Micro R820T tuner
Supported gain values (29): 0.0 0.9 1.4 2.7 3.7 7.7
8.7 12.5 14.4 15.7 16.6 19.7 20.7 22.9 25.4 28.0
29.7 32.8 33.8 36.4 37.2 38.6 40.2 42.1 43.4 43.9
44.5 48.0 49.6
No E4000 tuner found, aborting.
pi@raspberrypi ~ $

Note que o meu dongle é um RTL2832U com o tuner R820T, ao executar o “rtl_test -t” é também exibida uma lista com os ganhos suportados pelo dongle.

Como a frequência utilizada pelo ADS-B é de 1090MHz, uma antena boa seria uma discone ou uma dipolo confeccionada com as dimensões corretas para esta frequência. Como ainda não tenho os plugs adequados tive que utiliar a antena que veio junto com o dongle, que mesmo sendo de baixa qualidade e não relacionada com a frequência do ADS-B ainda consegue receber os sinais (de fato, mesmo desconectando a antena eu consigo receber os frames das aeronaves), mesmo com o dongle em um ambiente fechado (minha casa fica há uns 6-10km do aeroporto de Porto Alegre / RS).

Para receber e demodular os frames eu utilizei o utilitário “rtl_adsb” (ele vem junto com o pacote do rtl-sdr). Ao ser executado, ele ajustará o tuner para a frequência do transponder das aeronaves em 1090MHz e logo após fará a demodulação dos frames para o formato hexadecimal:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
pi@raspberrypi ~ $ rtl_adsb
Found 1 device(s):
  0:  Realtek, RTL2838UHIDIR, SN: 00000013

Using device 0: ezcap USB 2.0 DVB-T/DAB/FM dongle
Found Rafael Micro R820T tuner
Tuner gain set to automatic.
Tuned to 1090000000 Hz.
Sampling at 2000000 Hz.
Exact sample rate is: 2000000.052982 Hz
*825566cf477b3124c64b17e74b15;
*e6c7d7fdb34c855db6972204ea14;
*d1e27bb95df2454ca547c87718a2;
*906bc5a59b5c5053226fb94f3460;
*a6f51dc76353efeeabfbfe6946cb;
*e78ed3aaf547fcd87e8e4f41cea3;
*d26a5ecacdb7051f2ebe18efc613;

Cada frame inicia sempre com um asterisco na frente e cada um destes frames foi enviado diretamente por um transponder de um avião ou é uma requisição terra->ar de alguma base em terra. O que precisamos agora é fazer a decodificação destes frames para poder extrair o tipo do frame, latitude, longitude, callsign, origem e destino do avião, etc. Não existem hoje muitas alternativas open-source para fazer o plot dos aviões em um mapa, o melhor que eu encontrei foi o Virtual Radar Server que é open-source e roda também em Linux, além de ter uma interface web bem amigável e plotar as aeronaves usando o Google Maps.
Para fazer com que o Virtual Radar Server conecte no Raspberry Pi, primeiramente precisamos fazer o streaming via TCP dos frames que o “rtl_adsb” está recebendo no Raspberry Pi, o que dá para resolver utilizando o “netcat” mesmo:

1
rtl_adsb | netcat -lp 8080

Neste caso ele irá escutar na porta 8080 e enviar os frames ADS-B para o cliente que conectar nele. Logo após é só especificar o IP/porta do Raspberry Pi no Virtual Radar Server e se tudo der certo e você conseguir receber algum frame de algum transponder, você verá uma tela parecida com esta mostrada abaixo com o rastreamento dos aviões em tempo real:

No screenshot você pode ver um avião da GOL (GOL1446) se dirigindo ao aeroporto (logo a frente onde diz São João) a uma velocidade de 250km/h e descendo a uma velocidade de 195 metros por minuto. Geralmente os aviões enviam o broadcast da mensagem com a posição a cada segundo, tudo vai depender da sua antena, receptor e condições de difusão. Com a minha antena pequena e no meio de prédios eu consegui capturar o broadcast de aviõe em até uns 80km de distância, espero melhorar isto com uma outra antena, só preciso achar o material agora =)

- Christian S. Perone