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

IMG_2481_.jpg ・ Raspberry Pi4 ・ USB audio conversion cable ・ Audio cable ・ AirspyMini (Software Radio) ・ 50MHz band antennaCoaxial 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.

IMG_2474.jpg 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. Screenshot of プレビュー (2019-12-16 22-37-38).png

・ 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.

Screenshot of プレビュー (2019-12-17 0-55-17).png

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

Desktop.png

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. spike

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

Screenshot of プレビュー (2019-12-17 11-10-22).png

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. Screenshot of Google Chrome (2019-12-17 11-34-13).png

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

Recommended Posts

Observe the Geminids meteor shower with Raspberry Pi 4
GPGPU with Raspberry Pi
DigitalSignage with Raspberry Pi
Using the digital illuminance sensor TSL2561 with Raspberry Pi
Display images taken with the Raspberry Pi camera module
Try to visualize the room with Raspberry Pi, part 1
Take the value of SwitchBot thermo-hygrometer with Raspberry Pi
Log the value of SwitchBot thermo-hygrometer with Raspberry Pi
Mutter plants with Raspberry Pi
Play with the Raspberry Pi Zero WH camera module Part 1
I tried using the DS18B20 temperature sensor with Raspberry Pi
[Raspberry Pi] Stepping motor control with Raspberry Pi
Use vl53l0x with Raspberry Pi (python)
Serial communication with Raspberry Pi + PySerial
OS setup with Raspberry Pi Imager
VPN server construction with Raspberry Pi
Try moving 3 servos with Raspberry Pi
Using a webcam with Raspberry Pi
Control the motor with a motor driver using python on Raspberry Pi 3!
I tried to automate the watering of the planter with Raspberry Pi
I learned how the infrared remote control works with Raspberry Pi
Periodically log the value of Omron environment sensor with Raspberry Pi
Measure SIM signal strength with Raspberry Pi
Pet monitoring with Rekognition and Raspberry pi
A memo to simply use the illuminance sensor TSL2561 with Raspberry Pi 2
Ask for Pi with the bc command
Hello World with Raspberry Pi + Minecraft Pi Edition
Build a Tensorflow environment with Raspberry Pi [2020]
Use python on Raspberry Pi 3 to light the LED with switch control!
Get BITCOIN LTP information with Raspberry PI
Try fishing for smelt with Raspberry Pi
Programming normally with Node-RED programming on Raspberry Pi 3
Use the Grove sensor on the Raspberry Pi
Improved motion sensor made with Raspberry Pi
Try Object detection with Raspberry Pi 4 + Coral
Power SG-90 servo motor with raspberry pi
Working with sensors on Mathematica on Raspberry Pi
Use PIR motion sensor with raspberry Pi
Make a wash-drying timer with a Raspberry Pi
Infer Custom Vision model with Raspberry Pi
Operate an oscilloscope with a Raspberry Pi
Logging the value of Omron environment sensor with Raspberry Pi (USB type)
Create a car meter with raspberry pi
Inkbird IBS-TH1 value logged with Raspberry Pi
Working with GPS on Raspberry Pi 3 Python
Make a thermometer with Raspberry Pi and make it visible on the browser Part 3
I tweeted the illuminance of the room with Raspberry Pi, Arduino and optical sensor
Periodically notify the processing status of Raspberry Pi with python → Google Spreadsheet → LINE
Try to use up the Raspberry Pi 2's 4-core CPU with Parallel Python
Discord bot with python raspberry pi zero with [Notes]
Media programming with Raspberry Pi (preparation for audio)
I tried L-Chika with Raspberry Pi 4 (Python edition)
MQTT RC car with Arduino and Raspberry Pi
Power on / off your PC with raspberry pi
Use Majoca Iris elongated LCD with Raspberry Pi
CSV output of pulse data with Raspberry Pi (CSV output)
Get CPU information of Raspberry Pi with Python
Sound the buzzer using python on Raspberry Pi 3!
Play with your Ubuntu desktop on your Raspberry Pi 4
Get temperature and humidity with DHT11 and Raspberry Pi
Stock investment analysis app made with Raspberry Pi