samedi 8 octobre 2016

Reflex.ino

Petit exemple pour illustrer la facilité de prise en main d'un Arduino (nano v3.0) (avec un microcontrôleur Atmel ATmega328P). Après en avoir acquis un auprès d'un revendeur officiel ou d'un vendeur de clones chinois, installer le package Arduino sur Ubuntu avec votre gestionnaire de paquets favori (par exemple, Synaptic ou 'sudo apt-get install arduino'). Connectez l'Arduino sur un port USB, lancez l'IDE Arduino sur le PC (via les menus ou 'arduino' en ligne de commandes). Configurez le bidule ([Outils]->[Type de cartes]->[Nano v3.0] et [Outils]->[Port série]..., quelque chose comme /dev/ttyUSB0, /dev/ttyS0,... (regarder la fin de 'dmesg(1)' après avoir connecté le port USB).

Tester si cela fonctionne : [Fichier]->[Exemples]->[0.1 Basic]->[Blink]. En principe, c'est l'application qui tourne déjà quand on achète un Arduino, après peu de temps, une LED clignote ('blink' en anglais) dès que l'on met la carte sous tension. Il faut donc modifier un peu les timings pour observer un changement. Par exemple, remplacer le 1000 millisecondes par 200. Compiler et télécharger avec les boutons de l'IDE ('V' et '->'). Et voilà!...

Maintenant, ce serait bien de faire quelque chose d'original... Par exemple, un truc simple, c'est de tester nos 'réflexes' ('temps de réaction' est plus correct). Le scénario est le suivant : on va allumer la LED du blink après un temps aléatoire et mesurer le temps qu'il faut pour presser un bouton. On a besoin de connecter un bouton à la carte. On met le bouton entre le port 7 (par exemple) et la masse (GND). Pour éviter que le port ne soit 'flottant' (indéterminé), en plus de le programmer en entrée, on y écrit un état haut pour activer la pull-up intégrée au port.

Note : pour pouvoir connecter des fils au 'nano' ou l'utiliser sur un breadboard, il faut souder...

Fritzing (également disponible sous Ubuntu) nous permet de réaliser le schéma suivant :
On voit que c'est la simplicité même : deux fils, un bouton. On utilise les pins D7 (port 7) et GND (masse). La LED est celle qui équipe déjà la carte.

Reste à écrire le programme... Un petit tour sur https://www.arduino.cc/en/Reference/HomePage ou sur les exemples locaux pour savoir ce que nous avons à disposition. Typiquement, 'Serial.begin()' pour configurer l'impression via le lien série/USB sur une console sur le PC, random() pour générer un temps aléatoire et 'Serial.println()' pour afficher le temps.

/*
**  Reflex
**        la LED s'allume après un délai variable
**        on mesure le temps qu'il faut pour appuyer sur un bouton
*/
 
int led = 13;        /* port de la LED sur des cartes Arduino */
int bouton=7;        /* port connecté au bouton connecté à la masse */

void setup()
        {                
        pinMode(led, OUTPUT);
        pinMode(bouton, INPUT);
        digitalWrite(bouton, HIGH); /* activation de la pull-up */
        Serial.begin(9600);
        }

void loop()
        {
        int cnt;        /* ms entre l'allumage et le clic */

        digitalWrite(led, LOW);    /* démarrage LED éteinte */
        delay(random(2000,5000));  /* delai aléatoire 2..3 seconds */
        digitalWrite(led, HIGH);   /* allumage de la LED */
 
        for (cnt = 0; cnt < 2000; cnt++)
                { /* comptage en millisecondes */
                delay(1);
                if (digitalRead(bouton) == LOW)
                        break;     /* interruption de la boucle si le bouton est pressé */
                }
        Serial.println(cnt);       /* affichage du nombre de millisecondes */
        }


Reste à tester le programme. On le compile, on le charge. On a maintenant besoin d'une console, un terminal connecté au microcontrôleur. [Outils]->[Moniteur série]. On vérifie qu'il est sur le bon /dev/ttyUSB0 (ou autre; c'est le même que celui qui sert à programmer l'Arduino et que la vitesse est celle programmée dans le programme (ici, 9600 bauds). Il ne nous reste plus qu'à appuyer sur le bouton lorsque la lumière s'allume.



On peut aussi utiliser un autre moniteur comme 'minicom(1)' ou 'screen(1)', enregistrer ces résultats et faire des statistiques, calculer des moyennes, des écarts-types, faire des graphes, des histogrammes,... Le faire à différents moments de la journée, après avoir bu du café ou du vin, en téléphonant histoire de vérifier tout ce que l'on nous raconte sur la sécurité routière. On peut aussi tester si notre vigilance diminue après 2 minutes ou deux heures, quelle est notre endurance. Comparer les résultats avec ceux d'autres personnes, des amis, des membres de notre famille. Est-ce que l'entraînement apporte quelque chose? Est-ce que la situation se dégrade avec le temps? On est peut-être plus fringant à 20 ans qu'à 60... Et, est-ce que, si on déplace le bouton vers le pied, le temps s'allonge? (Dans les voitures le frein est aux pieds et l'influx nerveux a plus de chemin à parcourir pour aller de la tête aux pieds que de la tête à la main; les grands sont-ils moins réactifs?) Bref, on peut faire plein de choses avec ces deux bouts de fils et ces 20 lignes de code. Et, dans ce cas précis, je ne connais pas de produit du commerce accessible qui me permettent d'obtenir pareil résultat. L'imagination est la seule limite. (NB: dans ce cas-ci, une application sur PC ou smartphone pourrait le faire, mais ici, on peut remplacer le bouton par un capteur ou augmenter sensiblement la précision)

Le programme est, bien sûr, primitif. Il faut savoir que quand il affiche 2000, c'est qu'on n'a pas appuyé sur le bouton et que s'il affiche 0, c'est que l'utilisateur a appuyé avant d'avoir vu la LED s'allumer... Et, pour être sûr, il faudrait peut-être aussi l'étalonner.

Comme test extrême, vous pouvez imiter Gaston Lagaffe qui teste un nouveau type d’alcootest qu'il vient d'inventer. Il arrive au bureau avec plein de bouteilles d'alcool et augmente la dose afin que son invention, une cigarette-alcootest, vire au rouge... (planche 658; album 10 1977 "Le Géant de la Gaffe" page 5)

jeudi 19 novembre 2015

Echo sur Arduino nano sans IDE

En utilisant le strict minimum, à partir de la datasheet de l'ATmega368 :
#include <avr/io.h>
#define F_CPU 16000000UL
#include <util/delay.h>
/*
** UART
*/
void uart_init(void)
    {
    UBRR0H = 0x00;
    UBRR0L = 103;  /* F_CPU/16/baud-1 for 9600 bps*/
    UCSR0C = (1<<USBS0) | (3<<UCSZ00); /* 8N1 */
    UCSR0B = (1<<RXEN0) | (1<<TXEN0); /* enable */
    }
char my_getc()
    {
    while (!(UCSR0A & (1<<RXC0)))
        ;
    return(UDR0);  /* should handle errors */
    }
void my_putc(char c)
    {
    while (!(UCSR0A & (1<<UDRE0)))
        ;
    UDR0 = c;
    }
/*-----------------------------------------------*/
int main(void)
    {
    uart_init();

    for (;;)
        my_putc(my_getc());
    }
Avec comme commandes :
$ avr-gcc  -mmcu=atmega328 -Os -o echo.elf echo.c
$ avrdude -c arduino -p atmega328P -P /dev/ttyUSB0 -b 57600 -U flash:w:echo.elf
On peut tester avec
$ minicom -D /dev/ttyUSB0 -b 9600

AVR libc

Intrigué par le fait que mes putc() et getc() entraient en conflit avec des 'built-in', j'ai fini par me rendre compte de l'existence d'une véritable libc pour AVR... Et c'est quand même vachement mieux fichu que l'API Arduino. Ouf! Notamment du côté de stdio. Adapté, le programme devient quelque chose comme :
#include <avr/io.h>
#define F_CPU 16000000UL
#include <util/delay.h>
#include <stdio.h>
#include <string.h>
/*
** UART
*/
void uart_init(void)
    {
    UBRR0H = 0x00;
    UBRR0L = 103;  /* 9600 bps | F_CPU/16/baud-1 */
    UCSR0C = (1<<USBS0) | (3<<UCSZ00); /* 8N1 */
    UCSR0B = (1<<RXEN0) | (1<<TXEN0);
    }
int uart_getchar(FILE *stream)
    {
    while (!(UCSR0A & (1<<RXC0)))
        ;
    return(UDR0);
    }
int uart_putchar(char c, FILE *stream)
    {
    while (!(UCSR0A & (1<<UDRE0)))
        ;
    UDR0 = c;
    }
/*-----------------------------------------------*/
/* line input with minimal editing */
char *getline(char *cp, int n)
    {
    int    i=0;
    char    c;

    for(;;)
        {
        c = fgetc(stdin);
        if (c=='\n' || c=='\r' || i==(n-1))
            {
            cp[i] = '\0';
            return(cp);
            }
        if (c=='\b' && i>0)
            {
            fputc(c, stdout);
            fputc(' ', stdout);
            i -= 1;
            }
        else    cp[i++] = c;
        fputc(c, stdout);
        }
    }
int main(void)
    {
    char buf[32];

    uart_init();
    fdevopen(uart_putchar, uart_getchar);
    printf("Hello World\n\r");
    for (;;)
        {
        printf("\n\rReady> ");
        getline(buf, sizeof(buf));
        printf("\n\rstrlen('%s') is %d.\n\r", buf, strlen(buf));
        }
    }

lundi 6 juillet 2015

Digispark - ATTiny85

Digistump propose un petit gadget USB avec un AVR ATtiny85 directement relié au bus USB, le Digispark. Le micro-contrôleur se charge d'interpréter les signaux USB-1 (Étant OpenSource, on en trouve des clones).

Dans la même veine que l'article précédent, voici comment le programmer simplement.

Sur Ubuntu, il faut probablement les pré-requis suivants :
$ sudo apt-get install gcc-avr binutils-avr avr-libc
$ git clone https://github.com/micronucleus/micronucleus.git
$ ...
(en fait, je ne suis pas tout-à-fait sûr parce que j'ai aussi fait un 'apt-get install arduino'...)
Ensuite, le même petit programme 'C' :
#include <avr/io.h>
#define F_CPU 16000000UL
#include <util/delay.h>

int main(void)
    {
    DDRB  |= _BV(PB1);

    for (;;)
        {
        _delay_ms(1000);
        PORTB ^= _BV(PB1);
        }
    }
Il ne reste plus qu'à préparer le fichier .hex et le télécharger avec 'micronucleus' (micronucleus comporte deux parties : le firmware installé sur le chip en fin de FLASH (qui utilise V-USB) et la commande 'micronucleus' qui tourne sous Linux et utilise libusb).
$ avr-gcc  -mmcu=attiny85 -Os -o blink.elf blink.c
$ avr-objcopy -j .text -j .data -O ihex blink.elf blink.hex
$ micronucleus --run blink.hex
Il faut en fait lancer la commande 'micronucleus' puis introduire le bidule dans un port USB. Le flashage s'effectue et le programme commence à s'exécuter. La LED du Digispark (connectée à PB1) se met à clignoter (j'avais utilisé PB3 mais c'est probablement une mauvaise idée parce le port est utilisé pour communiquer sur l'USB). Notez que F_CPU est à 16 MHz. Cette vitesse est tirée de l'USB par le firmware micronucleus.

Le programme est compilé, par défaut, pour s'exécuter en début de FLASH. Le bootloader (firmware micronucleus) se situe en fin de FLASH; il s'exécute lors d'un 'power-on-reset', c'est lui qui flashe le programme au bon endroit et qui lui transfère l'exécution après cinq secondes. Il est conçu en sorte de ne pas s'écrire dessus, il n'y a donc aucun risque de 'briquer' le Digispark.

Maintenant, comment faire pour créer un firmware communiquant en USB? Soit partir de la documentation V-USB et/ou partir d'un firmware qui fonctionne comme Little Wire (ou une des nombreuses applications signalées sur le site de V-USB).

jeudi 7 mai 2015

ATTiny13 blink

Le but est de faire clignoter une LED avec un nombre réduit de composants, un environnement réduit et un programme simple... Dans ce cas-ci, le développement a lieu sur Ubuntu 14.04 mais cela devrait pouvoir se faire sur n'importe quel Linux, Windows,...

Arduino, c'est magique. Un peu trop magique... C'est quand même mieux de savoir ce que l'on fait et de s'approcher de la réalité des choses. Et donc, par exemple, d'utiliser un programmeur ISP (qui utilise un bus SPI...) pour programmer l'AVR directement plutôt que de passer par le bootloader installé dans l'Arduino. Au niveau du prix, quand on veut le faire soi-même, cela revient à peu près au même. Un clone de USBasp revient au même prix qu'une interface USB-série, quelque chose comme 2 euros. En plus, on n'est plus limité à l'ATMega368 (ou 168), on peut utiliser toute la gamme AVR, comme ici avec l'ATTiny13. On peut aussi se servir d'un Arduino comme programmateur avec ArduinoISP.



Il y a donc quelques fils à connecter : le 5 (ou 3?) volts (VCC), la terre (GND), les trois fils du SPI (MOSI, MISO, CLK) et le 'chip select' (c.-à-d. /RESET).

Le programme pour faire clignoter une LED, sans utiliser l'API Arduino (mais avr-libc), se réduit à
#include <avr/io.h>
#include <util/delay.h>

int main(void)
    {
    DDRB  |= _BV(PB3);
 
    for (;;)
        {
        _delay_ms(200);
        PORTB ^= _BV(PB3);
        }
    }
On programme le port3 en sortie et on joue sur le bit data correspondant pour faire clignoter la LED. (Notons que le compilateur va se plaindre que F_CPU n'est pas défini dans delay.h (il prend alors la valeur par défaut, 1 MHz); on pourrait ajouter un #define F_CPU 1000000UL avant l'include ou passer F_CPU au compilateur dans le Makefile)

Le Makefile, sans bell ni whistle, comprend 3 étapes : compiler avec avr-gcc, transformer le binaire obtenu en fichier texte (iHEX) et le téléverser avec avrdude. (On utilise les paquets gcc-avr, binutils-avr, avr-libc, avrdude; la plupart sont déjà requis par le paquet arduino)

blink.elf: blink.c
        avr-gcc  -std=c99 -Wall -Os -mmcu=attiny13 -o blink.elf blink.c

blink.hex: blink.elf
        avr-objcopy -j .text -j .data -O ihex blink.elf blink.hex

flash: blink.hex
        avrdude -c usbasp-clone -p t13 -U flash:w:blink.hex

Pour peu que l'ont ait créé un /etc/udev/rules.d/99-USBasp.rules avec

SUBSYSTEM=="usb", ATTR{product}=="USBasp", ATTR{idProduct}=="05dc", ATTRS{idVendor}=="16c0", MODE="0666"

dedans, on peut simplement faire
$ make flash
et obtenir une LED qui clignote sur PB3.



Une fois enlevé le cordon ombilical de programmation, il ne reste que l'alimentation, le microcontrôleur, une résistance de 10K entre /RESET et VCC, un condensateur de 100nF entre VCC et GND et une 'blinkenlight' entre PB3 et GND (dans ce cas-ci une LED et une résistance de 1K (c'est beaucoup mais selon AVR Programming, comme ça cela fonctionne aussi avec 12 volts).

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é').

vendredi 2 octobre 2009

Loi de Zipf sous Linux

Vers 1935, George Kingsley Zipf entreprit d'analyser la fréquence des mots dans Ulysses de James Joyce...

Après avoir classé les mots par leur fréquence d'apparition, il se serait aperçu que si le mot le plus fréquent apparaissait N fois, le dixième mot le plus fréquent apparaissait N/λ, le centième mot, N/λ², etc...

Il est aujourd'hui facile de vérifier ses travaux.

D'abord, il faut aller chercher Ulysses sur le Projet Gutenberg et en traitant le texte avec des outils standards Unix/Linux, on obtient assez facilement le résultat attendu.

Il suffit, en effet, d'envoyer le texte dans une série de filtres :

cat 4300-8.txt |

La première chose à faire est de transformer toutes les majuscules en minuscules pour que 'The' soit équivalent à 'the' :

|tr [A-Z] [a-z]

Ensuite, on sépare tous les mots pour en mettre un sur chaque ligne. Le plus simple est de substituer tous les caractères non alphabétiques par un retour à la ligne :

|sed 's/[^a-z]/\n/g'

Tant qu'à faire, supprimons toutes les lignes vides créées :

|awk '/[a-z]/{print $1;}'

Nous avons maintenant une liste de mots qu'il convient de trier :

|sort

Et compter chaque mot :

|uniq -c

Nous ne nous intéressons pas aux mots en eux-mêmes, mais au nombre de fois qu'ils apparaissent :

|awk '{print $1;}'

On va trier cette liste de fréquences par ordre descendant pour avoir la plus haute fréquence au premier rang (première ligne) :

|sort -rn

Et numéroter les lignes pour obtenir deux colonnes : le rang et la fréquence :

|pr -n -t

En résumé :

time cat 4300-8.txt|tr [A-Z] [a-z]|sed 's/[^a-z]/\n/g'|awk '/[a-z]/{print $1;}'|sort|uniq -c|awk '{print $1;}'|sort -rn|pr -n -t > zipf.data

Cela a pris, en gros, 9 secondes, bien moins que ce que le pauvre monsieur Zipf a dû y consacrer.

Si on trace le graphe log/log résultant avec Gnuplot

set logscale
set xlabel 'rang'
set ylabel 'frequence'
set title 'Loi de Zipf (Ulysses)'
set term png
set output 'UlyssesZipf.png'
plot 'zipf.data'

On obtient :




Un autre truc amusant (qui n'a probablement rien à voir avec Zipf) est la génération d'un graphique log/log qui montre 'le nombre 'y' de mots présents 'x' fois'.

À partir de la liste des fréquences, je compte combien de mots apparaissent un fois, deux fois,...

|sort -n|uniq -c

Et, comme on veut obtenir un graphique : nombre 'y de fois' qu'un mot apparaît 'x fois', il faut inverser les colonnes avant de l'injecter dans Gnuplot :

|awk '{print $2" "$1}'

En résumé :

time cat 4300-8.txt|tr [A-Z] [a-z]|sed 's/[^a-z]/\n/g'|awk '/[a-z]/{print $1;}'|sort|uniq -c|awk '{print $1;}'|sort -n|uniq -c|awk '{print $2" "$1;}' > data

Il ne reste plus qu'à 'gnuploter' tout ça :

set logscale
set term png
set output 'UlyssesZipf.png'
plot 'data'

Et voilà!

On 'voit bien' que sur un graphique log/log,
le nombre 'y' de fois qu'un mot apparaît 'x' fois
suit une loi de puissance.

mardi 9 juin 2009

Boule


Pour faire une boule, il suffit de frotter deux corps en rotation. (?)

Catadioptre


Une pile de cubes chromés formant un 'catadioptre'...

dimanche 7 juin 2009

Coin de cube avec POV-Ray


La boule rouge est placée derrière la caméra, le ciel est bleu, le sol est sable; on tourne autour d'un cube chromé dont on a retiré un 'coin de cube'... On remarquera que la boule (caméra) restent bien au centre quand elle se trouve dans le coin de cube.

#include "textures.inc"
#include "colors.inc"

global_settings {
assumed_gamma 1.5
noise_generator 2
}

light_source {
<-10, 100, -10>, rgb <1, 1, 1>
}
light_source {
<30, 20, -80>, rgb <1, 1, 1>
}
background { color rgb <0.4, 0.4, 0.6> }

camera {
location <8, 4, -4>
look_at <0, 0, 0>
rotate <0, 0, 15>
rotate <0, -15+clock*90, 0>
}
sphere {
<32, 16, -16>, 1
pigment { color Red}
rotate <0, 0, 15>
rotate <0, -15+clock*90, 0>
}
#declare D = 0.00001; // just a little bit !
difference {
box {
// big cube
<0, 0, 0>, <3, 3, -3>
texture{Polished_Chrome
finish {reflection 0.9}}
}
box {
// creux
<2, 2, -2>, <3+D, 3+D, -3-D>
texture{Polished_Chrome
finish {reflection 0.9}}
}
}
plane{ <0,1,0>, 0
texture{
pigment {color rgb <0.85,0.6,0.4>}
} // end of texture
translate <0, -2, 0>
} // end of plane


----

; Persistence Of Vision raytracer version 3.5 CubeCorner.ini file.

Antialias=On
Antialias_Threshold=0.3
Antialias_Depth=3

Input_File_Name=CubeCorner.pov

Initial_Frame=1
Final_Frame=30
Initial_Clock=0
Final_Clock=1

Cyclic_Animation=on
Pause_when_Done=off

----

povray -H200 -W200 CubeCorner.ini
convert -delay 50 -loop 0 -dispose Previous CubeCorner*.png CubeCornerMirror.gif

lundi 4 mai 2009

Faire un baudrier avec une sangle (et un peu de ferraille)



Il faut environ 2.50 mètres de sangle plate de 45 mm (noire sur le dessin), deux 'boucles' triangulaires (rouge) et deux 'boucles' carrées avec deux trous (bleu). Il est/était connu, à Bruxelles, sous le nom 'baudrier Lecomte' (du nom d'une famille de vendeurs de matériel d'escalade/spéléo/...). Lorsque la sangle est usée, il est facile et économique de la remplacer. Le confort n'est pas exceptionnel, mais il est simple, solide et efficace. Le schéma ci-dessus donne le parcours de la sangle dans les boucles. Le bas de la sangle passe derrière le dos; on passe les jambes dans les deux grandes boucles (reliée par une ficelle (verte)). Les 'triangles rouges' sont maintenus contre les plaquettes bleue par la sangle (la boucle est la plus courte possible); les 'jambières' passent en dessous. Le 'delta' (mousqueton en acier, triangulaire, à vis (non représenté)) ferme les deux 'triangles rouges'.

dimanche 14 décembre 2008

Mal lunée...

Il arrive que la Lune soit mal représentée. Par exemple, on ne peut pas avoir de derniers croissants le soir, c'est une Lune du matin... Ce type d'erreur est fréquente dans les épisodes des Simpsons. Par exemple, dans l'épisode GABF05, où l'on voit Homer dire bonsoir à Bart au fond du jardin parce qu'il est l'objet d'une mesure d'éloignement de sa soeur Lisa.
(aussi dans GABF12, GABF20, CABF20, DABF01, DABF02, DABF8, DABF11, DABF16, DABF18, DABF21... dans le film la Lune change de phase d'un moment à l'autre lors du voyage de retour d'Alaska; dans l'épisode 'Qui veut tuer Homer?' (EABF01), on a un dernier croissant quand Tahïti Bob sonne à la porte des Simpson et une pleine Lune lorsqu'il se retrouve dans la chambre de Bart). De même dans HABF02 :



Mais, cela arrive aussi sur les boîtes aux lettres de la poste belge :



Ou, sur une affiche des Grignoux :

où la photo est prise en direction de l'ouest...

Sur la météo de la télévision locale rtc...