[PYTHON] Observe the Geminids meteor shower with Raspberry Pi 4

Overview

Since I bought Raspberry Pi 4, I will talk about making a shooting star radio observation system with Raspberry Pi. Also, since 12/12 ~ 15 is the time of the Gemini meteor shower every year, I tried to see if it could actually be observed.

It will be a long story because it will be an explanation from the principle.

Observation principle

Shooting stars appear when dust floating in space jumps into the Earth's atmosphere, collides violently, and emits plasma. The gas called the ionization column generated at that time has the property of reflecting radio waves, and by using this phenomenon, it is possible to know that a meteor has flowed.

As shown in the figure above, radio waves from a distance that normally do not reach physically are reflected and reach the atmosphere only at the moment when the meteor flows.

system

This time, the beacon wave of 53.755MHz from Fukui Prefectural University will be received by SDR (Software Radio). The radio waves coming in from the antenna are converted into sound by SDR, the sound is output from USB audio, then put into the microphone input as it is, and FFT is performed to analyze the meteor echo.

Things necessary

・ Raspberry Pi4 ・ USB audio conversion cable ・ Audio cable ・ AirspyMini (Software Radio) ・ 50MHz band antenna ・ Coaxial cable ・ Various connectors (In the case of the above antenna, it becomes an M type connector, so M-SMA type conversion is required)

A total of about 30,000 to 40,000 yen. If you use the analog receiver as in the past, you will need an analog receiver and a PC, so the cost will be lower because it is SDR + Raspberry Pi. If you think that a shooting star is a dream system that makes your wish come true, 30,000 yen should not be so expensive.

You can also put it in a box and operate it, so it won't hurt your heart even if you leave it outdoors in the cold.

Installation

All you have to do is assemble and connect.

Instead of worrying about the frosty eyes of the neighborhood, just tie the antenna to a high place and point it at the target (Fukui) sky.

The shorter element is the forward direction.

setup

I will set up the Raspberry Pi.

** 1. Enable SDR ** This time, we will use the SDR software called gqrx to run it on the GUI.

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

#Deployment
$tar xvf gqrx-2.6-rpi3-2.tar.xz

$cd gqrx-2.6-rpi3-2/

$./setup_gqrx.sh

#Run
$sh run_gqrx.sh

At first, the following Configure screen will appear.

・ I / Q input is "AirSpy AIRSPY" ・ Input rate is 3000000 Then, the default settings should be fine.

Let's try to see if you can hear something according to the frequency of the FM radio.

** 2. Installation of library used for sound analysis **

#Graph drawing
$sudo apt install python3-scipy
$sudo pip3 install drawnow

#Audio system
$sudo apt-get install python3-pyaudio
$sudo apt-get install lame

Perform SDR

SDR tuning requires a little finer settings.

Open Receiver Options from the settings screen on the right ・ Filter width is Narrow ・ Filter shape is Sharp ・ Mode is USB

To Mode refers to the USB modulation method, which cuts and modulates frequency components lower than the frequency set in reception.

Next, regarding the frequency specification on the upper left, it is necessary to shift the reception frequency slightly from the transmission frequency in order to produce sound.

When combined at the 53.754MHz frequency, which is 1000Hz lower than the transmitting frequency of 53.755MHz The sound of 1000Hz will come out.

** [Transmission frequency-Reception frequency = Sound frequency] **

Analyze meteor echo

If the SDR works properly, when a meteor flows, a high-pitched sound of about 1000Hz will be output. When the reflection of radio waves is long, the sound may continue to be heard for 30 seconds or more.

↓ Like this. Be careful of the volume! https://soundcloud.com/user-395229817/20191218-175747a

We will analyze this sound with 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


#Frequency to analyze(Hz)
FREQ_LOW = 800
FREQ_HIGH = 1200

#Maximum volume
VOL_LIMIT = 2800


#Detection threshold 0~255
THRESHOLD = 70

#Recording time at the time of detection(s)
RECORD_SECONDS = 30

#Recording data save destination
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.')

#Whether to put out a graph
dataplot = False
parser = argparse.ArgumentParser(description="data plot enable")
parser.add_argument("-p",
                    metavar="--p",
                    type=str2bool,
                    nargs="+",
                    help="data plot enable"
                    )


#Whether to save MP3 when a meteor is detected
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

    #Detection flag
    detect = False

    #Recording flag
    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:
            #5 white noise~10%Adjusted to 25 or less
            vol[idx] = v.round()

            # VOL_Cut with LIMIT
            if VOL_LIMIT <= vol[idx]:
                vol[idx] = VOL_LIMIT

            #Normalized at 255
            vol[idx] = int(round(vol[idx] / VOL_LIMIT * 255, 0))

            #Maximum value
            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()

        #Threshold comparison
        if THRESHOLD <= self.maxVol:
            self.detect = True
        else:
            self.detect = False

        #Start recording
        if self.isRecording == False and self.detect:
            self.isRecording = True
            self.startRecording()

        #Processing during recording
        if self.isRecording:
            self.frames.append(in_data)

            #Processing at the end of recording
            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")

        #MP3 conversion
        self.convertMp3()

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

        #WAV removal
        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()

As a processing flow

** 1. FFT every 0.5 seconds 2. Extract data from 800 to 1200Hz 3. Determine if the volume in the above range exceeds the threshold 4. If it exceeds, save it as wav for 30 seconds 5. Convert wav to mp3 **

With that feeling, an mp3 file will be created each time a shooting star flows.

Threshold adjustment

Display the spectrogram to adjust the meteor detection threshold.

#Spectrogram drawing
# p → True
python3 /home/pi/streamfft.py -p True -s False

Since the FFT is performed every 0.5 seconds, the frequency resolution is 2Hz. Since the analysis is performed for 400Hz from 800 to 1200Hz, 200 points of data can be obtained by one FFT. I'm sorry it's unfriendly that it's not labeled, but the X-axis is 800Hz on the left and 1200Hz on the right. The Y axis is a value normalized from 0 to 255.

In the above figure, the volume VOL_LIMIT is adjusted so that the noise floor is about 10 out of 255 with only white noise contained. After making this adjustment, if the detection threshold is exceeded, it is a meteor! It is a rough judgment.

↓ When a meteor radio wave comes in, it spikes as shown in the figure below.

Let's go to the observation!

#Observation command
#s → Set to True and save the data
python3 /home/pi/streamfft.py -p False -s True

2019 Geminids Meteor Group Observation Results

The observation started at 2 o'clock on 12/14 and ended at 10 o'clock on 12/16, so it took about 56 hours. It seems that the beacon wave in Fukui stopped for about 12 hours, and I couldn't get any data during that time.

There were about 400 pieces of data, and apparently another noise-like thing overlapped or I played something that was not a meteor echo, and the impression was about 200 pieces. I couldn't get a good feeling because the AGC condition of SDR and the fixed threshold judgment did not mesh very well, or there was a data loss time.

This is the data of shooting star observation stations installed in various parts of Japan about two years ago. The value FFTed by Raspberry Pi is sent to AWS and the meteor echo is detected in real time on the server. Originally, the characteristic of the Gemini meteor shower is that the number of meteors gradually increases from about 12/12 and then drops sharply after the peak of 12 / 14-15.

Summary

The performance of RaspberryPi4 has improved significantly, and it is now possible to perform echo analysis while SDR. RaspberryPi3 is full of SDRs ...

In addition, this time I introduced it with a simple judgment of whether or not the threshold was exceeded, but when actually observing, various noises will come in, so I think it would be even more interesting to add advanced judgment processing to repel it.

ex. Doppler is applied to the echo reflected by the airplane. It is judged by the frequency displacement and repelled. Toka Toka