Reassembler module for satellite↔receiver distance and Doppler shift

[maehw]

Hi,

thanks for your great work in gr-iridium and iridium-toolkit. I appreciate the talks you held at various chaos events, got interested and finally got me some RX equipment. I have questions about reception quality, but I think I'll open another issue for that.

I've not dived too deep yet, but I've had some thoughts about Doppler shift for improved burst/frame decoding. I know that calculating Doppler shift here is a little late because it would actually have to be fed back to gr-iridium for the QPSK demodulator to consider for additional phase shifted symbol... it's probably better to use estimations from the TLE data, right? But I got interested in calculating it anyhow. And I'd be interested if you can see rx signal level correlate to satellite↔receiver distance.

The idea:

1. Get satellite position from ring alert frame (lat/long/alt; for altitudes > 100 km)
2. Calculate relative distance between (stationary) receiver and (moving) satellite vehicle
3. Calculate velocity of satellite vehicle as a change of relative distance over time
4. Calculate Doppler shift from velocity

If anyone finds this useful, this may become a PR. For now, it's the plain Python code (stored under iridium-toolkit/iridiumtk/reassembler/satdist.py on my machine):

#!/usr/bin/env python3
# vim: set ts=4 sw=4 tw=0 et pm=:
import math
import sys
import re
import os
from math import radians, sin, cos, acos, sqrt
from util import dt
from configparser import ConfigParser

from .base import *
from ..config import config, outfile

class SatDistance(Reassemble):
    def __init__(self):
        self.topic="IRA"
        self.last_dist = {}

        config = ConfigParser()
        config.read(['locations.ini', os.path.join(os.path.dirname(__file__), 'locations.ini')])

        location = "default"
        if location not in config:
            print("Location %s not defined" % location, file=sys.stderr)
            print("Available locations: ", ", ".join(config.sections()), file=sys.stderr)
            sys.exit(1)
        if 'lat' in config[location] and 'lon' in config[location]:
            self.rx_lat = config.getfloat(location, 'lat')
            self.rx_lon = config.getfloat(location, 'lon')
        else:
            print("Latitude and/or longitude not available")
            sys.exit(1)
    def filter(self,line):
        q=super().filter(line)
        if q==None: return None
        if q.typ=="IRA:":
            p=re.compile(r'sat:(\d+) beam:(\d+) (?:(?:aps|xyz)=\(([+-]?[0-9]+),([+-]?[0-9]+),([+-]?[0-9]+)\) )?pos=\(([+-][0-9.]+)/([+-][0-9.]+)\) alt=(-?[0-9]+) .* bc_sb:\d+(?: (.*))?')
            m=p.search(q.data)
            if(not m):
                print("Couldn't parse IRA: ",q.data, end=' ', file=sys.stderr)
            else:
                q.sat=  int(m.group(1))
                q.beam= int(m.group(2))
                if m.group(3) is not None:
                    q.xyz= [4*int(m.group(3)), 4*int(m.group(4)), 4*int(m.group(5))]
                q.lat=float(m.group(6))
                q.lon=float(m.group(7))
                q.alt=  int(m.group(8))
                if q.alt < 100: return None
                q.dist = self.get_distance(q)
                return q
    def process(self,q):
        q.enrich()
        strtime = dt.epoch(q.time).isoformat(timespec='centiseconds')
        v_light = 299792.458  # speed of light in vacuum; kilometers per second
        f_ra_nominal = 1.626270e9  # tx frequency of ring alert channel in Hz; taken from docs/channels.svg
        str = "%s %03d %s %6.2f %6.2f %3d %4d" % (strtime, q.sat, q.frequency, q.lat, q.lon, q.alt, int(q.dist))
        if q.sat in self.last_dist:
            delta_time = q.time - self.last_dist[q.sat]['time']
            delta_dist = q.dist - self.last_dist[q.sat]['dist']
            if delta_time > 1.0 and abs(delta_dist) >= 1.0:
                velocity = delta_dist/delta_time  # kilometers per second
                rel_doppler_shift = -velocity/v_light
                delta_f_pred = rel_doppler_shift * f_ra_nominal  # predicted absolute frequency shift in Hz
                delta_f_real = f_ra_nominal - q.frequency  # observed absolute frequency shift in Hz
                #str += " %.2f %.2f %.2f %s %.0f %.0f" % (delta_time, delta_dist, velocity, rel_doppler_shift, delta_f_pred, delta_f_real)
                str += " %.0f %.0f" % (delta_f_pred, delta_f_real)
                self.last_dist[q.sat] = {'time': q.time, 'dist': q.dist}
        else:
            self.last_dist[q.sat] = {'time': q.time, 'dist': q.dist}
        return [str]
    def consume(self,q):
        print(q, file=outfile)
    def get_distance(self,q):
        # assumptions for calculation: earth is a sphere not an ellipsoid;
        # it currently also does not take into account the receiver's altitude
        sat_lat = radians(q.lat)
        sat_lon = radians(q.lon)
        rx_lat = radians(self.rx_lat)
        rx_lon = radians(self.rx_lon)
        radius_earth = 6371  # mean earth radius in kilometers
        # calculate central angle between rx position and satellite position projected on earth
        central_angle = acos(sin(sat_lat)*sin(rx_lat) + cos(sat_lat)*cos(rx_lat)*cos(sat_lon-rx_lon))
        # calculate length of a virtual tunnel between rx position and satellite position projected on earth (cathetus)
        tunnel_length = 2 * radius_earth * sin(central_angle/2)
        # finally calculate the distance using Pythagorean theorem from the two known cathethi
        dist = sqrt(q.alt ** 2 + tunnel_length ** 2)
        return dist

modes=[
["satdist",       SatDistance,  ],
]

It can then be called python3 reassembler.py -i <something.parsed> -m satdist. I've added a pipe like | grep " 048 " to allow to filter for satellite no 48.

I first calculate the central angle between the rx position and satellite position projected on earth's surface ("spherical law of cosines", https://en.wikipedia.org/wiki/Great-circle_distance#Formulae). The rx position is read from the default section of 'locations.ini' (as this is already used by some other module), the satellite vehicle position is taken from the IRA frames.

When you use the commented line str += " %.2f %.2f %.2f %s %.0f %.0f" % (delta_time, delta_dist, velocity, rel_doppler_shift, delta_f_pred, delta_f_real) instead of str += " %.0f %.0f" % (delta_f_pred, delta_f_real), you'lle see some more values.

The command line output shows the estimated Doppler shift in Hz as calculated by my reassembler module and the deviation of the IRA rx frequency from the nominal IRA frequency of 1.626270 GHz (as per gr-iridium docs/channels.svg). TL;DR: Saldly, the numbers do not match and I don't really know why. I've only some ideas.

I'd expect the Doppler shift to become lower as long as the satellite vehice comes closer to the rx position and then to become higher again -- hence, the numbers to change in that way as well as the sign.

• Is the nominal ring alert tx channel frequency for every satellite at 1.626270 GHz as shown in the Iridium channels graphic? If this is the case, I see rx frequencies in the parsed bits of IRA frames deviate up to 56 kHz which exceeds the 36 kHz mentioned here.
• Is the position of the SV sent by itself not precise enough? (I've seen multiple timestamps for the same lat/long value being sent.)
• Are other numbers not precise enough? E.g. latitude/longitude being rounded to 2 decimal places? Ignoring altitude for the rx position?
• Is the approximation of earth beaing a sphere an not an ellipsoid to coarse?
• How "good" is the SV altitude? I e.g. see values jumping between 791..799 km -- or is this due to the ellipsoid?

Cheers

[Sec42]

A couple of notes. In our first iridium talk we calculated the estimated doppler shift from TLEs and compared them to the received frequencies in this video at ~11:00 - It seems we never published that code though.

Calculating from the SV positions is a neat idea, although made difficult by the fact that the location is sent with 4km granularity in x/y/z, so relatively imprecise. This means if you attempt to calculate the apparent speed of the satellite by two positions too close together, you will get wildly inaccurate results. You could attempt to compensate for this by taking the positions with a longer time interval between them, but that will of course introduce other inaccuracies.

Another alternative is to take only one position and calculate the apparent speed from the real speed and direction of the SV similar to what is done in beam-reception-plotter.py. (This code is not concerned with speed, but you should be able to derive the speed from the SV height)

For some of your questions:
The maximum doppler shift of 36kHz mentioned in the source is in both directions of the center frequency, for a total of 72kHz.
The altitude you see is different mostly depending on which plane a give satellite is on (and I think the orbits are also slightly elliptical)

[maehw]

Thanks for sharing your thoughts.

I might plot the calculated Doppler shift next to the IRA frequency - this should allow to visually inspect the differences.

I think it was a good way to initially get started with the Iridium toolkit, especially with reassembling and IRA frames.

beam-reception-plotter.py seems very advanced and currently is above my knowledge/understanding.

As this has all been more or less playing around... the overall question here is: does gr-iridium already take Doppler shift into account yet (I think you mentioned it in an "old talk" that it's not the case, but this may have changed). Do you see potention for improving demodulation? Or are bursts "short enough" anyway so that Doppler shift can be ignored because they do not lead to frequency drift/ phase offset in a relevant amount?