Massive update. NOAA should work.

.gitignore

@@ -1 +1,5 @@
+
+__pycache__/
+.vscode/
 config.ini
+.gitignore

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)

config.ini

@@ -1,9 +1,17 @@
 [RADIO]
-HOSTNAME = your.host.name
+HOSTNAME = office.toyaga.eu
-PORT = 1234 
+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

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)

src/classes/NOAADecoder.py

@@ -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

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
+
+        

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"]:
+    if satellite["passing_soon"] or satellite["passing_now"]:
-        info(f"Passing soon: {satellite['name']}")
+        if not satellite["passing_now"]:
-        sleeper = satellite["passing_soon"]
+            info(f"Passing soon: {satellite['name']}")
-        while sleeper > 0:
+            sleeper = satellite["passing_soon"]
-            if sleeper % 100 == 0:
+            while sleeper > 0:
-                warning(f"Passing soon: {satellite['name']}, {sleeper} seconds left")
+                if sleeper % 100 == 0:
-            elif sleeper <= 10:
+                    warning(f"Passing soon: {satellite['name']}, {sleeper} seconds left")
-                warning(f"{sleeper}:{satellite['name']}")
+                elif sleeper <= 10:
-            sleeper -= 1
+                    warning(f"{sleeper}:{satellite['name']}")
-            sleep(1)
+                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()
         
@@ -101,11 +125,15 @@ def main():
         # When everything is finished, sleep for 5 seconds and restart      
         sleep(5)
          
-        ## TODO: 
+    #     ## TODO: 
-        # Record audio
+    #     # Record audio
-        # Generate image
+    #     # Generate image
-        # Upload to discord
+    #     # Upload to discord
-        # Glue them together
+    #     # Glue them together
 
 if __name__ == "__main__":
-    main()
+    try:
+        main()
+    except KeyboardInterrupt:
+        print("Exiting...")
+        sys.exit(0)