From dd4c71fa72007568e59b303e56f996c5931ea376 Mon Sep 17 00:00:00 2001 From: Boyan Date: Fri, 30 Jun 2023 22:11:34 +0300 Subject: [PATCH] Massive update. NOAA should work. --- .gitignore | 4 ++ README.md | 4 +- config.ini | 14 ++++- src/classes/APT.py | 107 +++++++++++++++++++++++++++++++++++ src/classes/NOAADecoder.py | 65 --------------------- src/classes/RtlTcpClient.py | 109 ++++++++++++++++++++++++++++++++++++ src/main.py | 74 ++++++++++++++++-------- 7 files changed, 284 insertions(+), 93 deletions(-) create mode 100644 src/classes/APT.py delete mode 100644 src/classes/NOAADecoder.py create mode 100644 src/classes/RtlTcpClient.py diff --git a/.gitignore b/.gitignore index 2fa7ce7..c247281 100644 --- a/.gitignore +++ b/.gitignore @@ -1 +1,5 @@ + +__pycache__/ +.vscode/ config.ini +.gitignore diff --git a/README.md b/README.md index a3f5032..07afec2 100644 --- a/README.md +++ b/README.md @@ -1,7 +1,7 @@ # RTL-TCP automatic weather satellite image -This project is left for dead until I find a suitable library which can interact with rtl_tcp. + ## Capture, decode and send -Using ~~`pyrtlsdr`~~, interact with an RTL-TCP server. Tell it when and where to listen for weather satellites. Record said satellites. Locally decode the signals into images. Post the images to Discord. +Using [RtlTcpClient](src/classes/RtlTcpClient.py), interact with an RTL-TCP server. Tell it when and where to listen for weather satellites. Record said satellites. Locally decode the signals into images. Post the images to Discord. ### Ideas(difficulty level) * Clearness % algorithm(5/5) diff --git a/config.ini b/config.ini index fe677b1..0cfbe58 100644 --- a/config.ini +++ b/config.ini @@ -1,9 +1,17 @@ [RADIO] -HOSTNAME = your.host.name -PORT = 1234 +HOSTNAME = office.toyaga.eu +FREQUENCY = 137.1 +FREQCORRECTION = 0 +SAMPLERATE = 208000 +PORT = 51234 + +[LOCATION] +LON = 23.333300 +LAT = 42.682120 +ALT = 550 [DISCORD] -TOKEN = your.discord.token +WEBHOOK_URL = https://discordapp.com/api/webhooks/1234567890/abcdefghijklmnopqrstuvwxyz [LOG] LOGFILE = radio.log diff --git a/src/classes/APT.py b/src/classes/APT.py new file mode 100644 index 0000000..a657105 --- /dev/null +++ b/src/classes/APT.py @@ -0,0 +1,107 @@ +# Taken from: https://github.com/zacstewart/apt-decoder/ +# Slightly modified + +import numpy +import scipy.io.wavfile +import scipy.signal +import sys +from PIL import Image + +class APT(object): + + RATE = 20800 + NOAA_LINE_LENGTH = 2080 + + def __init__(self, signal) -> None: + + # def __init__(self, filename): + # (rate, self.signal) = scipy.io.wavfile.read(filename) + # if rate != self.RATE: + # raise Exception("Resample audio file to {}".format(self.RATE)) + + + # Keep only one channel if audio is stereo + # if self.signal.ndim > 1: + # self.signal = self.signal[:, 0] + + self.signal = signal + + truncate = self.RATE * int(len(self.signal) // self.RATE) + self.signal = self.signal[:truncate] + print(self.signal) + + def decode(self, outfile=None): + hilbert = scipy.signal.hilbert(self.signal) + filtered = scipy.signal.medfilt(numpy.abs(hilbert), 5) + reshaped = filtered.reshape(len(filtered) // 5, 5) + digitized = self._digitize(reshaped[:, 2]) + matrix = self._reshape(digitized) + image = Image.fromarray(matrix) + if not outfile is None: + image.save(outfile) + image.show() + return matrix + + def _digitize(self, signal, plow=0.5, phigh=99.5): + ''' + Convert signal to numbers between 0 and 255. + ''' + (low, high) = numpy.percentile(signal, (plow, phigh)) + delta = high - low + data = numpy.round(255 * (signal - low) / delta) + data[data < 0] = 0 + data[data > 255] = 255 + return data.astype(numpy.uint8) + + def _reshape(self, signal): + ''' + Find sync frames and reshape the 1D signal into a 2D image. + + Finds the sync A frame by looking at the maximum values of the cross + correlation between the signal and a hardcoded sync A frame. + + The expected distance between sync A frames is 2080 samples, but with + small variations because of Doppler effect. + ''' + # sync frame to find: seven impulses and some black pixels (some lines + # have something like 8 black pixels and then white ones) + syncA = [0, 128, 255, 128]*7 + [0]*7 + + # list of maximum correlations found: (index, value) + peaks = [(0, 0)] + + # minimum distance between peaks + mindistance = 2000 + + # need to shift the values down to get meaningful correlation values + signalshifted = [x-128 for x in signal] + syncA = [x-128 for x in syncA] + for i in range(len(signal)-len(syncA)): + corr = numpy.dot(syncA, signalshifted[i : i+len(syncA)]) + + # if previous peak is too far, keep it and add this value to the + # list as a new peak + if i - peaks[-1][0] > mindistance: + peaks.append((i, corr)) + + # else if this value is bigger than the previous maximum, set this + # one + elif corr > peaks[-1][1]: + peaks[-1] = (i, corr) + + # create image matrix starting each line on the peaks found + matrix = [] + for i in range(len(peaks) - 1): + matrix.append(signal[peaks[i][0] : peaks[i][0] + 2080]) + + return numpy.array(matrix) + + +if __name__ == '__main__': + apt = APT(sys.argv[1]) + + if len(sys.argv) > 2: + outfile = sys.argv[2] + else: + outfile = None + apt.decode(outfile) diff --git a/src/classes/NOAADecoder.py b/src/classes/NOAADecoder.py deleted file mode 100644 index acfb32e..0000000 --- a/src/classes/NOAADecoder.py +++ /dev/null @@ -1,65 +0,0 @@ -# Based on this article: -# https://medium.com/swlh/decoding-noaa-satellite-images-using-50-lines-of-code-3c5d1d0a08da - -import scipy.io.wavfile as wav -import scipy.signal as signal -import numpy as np -import matplotlib.pyplot as plt -from PIL import Image -from tqdm import tqdm - -class NOAADecoder: - def __init__(self, iq_data=None, filename=None) -> None: - # Check if iq_data or filename are set correctly - if iq_data and filename: - raise ValueError("Both iq_data and filename are set.") - if not iq_data and not filename: - raise ValueError("Neither iq_data nor filename are set.") - - # Set iq_data or filename - if filename: - self.filename=filename - self.sample_rate, self.data = wav.read(f'{self.filename}') - if iq_data: - self.sample_rate = 2048000 - self.data = self.iq_data - - def hilbert(self, data): - # Hilbert transform - analytical_signal = signal.hilbert(data) - amplitude_envelope = np.abs(analytical_signal) - return amplitude_envelope - - def decode(self, resample=4): - - # Resample data - data = self.data[::resample] - sample_rate = self.sample_rate//resample - - # Demodulate - data_am = self.hilbert(data) - frame_width = int(0.5*sample_rate) - - # Image processing - w, h = frame_width, data_am.shape[0]//frame_width - image = Image.new('L', (w, h)) - px, py = 0, 0 - - for p in tqdm(range(data_am.shape[0])): - if data_am[p][0] != data_am[p][1]: - raise BufferError("Mismatch in data array.") - lum = int(data_am[p][0]//32 - 32) - if lum < 0: lum = 0 - if lum > 255: lum = 255 - image.putpixel((px, py),lum) - px += 1 - if px >= w: - px = 0 - py += 1 - if py >= h: - break - image = image.resize((w, 4*h)) - plt.imshow(image, cmap='gray') - # Save image - image.save('noaa.png') - return \ No newline at end of file diff --git a/src/classes/RtlTcpClient.py b/src/classes/RtlTcpClient.py new file mode 100644 index 0000000..802b50f --- /dev/null +++ b/src/classes/RtlTcpClient.py @@ -0,0 +1,109 @@ +#!/usr/bin/env python3 + + +import socket +import struct +import time +import numpy as np +from scipy import signal +from scipy.io import wavfile +import sys +from math import pi + +SET_FREQUENCY = 0x01 +SET_SAMPLERATE = 0x02 +SET_GAINMODE = 0x03 +SET_GAIN = 0x04 +SET_FREQENCYCORRECTION = 0x05 + +class RtlTCPClient: + def __init__(self): + # Connect to a remote rtl_tcp server + # TODO: Do it with config file + self.remote_host = "office.toyaga.eu" + self.remote_port = 51234 + self.conn = socket.socket(socket.AF_INET, socket.SOCK_STREAM) + self.conn.connect((self.remote_host, self.remote_port)) + self.send_command(SET_SAMPLERATE, 20800) + + def tune(self, freq): + # Tune to frequency + self.send_command(SET_FREQUENCY, int(freq)) + + def send_command(self, command, parameter): + # Send command to server + cmd = struct.pack(">BI", command, parameter) + self.conn.send(cmd) + + def receive_data(self, duration): + # Receive data from server + raw_data = b"" + data_amount = 0 + second = time.time() + deadline = second + duration + while True: + + try: + if time.time() - second > 1: + second = time.time() + data_amount = 0 + data = self.conn.recv(1024) + data_amount += len(data) + if not data: + raise ConnectionError("Connection closed by server") + if time.time() >= deadline: + break + + raw_data += data + except socket.timeout: + break + + # Close the socket connection + self.conn.close() + + # Turn the raw data into a numpy array + iq = np.frombuffer(raw_data, dtype=np.uint16) + + # Unpack a sequence of bytes to a sequence of normalized complex numbers + # iq = (iq - 127.5) / 128 + # iq = iq.view(np.complex64) + # # Remove the leading 12 bytes that contain the header + iq = iq[12:] + return iq + def demodulate_FM(self, iq): + angles = np.angle(iq) + + # Determine phase rotation between samples + # (Output one element less, that's we always save last sample + # in remaining_data) + rotations = np.ediff1d(angles) + + # Wrap rotations >= +/-180ยบ + rotations = (rotations + np.pi) % (2 * np.pi) - np.pi + + # Convert rotations to baseband signal + output_raw = np.multiply(rotations, 0.99 / (pi * 200000 / (20800 / 2))) + output_raw = np.clip(output_raw, -0.999, +0.999) + + # Scale to signed 16-bit int + output_raw = np.multiply(output_raw, 32767) + output_raw = output_raw.astype(int) + + # Missing: low-pass filter and deemphasis filter + # (result may be noisy) + + # Output as raw 16-bit, 1 channel audio + bits = struct.pack(('<%dh' % len(output_raw)), *output_raw) + + # Write to raw file + with open('fm.raw', 'wb') as f: + f.write(bits) + + # return iq + + + def demodulate_AM(self, iq): + # https://stackoverflow.com/questions/61106688/rtl-sdr-iq-am-demodulation + pass + + diff --git a/src/main.py b/src/main.py index 7f3d2f5..1402bc5 100644 --- a/src/main.py +++ b/src/main.py @@ -2,10 +2,14 @@ from rtlsdr.rtlsdrtcp import RtlSdrTcpClient from configparser import ConfigParser import predict from json import loads -from logging import getLogger, basicConfig, INFO, info, warning +from logging import getLogger, basicConfig, INFO, info, DEBUG, debug, warning from datetime import datetime from multiprocessing.pool import Pool from time import sleep +import sys + +from classes.RtlTcpClient import RtlTCPClient +from classes.APT import APT basicConfig(level=INFO) logger = getLogger(__name__) @@ -33,7 +37,8 @@ def predict_satellite_passes(): "start": datetime.now().fromtimestamp(transit.start).isoformat(), "duration": round(transit.duration()/60, 2) , "max_elevation": round(transit.peak()["elevation"], 2), - "passing_soon": delta.seconds if delta.seconds < 10*60 else False + "passing_soon": delta.seconds if delta.seconds < 10*60 else False, + "passing_now": delta.seconds if delta.seconds <= 0 else False }) @@ -48,34 +53,45 @@ def generate_image(): pass -def record_audio(): +def setup(): # Connect to a remote rtl_tcp server sdr = RtlSdrTcpClient(hostname=config.get("RADIO", "HOSTNAME"), port=config.getint("RADIO", "PORT")) # Configure the SDR sdr.set_sample_rate(config.getint("RADIO", "SAMPLERATE")) - #TODO: Add bandwidth - sdr.set_center_freq(config.getfloat("RADIO", "FREQUENCY")*10**6) sdr.set_gain('auto') # Change if you know what you're doing + #TODO: Add bandwidth def satellite_tracking(satellite): # Brings all the functions together - if satellite["passing_soon"]: - info(f"Passing soon: {satellite['name']}") - sleeper = satellite["passing_soon"] - while sleeper > 0: - if sleeper % 100 == 0: - warning(f"Passing soon: {satellite['name']}, {sleeper} seconds left") - elif sleeper <= 10: - warning(f"{sleeper}:{satellite['name']}") - sleeper -= 1 - sleep(1) + if satellite["passing_soon"] or satellite["passing_now"]: + if not satellite["passing_now"]: + info(f"Passing soon: {satellite['name']}") + sleeper = satellite["passing_soon"] + while sleeper > 0: + if sleeper % 100 == 0: + warning(f"Passing soon: {satellite['name']}, {sleeper} seconds left") + elif sleeper <= 10: + warning(f"{sleeper}:{satellite['name']}") + sleeper -= 1 + sleep(1) + info(f"Passing now: {satellite['name']}") + + if "NOAA" in satellite["name"]: + tcp = RtlTCPClient() + tcp.tune(satellite["freq"]) + data = tcp.receive_data(satellite["duration"]*60) + + apt = APT(data) + apt.decode(f"{satellite['name']}_{datetime.now().isoformat()}.png") # Everything else goes here + else: sleeper = datetime.fromisoformat(satellite["start"]).timestamp() - datetime.now().timestamp() info(f"{satellite['name']} is not passing soon, sleeping until then :) (sleeping for {str(round(sleeper/(60*60),1)) + ' hours' if sleeper>60*60 else str(round(sleeper/60)) + ' minutes'})") + debug(satellite) sleep(sleeper) # Example of parallel running function @@ -87,9 +103,17 @@ def satellite_tracking(satellite): # print("Hey, I can run right here, forever!") # sleep(1) -def main(): +def debugging(freq:float, duration:float, name:str): + tcp = RtlTCPClient() + tcp.tune(freq) + data = tcp.receive_data(duration) + print("Done receiving") + apt = APT(data) + apt.decode(f"{name}_{datetime.now().isoformat()}.png") +def main(): # Predict passes and get ready to loop + # debugging(137.1e6, 10, "NOAA-19") while True: passes = predict_satellite_passes() @@ -100,12 +124,16 @@ def main(): # When everything is finished, sleep for 5 seconds and restart sleep(5) - - ## TODO: - # Record audio - # Generate image - # Upload to discord - # Glue them together + + # ## TODO: + # # Record audio + # # Generate image + # # Upload to discord + # # Glue them together if __name__ == "__main__": - main() + try: + main() + except KeyboardInterrupt: + print("Exiting...") + sys.exit(0) \ No newline at end of file