[PYTHON] Observez le groupe de météores Futago avec RaspberryPi4

Aperçu

Depuis que j'ai acheté Raspai 4, je parlerai de la réalisation d'un système d'observation des ondes radio d'étoile filante avec Raspai. Aussi, puisque le 12 au 15 décembre est le temps du groupe de météores de Futago chaque année, j'ai essayé de voir s'il pouvait réellement être observé.

Ce sera une longue histoire car elle sera expliquée à partir du principe.

Principe d'observation

Des étoiles filantes apparaissent lorsque la poussière flottant dans l'espace saute dans l'atmosphère terrestre, entre violemment en collision et émet du plasma. Le gaz appelé colonne d'ionisation généré à ce moment-là a la propriété de réfléchir les ondes radio, et en utilisant ce phénomène, il est possible de savoir qu'un météore a coulé.

Comme le montre la figure ci-dessus, les ondes radio d'une distance qui normalement n'atteignent pas physiquement sont réfléchies et n'atteignent dans l'atmosphère qu'au moment où le météore coule.

système

Cette fois, nous recevrons l'onde de balise 53,755MHz de l'Université préfectorale de Fukui avec SDR (radio logicielle). Les ondes radio provenant de l'antenne sont converties en son par SDR, le son est émis par l'audio USB, puis placé dans l'entrée du microphone tel quel, et la FFT est effectuée pour analyser l'écho du météore.

Les choses nécessaires

・ Raspberry Pi4 ・ Câble de conversion audio USB · Câble audio ・ AirspyMini (Software Radio) ・ Antenne bande 50 MHz ・ Câble coaxial ・ Divers connecteurs (Dans le cas de l'antenne ci-dessus, elle devient un connecteur de type M, donc une conversion de type M-SMA est nécessaire)

Un total d'environ 30 000 à 40 000 yens. Si la méthode utilise un récepteur analogique comme par le passé, un récepteur analogique et un PC sont nécessaires, le coût est donc probablement inférieur car il s'agit de SDR + rasp pie. Si vous pensez que les étoiles filantes sont un système de rêve qui réalisera votre souhait, 30000 yens ne devraient pas être si chers.

Vous pouvez également le mettre dans une boîte et l'utiliser pour ne pas vous blesser le cœur même si vous le laissez à l'extérieur dans le froid.

Installation

Tout ce que vous avez à faire est d'assembler et de vous connecter.

Au lieu de vous soucier des yeux glacés du quartier, attachez simplement l'antenne à un endroit élevé et pointez-la vers le ciel cible (Fukui).

L'élément le plus court est la direction avant.

installer

Je vais préparer la tarte aux râpes.

** 1. Activer SDR ** Cette fois, nous utiliserons gqrx, un logiciel SDR, pour l'exécuter sur l'interface graphique.

#téléchargement gqrx
$wget https://github.com/csete/gqrx/releases/download/v2.6/gqrx-2.6-rpi3-2.tar.xz

#Déploiement
$tar xvf gqrx-2.6-rpi3-2.tar.xz

$cd gqrx-2.6-rpi3-2/

$./setup_gqrx.sh

#Courir
$sh run_gqrx.sh

Au début, l'écran de configuration suivant apparaît.

・ L'entrée I / Q est "AirSpy AIRSPY" ・ Le taux d'entrée est de 3000000 Ensuite, les paramètres par défaut devraient convenir.

Essayons de voir si vous pouvez entendre quelque chose en fonction de la fréquence de la radio FM.

** 2. Installation de la bibliothèque utilisée pour l'analyse sonore **

#Dessin graphique
$sudo apt install python3-scipy
$sudo pip3 install drawnow

#Système audio
$sudo apt-get install python3-pyaudio
$sudo apt-get install lame

Effectuer SDR

Le réglage du SDR nécessite un petit réglage détaillé.

Ouvrez les Options du récepteur à partir de l'écran des paramètres sur la droite et ouvrez ・ La largeur du filtre est étroite ・ La forme du filtre est nette ・ Mode vers USB

À Le mode fait référence à la méthode de modulation USB, qui coupe et module les composantes de fréquence inférieures à la fréquence définie en réception.

Ensuite, en ce qui concerne la spécification de fréquence supérieure gauche, il est nécessaire de décaler légèrement la fréquence de réception de la fréquence d'émission afin de produire du son.

Lorsqu'il est combiné à la fréquence 53,754 MHz, qui est 1000 Hz inférieure à la fréquence côté émission de 53,755 MHz Le son de 1000 Hz sortira.

** [Fréquence de transmission-Fréquence de réception = Fréquence sonore] **

Analyser l'écho des météores

Si le SDR fonctionne correctement, lorsqu'un météore coule, un téléphone ♪ et un son aigu d'environ 1000 Hz seront émis. Lorsque la réflexion des ondes radio est longue, le son peut continuer à être entendu pendant 30 secondes ou plus.

↓ Comme ça. Attention au volume! https://soundcloud.com/user-395229817/20191218-175747a

Nous analyserons ce son avec Python.

streamfft.py

# -*- coding: utf-8 -*-

from subprocess import getoutput

import argparse
from scipy.fftpack import fft

import numpy as np
import itertools
import matplotlib
matplotlib.use("TkAgg")
import matplotlib.pyplot as plt
import pyaudio
import wave
from datetime import datetime, timedelta


#Fréquence à analyser(Hz)
FREQ_LOW = 800
FREQ_HIGH = 1200

#Volume maximum
VOL_LIMIT = 2800


#Seuil de détection 0~255
THRESHOLD = 70

#Heure d'enregistrement au moment de la détection(s)
RECORD_SECONDS = 30

#Destination de sauvegarde des données d'enregistrement
SAVE_PATH = '/PATHtoSAVE/'

def str2bool(v):
    if v.lower() in ('true', 't', 'y', '1'):
        return True
    elif v.lower() in ('false', 'f', 'n', '0'):
        return False
    else:
        raise argparse.ArgumentTypeError('Boolean value expected.')

#Que ce soit pour afficher un graphique
dataplot = False
parser = argparse.ArgumentParser(description="data plot enable")
parser.add_argument("-p",
                    metavar="--p",
                    type=str2bool,
                    nargs="+",
                    help="data plot enable"
                    )


#S'il faut enregistrer le MP3 lorsqu'un météore est détecté
datasave = False
parser.add_argument("-s",
                    metavar="--s",
                    type=str2bool,
                    nargs="+",
                    help="MP3 Save"
                    )

args = parser.parse_args()

dataplot = bool(args.p[0])
datasave = bool(args.s[0])
print("data plot enable : " + str(dataplot))
print("MP3 Save : " + str(datasave))


class SpectrumAnalyzer:
    FORMAT = pyaudio.paInt16
    CHANNELS = 1
    RATE = 44100
    CHUNK = 22050
    N = 22050

    #Indicateur de détection
    detect = False

    #Drapeau d'enregistrement
    isRecording = False

    data = []
    frames = []
    y = []
    maxVol = 0

    startTime = datetime.now()
    endTime = startTime + timedelta(seconds=RECORD_SECONDS)
    filename = startTime.strftime("%Y%m%d_%H%M%S")

    def __init__(self):

        self.audio = pyaudio.PyAudio()
        DEVICE_IDX = 0
        for i in range(self.audio.get_device_count()):
            dev = self.audio.get_device_info_by_index(i)
            if dev['name'].find('USB PnP Sound Device') != -1:
                DEVICE_IDX = i

        self.stream = self.audio.open(
            format = self.FORMAT,
            channels = self.CHANNELS,
            rate = self.RATE,
            input_device_index = DEVICE_IDX,
            input = True,
            output = False,
            frames_per_buffer = self.CHUNK,
            stream_callback = self.callback)

        self.stream.start_stream()
        print('SpectrumAnalyzer init')


    def close(self):
        self.stream.stop_stream()
        self.stream.close()
        self.audio.terminate()
        print("Quitting...")

    def callback(self, in_data, frame_count, time_info, status):
        yf = fft(np.frombuffer(in_data, dtype=np.int16)) / self.CHUNK
        self.y = np.abs(yf)[1:int(self.CHUNK/2)]

        freqRange = self.RATE / self.N
        lowIdx = int(FREQ_LOW / freqRange)
        highIdx = int(FREQ_HIGH / freqRange)
        vol = self.y[lowIdx:highIdx]

        self.maxVol = 0
        idx = 0
        for v in vol:
            #Le bruit blanc est de 5~10%Ajusté à 25 ou moins
            vol[idx] = v.round()

            # VOL_Couper avec LIMIT
            if VOL_LIMIT <= vol[idx]:
                vol[idx] = VOL_LIMIT

            #Normalisé à 255
            vol[idx] = int(round(vol[idx] / VOL_LIMIT * 255, 0))

            #Valeur maximum
            if (self.maxVol <= vol[idx]):
                self.maxVol = vol[idx]

            if (vol[idx] > 255):
                vol[idx] = 255

            idx = idx + 1


        if(dataplot):
            self.data = vol.tolist()
            self.drawGraph()

        #Comparaison des seuils
        if THRESHOLD <= self.maxVol:
            self.detect = True
        else:
            self.detect = False

        #Commencer l'enregistrement
        if self.isRecording == False and self.detect:
            self.isRecording = True
            self.startRecording()

        #Traitement pendant l'enregistrement
        if self.isRecording:
            self.frames.append(in_data)

            #Traitement à la fin de l'enregistrement
            if datetime.now() > self.endTime:
                self.isRecording = False
                self.stopRecording()


        return(None, pyaudio.paContinue)


    def drawGraph(self):
        plt.subplot(1,1,1)
        plt.clf()
        plt.plot(self.data)
        plt.draw()
        plt.pause(0.05)

    def startRecording(self):
        self.startTime = datetime.now()
        self.endTime = self.startTime + timedelta(seconds=RECORD_SECONDS)
        self.filename = self.startTime.strftime("%Y%m%d_%H%M%S")

        print(self.startTime.strftime("%Y%m%d_%H%M%S") + " Recording Start")


    def stopRecording(self):
        waveFile = wave.open(SAVE_PATH + self.filename + '.wav', 'wb')
        waveFile.setnchannels(self.CHANNELS)
        waveFile.setsampwidth(self.audio.get_sample_size(self.FORMAT))
        waveFile.setframerate(self.RATE)
        waveFile.writeframes(b''.join(self.frames))
        waveFile.close()

        self.frames = []

        print("Recording Stop")

        #Conversion MP3
        self.convertMp3()

    def convertMp3(self):
        checkConvert = getoutput("lame -b 128 " + SAVE_PATH + self.filename + '.wav' + " " + SAVE_PATH + self.filename + '.mp3')
        print(checkConvert)

        #Suppression du WAV
        srcDel = getoutput("rm -rf " + SAVE_PATH + self.filename + '.wav')
        print(srcDel)


def keyInterrupt():
    while True:
        try:
            pass
        except KeyboardInterrupt:
            spec.close()
            return

spec = None
if __name__ == "__main__":
    spec = SpectrumAnalyzer()
    keyInterrupt()

En tant que flux de traitement

** 1. FFT toutes les 0,5 seconde 2. Extraire les données de 800 à 1200 Hz 3. Déterminez si le volume dans la plage ci-dessus dépasse le seuil 4. S'il dépasse, enregistrez-le sous wav pendant 30 secondes 5. Convertir wav en mp3 **

Avec ce sentiment, un fichier mp3 sera créé à chaque fois qu'une étoile filante coule.

Ajustement du seuil

Affichez un spectrogramme pour ajuster le seuil de détection des météores.

#Dessin spectrogramme
# p → True
python3 /home/pi/streamfft.py -p True -s False

Étant donné que la FFT est effectuée toutes les 0,5 seconde, la résolution de fréquence est de 2 Hz. Puisque nous analysons 400Hz de 800 à 1200Hz, nous pouvons obtenir 200 points de données dans une FFT. Je suis désolé de ne pas avoir d'étiquette, mais l'axe X est de 800 Hz à gauche et de 1200 Hz à droite. L'axe Y est une valeur normalisée de 0 à 255.

Dans la figure ci-dessus, le volume VOL_LIMIT est ajusté de sorte que le plancher de bruit soit d'environ 10 sur 255 avec seulement du bruit blanc. Après avoir effectué ce réglage, si le seuil de détection est dépassé, c'est un météore! C'est un jugement approximatif.

↓ Lorsque l'onde radio du météore entre, elle augmente comme indiqué sur la figure ci-dessous.

Maintenant à observer!

#Commande d'observation
#s → Réglez sur True et enregistrez les données
python3 /home/pi/streamfft.py -p False -s True

Résultats de l'observation du groupe Météore Futago 2019

L'observation a commencé à 2 heures le 14/12 et s'est terminée à 10 heures le 16/12, donc cela a pris environ 56 heures. Il semble que la vague de balise de Fukui s'est arrêtée pendant environ 12 heures, et aucune donnée n'a pu être obtenue pendant cette période.

Il y avait environ 400 données, et apparemment une autre chose semblable à du bruit se chevauchait ou j'ai joué quelque chose qui n'était pas un écho de météore, et l'impression était d'environ 200. Je ne pouvais pas avoir une bonne impression parce que la condition AGC du SDR et le jugement de seuil fixe ne s'accordaient pas très bien, ou il y avait un temps de perte de données.

Ce sont les données des stations d'observation des étoiles filantes installées dans diverses régions du Japon il y a environ deux ans. La valeur FFTed par Raspeye est envoyée à AWS et l'écho du météore est détecté en temps réel sur le serveur. À l'origine, le nombre de météores a progressivement augmenté à partir d'environ 12/12, et la caractéristique du groupe de météores de la constellation Futago était qu'il a fortement chuté après le pic de 12 / 14-15.

Résumé

Les performances de RaspberryPi4 se sont considérablement améliorées, et il est désormais possible d'effectuer une analyse d'écho tout en SDR. RaspberryPi3 est plein de SDR ...

De plus, cette fois, je l'ai introduit avec un simple jugement quant à savoir si le seuil a été dépassé ou non, mais lors de l'observation réelle, divers bruits seront introduits, donc je pense qu'il serait encore plus intéressant d'ajouter un traitement de jugement avancé pour le repousser.

ex. Doppler est appliqué à l'écho réfléchi par l'avion. Il est jugé par le déplacement de fréquence et repoussé. Toka Toka