jeudi 31 juillet 2014

OLED





Un petit breakout avec un afficheur OLED 64x128 pixels piloté par un SSD1306 via i2c...  Chez Adafruit pour les puristes (ou pour quelques euros sur eBay pour les opportunistes).  Attention qu'il existe en SPI, en I2C,... et que c'est configuré en hardware (circuit imprimé et/ou pontages SMD) en plus du nombre de fils disponibles sur l'interface.

J'ai eu du mal à trouver un bout de code simple pour l'essayer sur mon Raspberry Pi.  Finalement, je me suis rabattu sur le code C++ proposé par JPB-HK sur la discussion 'I2C with Adafruit 128x32 I2C OLED' sur RaspberryPi.org (128x32i2c.tar.gz).

Comme je voulais un truc en 'C', le plus simple possible, j'ai taillé le code à la hache.  Ce qui donne finalement :

#include <sys/ioctl.h>
#include <unistd.h>
#include <fcntl.h>
#include <linux/i2c-dev.h>
#include <stdio.h>
#include <string.h>

typedef unsigned char uchar;
#define _GLCDFONT
#include "glcdfont.h"

#define I2C_PORT "/dev/i2c-1" // Using /dev/i2c-1 port
#define I2C_ADDR 0x3C         // Address for oled i2c, 2 possible, 0x3C and 0x3D

#define WIDTH    128          // 128 columns of 8bits deep
#define HEIGHT   4            // 4 rows, 8 bits high


uchar    oled_buf[HEIGHT * WIDTH];
int    i2c_fd;

void OLED_close()
    {
    if (i2c_fd)
        close(i2c_fd);
    }
void OLED_clear_buf()
    {
    memset(oled_buf, 0x00, (WIDTH * HEIGHT));
    }
void OLED_write_screen()
    {
    int     i, si = 0; /* si : screen index */
    uchar     buf[0x09];

    /* set screen pointer to (0,0) */
    OLED_i2c_write("\x00\x21\x00\x7f", 4);    // Set column address; start 0, end 127
    OLED_i2c_write("\x00\x22\x00\x03", 4);    // Set row address; start 0, end 3

    buf[0] = 0x40;           // Co = 0; D/C# = 1 for data
    while (si < (WIDTH * HEIGHT))
        {
        for (i = 0; i < 8; i++)
            buf[i+1] = oled_buf[si+i];
        OLED_i2c_write (buf, 9);
        si += 8;
        }
    }
void OLED_clear_screen()
    {
    OLED_clear_buf();
    OLED_write_screen();
    }
void OLED_write_text (char *str, int li, int co)
    {
    int    j;
    int    si;     /* screen index */

    li %= HEIGHT; co %= WIDTH;    // ensure in screen
    si = (li * WIDTH) + co;             // Set absolute start position in buffer

    while (*str)
        {
        for (j = 0; j < 5; j++)             // Loop through each 7bit segment of the character
            oled_buf[si++] = font[((*str)*5)+j];
        oled_buf[si++] = 0x00;                  // This puts a space between the characters on screen
        str++;
        }
    }
int OLED_init(void)
    {
    if ((i2c_fd = open(I2C_PORT, O_RDWR)) < 0) // Try to open /dev/i2c-x port
        {
        perror(I2C_PORT);
        return(-1);
        }

    if (ioctl(i2c_fd, I2C_SLAVE, I2C_ADDR) < 0) // Try to access the oled display
        {
        perror("i2c ioctl()");
        close(i2c_fd);
        i2c_fd = -1;
         return(-1);
        }
    /* commands are prefixed with 0x00 */
    OLED_i2c_write("\x00\xae", 2);        // Display off
    OLED_i2c_write("\x00\xd5\x80", 3);    // set display clock division
    OLED_i2c_write("\x00\xa8\x1f", 3);    // set multiplex
    OLED_i2c_write("\x00\xd3\x00", 3);    // set display offset
    OLED_i2c_write("\x00\x40", 2);        // set start line #0
    OLED_i2c_write("\x00\x8d\x14", 3);    // set charge pump
    OLED_i2c_write("\x00\x20\x00", 3);    // Memory mode
    OLED_i2c_write("\x00\xa1", 2);        // Segremap(0xA0 = reset, 0xA1 = 127 = SEG0)
    OLED_i2c_write("\x00\xc8", 2);        // Com scan dec (0xC0 = reset normal, 0xC8 = scan  from Com[n-1] - Com 0
    OLED_i2c_write("\x00\xda\x02", 3);    // Set com pins
    OLED_i2c_write("\x00\x81\x8f", 3);    // Set contrast
    OLED_i2c_write("\x00\xd9\xf1", 3);    // Set precharge
    OLED_i2c_write("\x00\xdb\x40", 3);    // Set Vcom select
    OLED_i2c_write("\x00\xa4", 2);        // Resume RAM content display
    OLED_i2c_write("\x00\xa6", 2);        // Normal display not inverted
    OLED_i2c_write("\x00\x21\x00\x7f", 4);    // Set column address; start 0, end 127
    OLED_i2c_write("\x00\x22\x00\x03", 4);    // Set row address; start 0, end 3
    return(OLED_i2c_write("\x00\xaf", 2));    // Display ON
    }
int OLED_i2c_write(uchar *ucp, int len)
    {
    if (write(i2c_fd, ucp, len) != len)
        {
        perror("OLED_i2c_write()");
        return(-1);
        }
    usleep(1);

    return(len);
    }

main()
    {
    OLED_init();
    OLED_clear_screen();
    OLED_write_text("Oufti!", 10, 2);
    OLED_write_screen();
    }

Cela pourrait être un peu plus propre mais pour le moment, c'est bon comme ça. Cela permet de constater que l'afficheur fonctionne sans se casser la tête. (Par exemple, OLED_write_text() devrait proposer x, y avant le texte et 'x' est en pixels et 'y' en ligne de texte...)

Le fichier se compile en faisant simplement 'gcc -o oled oled.c'.  Le fichier 'glcdfont.h', qui contient les caractères en bitmap, est quasi l'original (juste remplacé 'unsigned char' par 'uchar').  Le code de JPB-HK est dérivé de ce qu'il avait trouvé sur AdaFruit.com et était abondamment commenté.  En plus d'avoir écrémé le C++, j'ai fini par également enlever les commentaires qui, à mon sens alourdissaient le fichier inutilement.  En fait, il n'y a que trois choses importantes dans le code : l'initialisation (OLED_init()), l'écriture dans le frame-buffer local (OLED_write_text(string, x, y)) et l'envoi du frame-buffer vers l'afficheur (OLED_write_screen()).

L'initialisation est complexe.  Outre qu'il faille correctement configurer le RPi avant de pouvoir accéder aux périphériques sur le bus /dev/i2c-1 avec les i2c-tools (par exemple, 'sudo i2cdetect -y 1' nous apprend que le display est bien à l'adresse 0x3C), une petite vingtaine de commandes semblent nécessaires pour initialiser le chip.  Je ne suis pas sûr qu'on y arrive en lisant seulement la datasheet du SSD1306.  Une fois acquis ce 'Hello World!', elle est probablement indispensable pour explorer plus avant les fonctionnalités de ce composant qui semble offrir de nombreuses possibilités.

'sudo ./oled' affiche 'Oufti!' sur la troisième ligne de l'afficheur.

lundi 14 octobre 2013

ADS-B et RTL SDR

En cours de rédaction...
 Un 'squitter' de 56 bits fort et clair, bien isolé, au milieu du bruit de fond.

Petite exploration dans la réception des messages qu'émettent les avions sur 1090 MHz...  Grâce à certaines clés de réception de TV et radio numérique, il est possible de capter les messages ADS-B émis par les avions.  Ces messages de 56 ou 112 bits contiennent des informations sur leur identification, leur position, leur vitesse,...  La captation de ces messages et leur mise en réseau permettent à des sites comme FlightRadar24.com d'afficher la position de la plupart des avions en vol dans les régions couvertes.

Je ne vais m'intéresser ici qu'à la réception des 'squitters' et non au décodage des messages qui semble assez complexe, étant le fruit d'une longue évolution (on le trouve déjà, par exemple, dans ce document.pdf de 1974 (p.14)).  L'encodage contient des archaïsmes comme le code Gray pour l'encodage de l'altitude, des latitude/longitude,...

Intéressons-nous donc uniquement à la réception des messages...  On est à la limite du possible avec ce type de matériel.  Pas tellement pour une question de fréquence (certaines clés peuvent monter jusqu'à 2000 MHz) mais pour une question d'échantillonnage : l'ADS-B émis par les avions à 1090 MHz a un débit binaire de 1 Mbps avec une modulation en position d'impulsions.
Ces impulsions font 0.5 µS et ces clés USB ne nous autorisent pas plus de 2 millions d'échantillons par seconde.  Ce qui fait que, soit le message est bien positionné dans le temps, on l'échantillonne aux moments propices et on a un espoir de pouvoir le décoder; soit on échantillonne au mauvais moment et on le rate complètement avec aucun espoir de le récupérer.  Le théorème de Nyquist-Shanon nous apprends, en effet, que pour bien faire, il aurait fallu échantillonner à 4 millions d'échantillons par seconde pour bien voir les impulsions de 0.5 µS.  Mais, bon...  Comme les avions sont très bavards (du genre un message par avion et pas seconde), on en capte quand même beaucoup.

Les clés comportent principalement deux composants électroniques : un 'tuner' que l'on règle sur 1090 MHz et qui descend le signal à une fréquence intermédiaire que le fameux Realtek RTL2832U, démodulateur COFDM, qui échantillonne le signal et envoie les I/Q sur le port USB.  Le driver Osmocom (git://git.osmocom.org/rtl-sdr) permet de piloter l'ensemble.  Il existe un module pour GNU Radio (gr-air-modes) mais je préfère partir de dump1090 qui est plus simple et, en même temps, plus avancé.
Diagramme de constellation du 'squitter' *5D4CA27B919917; tel qu'il sort de la clé.  La construction de la table des magnitudes (maglut) par ModesInit() dans dump1090.c nous apprend qu'il s'agit d'octets non signés (0..255 représentant -1..1).  Pour la petite histoire, ce message de type DF-11 provient d'un B737-800 de Ryanair (code AOCI 4CA27B correspondant au 'callsign' RYR39TY) qui volait à 38 000 pieds dans la région (comme nous l'apprend un message de type DF-4 avec un CRC égal à 4ca27b (!) un peu plus tard...).  Il avait une 'bonne qualité de signal' (50 points entre les '*' et les '_' du préambule).

Le but du jeu est ici d'extraire le code minimum qui permet de capter les messages ADS-B.

Dans un premier temps, on s'assure que l'on peut initialiser la clé correctement.  Le court programme suivant initialise la clé et envoie tout ce qu'on y lit dans la 'sortie standard'.  En la redirigeant dans un fichier, on obtient un fichier binaire brut qui peut servir d'entrée à 'dump1090 --ifile monfichier.bin'.  Si celui-ci parvient à récupérer des messages, on est bon.  Les routines appelées par init() on été retrouvée dans le code de dump1090 avec l'aide de l'API décrit dans rtl-sdr.h.

#include <stdlib.h>
#include <stdio.h>
#include <errno.h>
#include <limits.h>
#include <signal.h>
#include <math.h>
#include "rtl-sdr.h"
  #include <string.h>
/*
** fast magnitude calculation mag = maglut[IQ]
*/
int error(char *str)
 {
 perror(str);
 return(-1);
 }
rtlsdr_dev_t *devp;
int running = 1;
int onint(int signum)
 {
 running = 0;
 }
int init(int sps, int freq)
 {
 int numgains;
 int gains[100];

 signal(SIGINT, onint);
 if (rtlsdr_open(&devp, 0) < 0)
  return(error("Can't open device!"));

 rtlsdr_set_tuner_gain_mode(devp, 1);  /* manual */

 numgains = rtlsdr_get_tuner_gains(devp, gains);
 rtlsdr_set_tuner_gain(devp, gains[numgains-1]); /* set max gain */
 fprintf(stderr, "Tuner gain = %d\n", rtlsdr_get_tuner_gain(devp));
 rtlsdr_set_freq_correction(devp, 0);  /* why not */
 rtlsdr_set_center_freq(devp, freq);
 rtlsdr_set_sample_rate(devp, sps);
 rtlsdr_reset_buffer(devp);

 return(0);
 }
#define BUFSIZE (32*1024)
int main(int argc, char *argv[])
 {
 ushort buf[BUFSIZE];
 int n_read=0;
 int len;
 int sps, freq;

  if (argc < 3)
  {
  fprintf(stderr, "rtl_grab: usage is 'rtl_grab sps freq'\n");
  exit(-1);
  }
 sps = atoi(argv[1]);
 freq = atoi(argv[2]);
 if (sps > 2024000) /* 2Msps */
  {
  fprintf(stderr, "sps too high\nrtl_grab: usage is 'rtl_grab sps freq'\n");
  exit(-1);
  }
 if (freq < 27000000 || freq > 2000000000)
  {
  fprintf(stderr, "freq too low or too hig\nrtl_grab: usage is 'rtl_grab sps freq'\n");
  exit(-1);
  }

 fprintf(stderr, "rtl_grab : %d samples/sec -- freq = %d Hz (-C to stop)\n", sps, freq);
 signal(SIGINT, onint);
 if (isatty(1))
  {
  fprintf(stderr, "Output should be redirected.\n");
  exit(-1);
  }
 if (init(sps, freq) < 0)
  exit(-1);

 while (running  && (len = rtlsdr_read_sync(devp, buf, sizeof(buf), &n_read)) >= 0)
  {
  if (write(1, buf, n_read) != n_read)
   {
   perror("Can't write.\n");
   break;
   }
  }
 rtlsdr_close(devp);
 return(0);
 }


Si je le laisse tourner pendant environ 8 secondes, ça me donne un fichier de 32 MB (2 octets * 2 millions d'échantillons par seconde * 8 secondes) que je peux utiliser avec '$ ./dump1090 --stats --ifile iq-130926.120740.bin', ce qui me donne :


16240 valid Mode-S preambles

26 DF-?? fields corrected for length
10 DF-?? fields corrected for type
219 demodulated with 0 errors
2 demodulated with 1 error
0 demodulated with 2 errors
0 demodulated with > 2 errors
221 with good crc
0 with bad crc
0 errors corrected
   0 with 1 bit error
   0 with 2 bit errors
221 total usable messages


Donc, jusque là, c'est bon, la routine init() fonctionne et rtlsdr_read_sync() me donne bien des données binaires brutes sur lesquelles je peux 'travailler'.

Le pas suivant, c'est de détecter le 'préambule', un 'pattern' '*_*____*_*' suivi de '___' qui précède un message ('*' = grande amplitude; '_' faible amplitude).  L'amplitude, c'est la racine carrée de la somme de I² + Q² avec les I/Q lus par  rtlsdr_read_sync().  Dans le code de dump1090, cela donne quelque chose comme :



for (p = magvec, i = 0; i < cnt;)
 { /*    *_*____*_*     0_2____7_9  */
 if ( p[0] < p[1] || p[2] < p[1]
  || p[2] < p[3]
  || p[3] > p[0] || p[4] > p[0] || p[5] > p[0] || p[6] > p[0]
  || p[7] < p[8] || p[8] > p[9]
  || p[9] < p[6])
  { /* not a valid preamble */
  p++; i++;
  continue;
  }
 h=(p[0]+p[2]+p[7]+p[9])/4;
 if (p[4] >= h || p[5] >= h)
  { /* not a valid preamble */
  p++; i++;
  continue;
  }
 if (p[11] >= h || p[12] >= h || p[13] >= h || p[14] >= h)
  { /* not a valid preamble */
  p++; i++;
  continue;
  }
 sstrength = p[0]+p[2]+p[7]+p[9] - (p[1]+p[3]+p[6]+p[8]);
 good_preamble_t0 += 1;


Mais, on remarquera que toutes ces précautions sont loin d'être suffisantes pour s'assurer qu'il s'agit bien d'un message ADS-B valide : en 8 secondes, dump1090 trouve 16240 préambules alors qu'il ne décode correctement que 221 messages valides (1/73!).  En fait, même dans un signal complètement aléatoire, on détecte une masse de préambules.

On va maintenant s'assurer que l'on reçoit 56 (ou 112) impulsions codées PPM (haut-bas ou bas-haut mais pas bas-bas ni haut-haut).  Les échantillons qui se suivent doivent être différents; si le premier est plus grand que le second, c'est un '1'; si le second est plus grand que le premier, on a un '0'; s'ils sont égaux, on a un problème...  À noter que l'égalité dépend probablement surtout du gain que l'on applique au signal : trop faible, on va avoir des zéro; trop fort des 'maximums' (sqrt(FF²+FF²).  C'est aussi possible pour les valeurs intermédiaires mais c'est alors la faute à 'pas de chance' (qui, sur 8 bits, n'est pas vraiment négligeable).  Dump1090 se lance dans certains paris pour quand même récupérer certains messages mais cela me paraît hasardeux (il n'en reste pas moins que son truc fonctionne et le mien, pas encore...)

On pourrait croire qu'il suffit maintenant de prendre 56 paires de magnitudes et de s'assurer que le CRC est nul mais il n'en est rien.  Il faut tenir compte du type de message DFxx, non seulement pour voir s'il ne s'agit pas d'un message de 112 bits, mais aussi parce que la valeur du CRC varie avec le type de message (!) (comme on peut le voir dans la routine decodeModesMessage() de mode_s.c; ce serait évidemment mieux de trouver les spécifications des messages...).


Le 'squitter' *5D4CA27B919917; en magnitude

for (byte = 0, err = 0; byte < 7; byte++)
 {
 x = 0;
 for (bit = 0; bit < 8; bit++)
  {
  x <<= 1;
  a = *mp++; b = *mp++;
  if (a > b)
   x |= 1;
  }
 msg[byte] = x;
 }
...À suivre...

mardi 10 janvier 2012

Rendez-vous avec la Station spatiale internationale

Des sites comme Heavens-Above donne les éphémérides d'observation en fonction du lieu. Parfois ISS est visible le soir, parfois le matin et parfois pas du tout. La visibilité suit un cycle immuable. Elle n'est observable confortablement que lorsque le Soleil est déjà couché/pas encore levé au lieu d'observation mais encore/déjà en plein soleil environ 400 km au dessus de notre tête à la manière d'un coucher/lever de soleil sur une montagne voisine. Comme son altitude varie au fil du temps (l'atmosphère la freine et elle est régulièrement 'remontée'), sa période orbitale varie également. Il est donc difficile de prévoir le moment exact de son passage longtemps à l'avance. Les logiciels qui font cela se basent sur des mesures récentes qui sont publiées régulièrement sous forme de TLE.

Mais, peut-on, à la manière des anciens qui prédisaient les éclipses, se faire une idée des périodes d'observations? En gros, alors que je prévois un voyage d'une semaine en Papouasie dans six mois, puis-je savoir aujourd'hui si j'aurai l'occasion de montrer un passage d'ISS lors d'un souper à mes amis papous?

Il semble que oui.

Il y a deux paramètres importants dans l'orbite d'ISS qui déterminent sa visibilité en soirée (ou avant l'aube) en fonction de l'époque de l'année : l'inclinaison du plan de l'orbite et la coordonnée du nœud ascendant. L'inclinaison du plan de l'orbite est remarquablement constant, il semble maintenu au centième de degré près. Le nœud ascendant, quant à lui, subit un mouvement de précession qui semble assez stable pour pouvoir effectuer des prédictions valables.

Les données historiques collectées par Jonathan McDowell permettent de s'en rendre compte.

Inclinaison du plan de l'orbite


En ce qui concerne l'inclinaison, une simple commande Unix permet d'extraire l'inclinaison de la liste de TLE :
awk  '/^2/{print $3;}' < ISS-S25544-historical.tle.txt

qui imprime la troisième colonne des lignes commençant par '2' de ISS-S25544-historical.tle.txt. Ce qui, graphiquement, donne :


On peut remarquer deux 'anomalies'. Au début du graphe, un 'réglage/modification de l'orbite' (?) effectué aux environs du 4 juin 1999 et un 'glitch' (?) vers le 15 mars 2007. À part cela, c'est remarquablement stable, entre 51.63 et 51.65° (un demi pour mille). En zoomant sur les données, on peut observer des corrections mineures :


Précession nodale


Le plan de l'orbite est animé d'un mouvement de précession. Le quatrième paramètre de la seconde ligne est l'ascension droite du nœud ascendant. Il convient ici d'être plus précis et de regarder comment ce paramètre évolue en fonction du temps. Pour simplifier, intéressons-nous à l'année 1999. Si on fait un graphique avec le temps en abscisse (on le trouve en ligne un des TLE) et l'ascension droite en ordonnée, on obtient le graphe suivant :

Ce qui semble encourageant.
Si je m'arrange pour éviter le 'modulo 360', j'obtiens quelque chose qui ressemble (de loin) à une belle droite.


Reste à vérifier que cela ne s'éloigne pas trop de la droite idéale.
(Je comptais utiliser GNU R mais je ne sais pas encore m'en servir)

Autre option : la théorie...

En fait, Michel Capderou, dans Satellites : orbites et missions donne une formule pour la vitesse de précession en fonction de l'inclinaison et du demi grand-axe (page 134).


où le résultat est en radians par seconde; J2 est le 'facteur d'ellipticité potentiel de la Terre' (0.0010826362); μ est la 'constante géocentrique de gravitation' (3.986004415E+14 m³/s²); R est le rayon équatorial de la Terre (6378136.3 m); a, le demi grand-axe (R+altitude; m); i l'inclinaison de l'orbite.

En fonction de l'altitude, pour ISS, inclinée à 51.64°, cela donne :

On peut voir qu'une variation de 40 km (ce qui plutôt beaucoup), on a une variation de vitesse d'un dixième de degré par jour. Ce qui est relativement peu : en 6 mois, cela ne fait 'que' 18°, ce qui correspond à un peu plus d'une heure de rotation terrestre. Comme les passages de ISS sont visibles pendant plusieurs heures après le coucher du Soleil, cela semble acceptable.

Sur l'image suivante, on peut observer l'évolution de l'orbite en rouge et du terminateur terrestre en bleu entre le 16 et le 28 mars.


Cliquez sur l'image pour la voir en entier (<esc> pour revenir)

En abscisse, l'heure solaire avec midi au centre et en ordonnée la latitude, le nord au dessus, l'équateur au centre. On voit bien la précession de 5° : il faut trois jours pour que l'orbite se décale d'une heure et le changement de saison (le Soleil disparaît au Pôle sud et apparaît au Pôle nord.

La page qui dessine cela est ici. Elle utilise les canvas de HTML5 et nécessite donc un navigateur récent. La page permet (en principe) de faire le graphique pour n'importe quel satellite à n'importe quelle date. Il faut fournir un TLE pas trop éloigné de la date en question (sinon le résultat vire au n'importe quoi).

Le graphique donne l'heure (locale) solaire de passage à une latitude donnée. En regardant du côté des intersections entre l'orbite et le terminateur, on peut évaluer les moments et les lieux où le satellite est observable (éclairé par le Soleil au dessus d'une région plongée dans la nuit). Dans ce cas-ci, on voit qu'il est observable le matin vers 50° de latitude nord et le soir vers 50° de latitude sud. ISS est observable dans un rayon d'environ 1000 kilomètres, ce qui se traduit par +/- 10° de latitude en vertical mais elle reste visible longtemps après le passage du terminateur (combien?; et plus en été qu'en hiver(?)).

Et donc, pour répondre à la question initiale, oui, s'il n'y a pas de changement de politique en ce qui concerne l'altitude de la station spatiale internationale, on peut prédire longtemps à l'avance si elle sera visible le matin ou le soir. Le graphe suivant montre les prévisions de la position du nœud ascendant début mars 2012 pour EGP et pour ISS à partir des TLE compris entre juin 2011 et début avril 2012 ( AD(9avr)= AD(t) + dn*(9avr-t) ).



On voit qu'une fois faite une modification importante de l'orbite d'ISS mi-juin 2011, son altitude est devenue suffisamment stable pour prévoir dès juillet 2011 que le nœud ascendant se trouverait vers 107° le 9 avril 2012 (EGP, moins sujet au freinage atmosphérique, donne 3° d'erreur à un an).

L'évolution de la vitesse de précession nodale de l'orbite de ISS montre qu'elle n'est pas négligeable (elle est liée à la variation de son altitude) mais cela permet quand même des prédictions raisonnables. Quand bien même utiliserait-on une précession journalière différente de 1/10 de degré, au bout de 150 jours, cela fait 15°, c'est-à-dire une heure de décalage entre la prévision du passage sous le plan de l'orbite et le passage effectif.



En se basant sur la vitesse de précession début juin 2011, on fait une erreur de 30°, 300 jours plus tard (9 avril 2012) parce que, vers le 15 juin 2011, ISS a été remontée d'une cinquantaine de kilomètres et qu'elle est ensuite restée plus ou moins à la même altitude.

mardi 8 novembre 2011

Tablette Android Xiron ATF3416


Achetée chez Blokker fin octobre 2011 en Belgique à 99.95€ (elle avait été soldée 90.00€ peu avant). Après avoir tourné en rond (entre un compte Google sans mobile et un mobile sans compte Google) sans succès pour essayer d'accéder à l'AndroidMarket, j'ai téléchargé Term.apk sur code.Google.com... Et cela permet d'accéder à un shell avec les commandes Android/Linux/Unix standard. 'Top' m'apprend qu'il y a 3*64 MB de RAM et 'cat /proc/cpuinfo' que le processeur est un ARM926EJ-S rev. 5 (v51) donnant 797.97 BogoMIPS. Et 'uname -a' que le kernel Linux est un 2.6.32.9 compilé le 5 août 2011 pour un armv5yej1.

mercredi 4 mai 2011

Code Parsons


'Profil mélodique' du dernier mouvement de Herz und Mund und Tat und Leben (J.S. Bach, Jésus que ma joie demeure).

Qui n'a jamais eu un air en tête sur lequel il est incapable de mettre un nom? Pour le non musicien, c'est d'autant plus frustrant qu'il est incapable de noter la mélodie, qu'il risque de l'oublier et ne peut la communiquer à d'autres qu'en la chantonnant, la sifflant, etc... Dans des circonstances qui ne sont jamais favorables. Le résultat est souvent désastreux...

Le code Parsons propose un remède à ce problème. Inventé par Denys Parsons en 1975, c'est une notation musicale accessible au non-initié. ...Enfin presque.

Le code correspondant à une mélodie est constitué d'une suite de quatre symboles :
- '*', la première note de la mélodie
- 'U', une note plus 'haute' (aigüe) que la précédente (Up)
- 'D', une note plus 'basse' que la précédente (Down)
- 'R', une note de la même hauteur que la précédente (Repeat)
Ceci, indépendamment du rythme, de l'écart entre les notes, etc... Cela ne décrit donc pas la mélodie avec précision, ce code est une espèce de 'signature'. Il est impossible de jouer la mélodie en décodant le code Parsons. À la limite, on peut essayer de deviner. Par contre, on peut vérifier assez facilement que le code Parsons correspond bien à la mélodie.

Denys Parsons avait compilé plus de 10 000 mélodies dans un ouvrage appelé Directory of Tunes and Musical Themes en 1975 (et ré-édité/amendé en 2010 chez Piatkus Books sous le nom The Directory of Classical Themes. Il s'agit d'une espèce de dictionnaire de mélodies triées par code.

À partir du code, ce bouquin ou le site Musipedia me permet de mettre un nom sur la mélodie. Donc, c'est génial. Sauf que, pour un non-initié de mon acabit, il est difficile de décider laquelle de deux notes qui se suivent est la plus grave ou même s'il s'agit de la même note. Sauf dans les cas les plus évidents. Et, encore... Par exemple, je suis incapable de trouver *UUDRUUDUUUDUU à partir de Frère Jacques.

C'est là que l'informatique vient à mon secours. Si j'enregistre mon sifflement, sous Linux, sndfile-spectrogram (de sndfile-tools) me permet de visualiser la hauteur des notes :


On voit que, sifflé, Jésus que ma joie demeure donne un beau sonagramme bien net et relativement facile à décoder. On trouve assez facilement uuudrudrududdduudud que la recherche par contour mélodique sur http://musipedia.org trouve facilement.

On voit que le sifflement se passe du côté de 1000 Hz (une octave plus haut que le La 440 Hz); il était donc inutile d'échantillonner à 22.05 Khz; échantillonné à 4 KHz, cela donne :



C'est bien, mais ce serait encore mieux d'écrire un petit programme qui génère le code à partir du sifflement (comme le fait le site si on active le micro en FLASH; mauvaise idée...). En fait, il 'suffit' de calculer la transformée de Fourier du signal audio capté par le micro. Reste à trouver le code pour accéder au micro et le code pour faire une FFT.

* Parrot.c, un petit programme de démonstration d'accès à l'audio sous Linux
* Discrete Fourier Transform à pied un tutorial sur les FFT comprenant un code simple.

Accès aux données audio sous Linux


Pour pouvoir analyser un son, il faut accéder à l'échantillonnage du signal. /dev/dsp nous donne accès aux données brutes du PCM. Parrot.c montre comment initialiser le périphérique pour y lire 8000 échantillons de 8 bits par seconde. On en dérive un petit programme (minimum) très simple qui, suivant qu'il s'appelle 'r' (record) ou 'p' (play), copie les octets échantillonnés de /dev/dsp dans la sortie standard (utilisation : $ r > mon_enregistrement.raw) ou 'joue' les octets échantillonnés de l'entrée standard dans /dev/dsp (utilisation $ p < mon_enregistrement.raw). Il n'en faut pas plus.
#include &ltstdlib.h>
#include &ltstdio.h>
#include <fcntl.h>
#include <sys/ioctl.h>
#include <linux/soundcard.h>

#define SAMPLE_RATE 8000 /* the sampling rate */
#define SAMPLE_SIZE 8 /* sample size: 8 or 16 bits */
#define CHANNELS 1 /* 1 = mono 2 = stereo */

typedef unsigned char uchar;
uchar audio_buf[1024];

int main(int argc, char *argv[])
{
int fd, n;
int arg; /* argument for ioctl calls */

if ((fd = open("/dev/dsp", O_RDWR)) < 0)
{
perror("open of /dev/dsp failed");
exit(1);
}
arg = SAMPLE_SIZE; /* sample size */
if (ioctl(fd, SOUND_PCM_WRITE_BITS, &arg) < 0)
perror("SOUND_PCM_WRITE_BITS ioctl failed");
if (arg != SAMPLE_SIZE)
perror("unable to set sample size");

arg = CHANNELS; /* mono or stereo */
if (ioctl(fd, SOUND_PCM_WRITE_CHANNELS, &arg) < 0)
perror("SOUND_PCM_WRITE_CHANNELS ioctl failed");
if (arg != CHANNELS)
perror("unable to set number of channels");

arg = SAMPLE_RATE; /* sampling rate */
if (ioctl(fd, SOUND_PCM_WRITE_RATE, &arg) < 0)
perror("SOUND_PCM_WRITE_RATE ioctl failed");

if (argv[0][0] == 'p')
{ /* play stdin into /dev/dsp */
while ((n = read(0, audio_buf, sizeof(audio_buf))) > 0)
write(fd, audio_buf, n);
}
else if (argv[0][0] == 'r')
{ /* record /dev/dsp into stdout */
while ((n = read(fd, audio_buf, sizeof(audio_buf))) > 0)
write(1, audio_buf, n);
}

exit(0);
}

On pourrait aussi se servir des outils standards, l'enregistreur de son et le 'couteau suisse' Sox (SOund eXchange) pour convertir les différents formats audio de/vers un 'raw' 8 bits à 8 KHz (ce qui est suffisant pour l'usage qu'on en fait). L'avantage du raw 8 bits, c'est que l'on n'a pas besoin de se préoccuper des en-tête, d'un autre codage,... Les outils font ça très bien; moi, tout ce que je veux faire, c'est analyser le signal audio pour voir si je peux détecter la différence de hauteurs (fréquences) entre deux notes sifflées dans le micro. (Si je veux communiquer les fichiers audio à quelqu'un qui n'a que des outils standard, j'utilise LAME pour produire un .mp3)

Je peux m'assurer que tout cela fonctionne en faisant un petit test (j'enregistre 'Frère Jacques', je transforme mon fichier .raw en un .wav pour pouvoir en faire un sonagramme avec sndfile-spectrogram puis, je 'joue' le morceau tel que je l'ai enregistré; l'enregistrement est interrompu par ^C (-C))):
$ ./r > FJ.raw
$ sox -v 2.0 -r 8k -b 8 -e unsigned -c 1 FJ.raw FJ.wav
$ sndfile-spectrogram FJ.wav 640 480 FJ.png
$ ./p < FJ.raw

Il faut noter que la gestion du son sur Linux n'est pas des plus simple. ALSA a remplacé OSS et le contrôle se fait via divers programmes comme alsamixer, alsactl, gnome-audio-contrl,... (et, par exemple, il est difficile de savoir où se trouve la pré-amplification du micro dans une Ubuntu 8.04)

Traitement du signal


'On peut montrer que tout signal périodique est décomposable en une somme de sinusoïdes' ...ou quelque chose comme ça. Quelles sinusoïdes?, avec quelles intensités?, comment les calculer?

FFT à pied explique que
* les fréquences sont [1..N/2] * fréquence_échantillonnage/N
* l'intensité est la somme des produits de la mesure par le sinus/cosinus à 'cet endroit'

C'est assez simple, mais il existe une méthode plus efficace qui élimine des calculs redondants : la FFT. C'est cette méthode qui est utilisée dans la pratique. On passe d'un temps calcul en N² en Nlog(N); pour N=1024, on passe de 1M à 10K (on boucle 128 fois moins (mais la boucle est un peu plus complexe)).

J'aurais sans doute pu utiliser KissFFT ou un autre librairie plus ou moins généraliste et optimisée. Mais, pour faire simple, j'utilise le code court et auto-suffisant de Discrete Fourier Transform à pied :
//
// Listing 1.4: The Discrete Fast Fourier Transform (FFT):
//

#define M_PI 3.14159265358979323846

void smbFft(float *fftBuffer, long fftFrameSize, long sign)
/*

FFT routine, (C)1996 S.M.Bernsee. Sign = -1 is FFT, 1 is iFFT (inverse)

Fills fftBuffer[0...2*fftFrameSize-1] with the Fourier transform of the time
domain data in fftBuffer[0...2*fftFrameSize-1]. The FFT array takes and returns
the cosine and sine parts in an interleaved manner, ie.
fftBuffer[0] = cosPart[0], fftBuffer[1] = sinPart[0], asf. fftFrameSize must
be a power of 2. It expects a complex input signal (see footnote 2), ie. when
working with 'common' audio signals our input signal has to be passed as
{in[0],0.,in[1],0.,in[2],0.,...} asf. In that case, the transform of the
frequencies of interest is in fftBuffer[0...fftFrameSize].

*/
{
float wr, wi, arg, *p1, *p2, temp;
float tr, ti, ur, ui, *p1r, *p1i, *p2r, *p2i;
long i, bitm, j, le, le2, k, logN;
logN = (long)(log(fftFrameSize)/log(2.)+.5);

for (i = 2; i < 2*fftFrameSize-2; i += 2) {
for (bitm = 2, j = 0; bitm < 2*fftFrameSize; bitm <<= 1) {
if (i & bitm) j++;
j <<= 1;
}
if (i < j) {
p1 = fftBuffer+i; p2 = fftBuffer+j;
temp = *p1; *(p1++) = *p2;
*(p2++) = temp; temp = *p1;
*p1 = *p2; *p2 = temp;
}
}

for (k = 0, le = 2; k < logN; k++) {
le <<= 1;
le2 = le>>1;
ur = 1.0;
ui = 0.0;
arg = M_PI / (le2>>1);
wr = cos(arg);
wi = sign*sin(arg);
for (j = 0; j < le2; j += 2) {
p1r = fftBuffer+j; p1i = p1r+1;
p2r = p1r+le2; p2i = p2r+1;
for (i = j; i < 2*fftFrameSize; i += le) {
tr = *p2r * ur - *p2i * ui;
ti = *p2r * ui + *p2i * ur;
*p2r = *p1r - tr; *p2i = *p1i - ti;
*p1r += tr; *p1i += ti;
p1r += le; p1i += le;
p2r += le; p2i += le;
}
tr = ur*wr - ui*wi;
ui = ur*wi + ui*wr;
ur = tr;
}
}
}

Ce code, très simple et évident (en ce sens qu'il ne s'agit 'que' de quelques boucles imbriquées et d'affectations), ne semble pas optimisable de manière flagrante. Ce n'est pas très long, c'est vraisemblablement le cas d'école de la FFT, l'utilisation est évidente : 3 paramètres et cela calcule la FFT .

Premier test


Comme premier test, on peut traiter 1024 échantillons d'un sifflement (ou deux) et calculer la magnitude du signal aux fréquences calculées par la FFT.
unsigned char audio_buf[1024];
float fft_buf[2*1024]; /* alternance réel+imaginaire */

void bufcpy(float *fp, uchar *ucp, int len)
{ /* mise en place des mesures audio_buf dans le fft_buffer */
while (--len >= 0)
{
*fp++ = (float)*ucp++;
*fp++ = 0.0; /* i-part */
}
}
void print_mag(float *fp, int len)
{
float r_f, i_f;

while (--len >= 0)
{
r_f = *fp++;
i_f = *fp++;
fprintf(stdout, "%f\n", sqrt(r_f*r_f+i_f*i_f));
}
}
main()
{ /* ~pseudo-code */
...
read(fd, audio_buf, sizeof(audio_buf));
bufcpy(fft_buf, audio_buf, 1024);
smbFft(fft_buf, 1024, -1);
print_mag(fft_buf, 1024);
}

On obtient quelque chose comme :


On repère bien deux pics à certaines fréquences (mais l'abscisse est en 'indices de buffer'; pour connaître la fréquence correspondante, il faut multiplier par 8000/512). L'indice 0 a été supprimé (il correspond au niveau de base (?)) et on remarquera que l'information utile se trouve entre 0 et 512, le reste étant un calcul miroir (c'est produit par la méthode de calcul (FFT); de toute façon on ne calcule pas de sinus de fréquence supérieure à la moitié de la fréquence d'échantillonnage (Théorème d'échantillonnage de Nyquist-Shannon)

Il faut noter que cette échelle est linéaire alors que notre perception musicale ne l'est pas; l'oreille est sensible aux rapports entre les notes : une octave correspond à un doublement de fréquence. L'octave la plus grave du piano va de 32.70 à 65.41 Hz (environ 33 Hz pour 12 notes (les 7 blanches et 5 noires du piano - la gamme tempérée); la plus aigüe de 2093,00 Hz à 4186,01 Hz (2093 Hz pour 12 notes) (ref.). On a donc intérêt à siffler 'aigu', il y aura plus de Herz entre les notes, elles seront plus faciles à individualiser.

Contour mélodique


En effectuant une boucle et en n'imprimant que l'indice de la fréquence maximum, on obtient une espèce de sonagramme simplifié, un contour mélodique.
int get_max(float *fp, int len)
{
int i, imax, n;
float max, mean, mag2;

mean = max = 0.0; n = 0;
for (i = 1; i < len; i++)
{
mag2 = fp[2*i] * fp[2*i] + fp[2*i+1] * fp[2*i+1];
if (mag2 > max)
{
imax = i;
max = mag2;
}
mean += mag2; n += 1;
}
if (max < (3.0 * mean/n))
imax = 0; /* no signal */
return(imax);
}
void print_max(float *fp, int len)
{
fprintf(stdout, "%d\n", get_max(fp, len));
}
main()
{
... ~ pseudo-code
for (;;)
{
read(fd, audio_buf, sizeof(audio_buf));
bufcpy(fft_buf, audio_buf, 1024);
smbFft(fft_buf, 1024, -1);
print_max(fft_buf, 512); /**/
}
}

Par exemple, pour Frère Jacques, on obtient quelque chose comme


On remarque que la cinquième note n'est pas exactement au niveau de la quatrième (alors qu'elle devrait...). Il faudra donc être un peu tolérant (ou siffler plus juste).

Cela s'annonce donc bien. On peut noter que l'on peut jouer sur deux paramètres : la fréquence d'échantillonnage et taille du buffer. Dans la configuration actuelle avec 8000 échantillons par seconde et un buffer de 1024 échantillons, on ne calcule que 8 maxima par seconde avec une précision d'environ 8 Hz. Il faut donc maintenir la note suffisamment longtemps pour que 'ce soit clair'. Si on réduit la taille du buffer, la discrimination des fréquences sera moins fine. Si on augmente la fréquence d'échantillonnage (...à déterminer)

Le code Parsons


En essayant de décoder la succession des notes par
void print_parsons(float *fp, int len)
{
static int last_one, sound;
int imax;

imax = get_max(fp, len);
if (imax < 100)
{
sound = 0;
return;
}
if (sound)
return; /* already processed */
if (last_one == 0)
fprintf(stdout, "*");
else {
if (imax > (last_one + 10))
fprintf(stdout, "u");
else if (imax < (last_one - 10))
fprintf(stdout, "d");
else fprintf(stdout, "r");
}
last_one = imax;
sound = 1;
fflush(stdout);
}

On obtient le code Parsons au fur et à mesure.

Contour mélodique schématique


C'est curieusement l'algorithme le plus sophistiqué...

On peut voir le code ici.

Entendre et voir de quoi l'on parle...


Une courte présentation multi-média (assez primitive) pour fixer les idées.




Elle a été réalisée avec Imagination 2.1.

dimanche 28 novembre 2010

Lecteurs MP3



Le Xiron MP950 ressemble très fort au Sweex MP470. 4GB tous les deux. Le Sweex joue aussi les fichiers FLAC et ogg. Le Xiron a un comportement un peu curieux quant à la navigation (je dirais franchement buggé), il faut chipoter pour aller là où l'on veut. Le Sweex ne semble pas passer automatiquement au répertoire suivant. (Dans les versions que j'ai).

Le soft est différent. Le hardware est différent :

Le Xiron :


Le Sweex :


Tous deux utilisent un microContrôleur AK2027C de Artek avec un core Z80 (simple datasheet.pdf) et une flash de 4 GB.

jeudi 28 janvier 2010

PicKit2...

Quelques expériences sous Linux (Ubuntu - 9.10 - Karmic Koala) avec le PICkit 2 Debug Express de Microchip. L'ensemble qui coûte moins de 50.00 € se compose du PICkit 2 avec son interface USB qui permet de programmer la carte d'expérimentation avec un PIC16F887, huit LEDs, un bouton poussoir, un potentiomètre,... Côté logiciel, j'utiliserai uniquement des logiciels libres : pk2cmd, sdcc et gputils.


En chantier. (c'est rien de le dire...)

Premiers tests. (voir aussi)

sudo ./pk2cmd -PPIC16F887 -R -T # start the pre-installed program
# program the device
sudo ./pk2cmd -PPIC16F887 -M -F/media/cdrom0/Install/Lessons/44-Pin\ Demo\ Board/02\ Blink/Blink.HEX
#
sudo ./pk2cmd -PPIC16F887 -M -F../PicKit/Install/Lessons/44-Pin\ Demo\ Board/02\ Blink/Blink.HEX
sudo ./pk2cmd -PPIC16F887 -R -T



Curieusement, le lendemain la LED 'busy' clignote lentement (slow busy blink) et le PicKit2 ne semble plus rien vouloir entendre... Le firmware semble être corrompu, la version n'est plus bonne... Heureusement, il suffit de le recharger (le fichier 'PK2V023200.hex' vient avec pk2cmd).

# ./pk2cmd -?V
Executable Version: 1.20.00
Device File Version: 1.55.00
OS Firmware Version: 118.90.90


Operation Succeeded
# ./pk2cmd /DPK2V023200.hex
Downloading OS...
Verifying new OS...
Resetting PICkit 2...
OS Update Successful.

Operation Succeeded
# ./pk2cmd -?V

Executable Version: 1.20.00
Device File Version: 1.55.00
OS Firmware Version: 2.32.00


Operation Succeeded
# ./pk2cmd -R -T -PPIC16F887

Operation Succeeded



En m'inspirant de ce blog, un petit programme...

#include <pic16f887.h>

/* inspired by http://jsmerritt.blogspot.com/2009/06/getting-started-with-pic.html */
/* ------------------ */
/* Configuration bits */
typedef unsigned int word;
word at _CONFIG1 CONFIG1 = _LVP_OFF & _FCMEN_OFF & _IESO_OFF &
_BOR_OFF & _CPD_OFF & _CP_OFF & _MCLRE_OFF &
_PWRTE_ON & _WDT_OFF & _INTRC_OSC_NOCLKOUT;
word at _CONFIG2 CONFIG2 = _WRT_OFF & _BOR21V;

void isr() interrupt 0
{
/* interrupt service routine */
}

void foo()
{
int i;

for (i = 0; i < 40; i++)
;
}
void main()
{
long i;
unsigned char state;

RP0 = 1;
TRISD = 0; // set PORTD as all outputs
RP0 = 0;

// counter leds
state = 1;
while (1)
{
for (i = 0; i < 2000; ++i)
foo();
state += 1;
PORTD = state;
}
}

Il n'y a plus qu'à...

$ sdcc -mpic14 -p16f887 count.c
$ sudo pk2cmd -PPIC16F887 -M -F count.hex
$ sudo pk2cmd -PPIC16F887 -R -T



Pour programmer la PIC18F4550 USB Development Board de Futurlec, il faut un adaptateur RJ12 (RJ11 à 6 fils) - ICSP (pdf). Il existe un kit à ~7.00 euros chez Microchip, mais si on ajoute les frais (port et taxes), cela devient hors de prix... Fabriquer l'adaptateur est facile ...si on a le matériel (notamment la pince à servir, le bout de câble et un connecteur (type 'header' en ligne pour jumpers).



On remarquera l'inversion 1<->6 (le câble est 'twisté').