FMTlab 2026-02-13¶

# Data from monitoring RWM:
TOTAL_DELAY_FMTLAB = 28.5e-3 # 28.5 ms ± 1.2 ms or so
# TRX frequency is too high by about -0.01402 ppm 
FREQ_DEVIATION = -0.014 # PPM
FILENAME = "2026-03-06_205119"
DIAL_FREQUENCY = 3563400
RF_FREQUENCY = DIAL_FREQUENCY * (1 + FREQ_DEVIATION * 1e-6)
RF_FREQUENCY

grab = GrabSound(FILENAME)
rx_qrg = RF_FREQUENCY
print(f"RX oscillator assumed to be running at {rx_qrg:.3f} Hz")

def plot_time_deviation(grab) -> pandas.DataFrame:
    deltas = [(grab.time_raw(tp.written_samples) - tp.t_ns)*1e-9 for tp in grab.timed_points[:-1]]
    deltas_without_outliers = [d if abs(d) < 0.5e-3 else None for d in deltas]
    outlier_count = 0
    for d in deltas_without_outliers:
        if d is None:
            outlier_count += 1
    
    df_delta = pandas.DataFrame(
        {
            "delta": deltas,
            f"delta without {outlier_count} outliers": deltas_without_outliers,
        },
        index = [tp.written_samples for tp in grab.timed_points[:-1]]
    )
    df_delta.plot(
        subplots=True,
        figsize=(12,6),
        grid=True,
        title="Precision of Sound Card",
        # legend=False,
        xlabel=f"Samples at {grab.samples_per_second} samples/second.",
        ylabel="Difference / s"
        # ylabel="How much ahead (+) or delayed (-) of average flow\ndid individual sample packages arrive, in seconds."
    )
    return df_delta
plot_time_deviation(grab)

def plot_smooth(grab) -> pandas.DataFrame:
    df_s_delta = pandas.DataFrame(
        {
            "smooth delta": [tpe.backprojected_sample_start*1e-9 for tpe in grab.timed_points_extended]
        },
        index=[tpe.written_samples for tpe in grab.timed_points_extended]
    )
    df_s_delta.plot(
        figsize=(12,6),
        grid=True,
        title="Precision of Sound Card",
        ylabel="Deviation / s",
        xlabel=f"Samples at {grab.samples_per_second} samples/second.",        
    )
    return df_s_delta

plot_smooth(grab), len(grab.timed_points)

def plot_delta_to_next(grab) -> pandas.DataFrame:
    pns = list(zip(grab.timed_points_extended[:-1],grab.timed_points_extended[1:]))
    
    df_delta_to_next = pandas.DataFrame(
        {
            "from_one_to_next": [
                (pn[1].backprojected_sample_start-pn[0].backprojected_sample_start)*1e-9 for pn in pns
            ]
        },
        index=[pn[0].written_samples for pn in pns]
    )
    df_delta_to_next.plot(
        figsize=(12,6),
        grid=True,
        title="From one smoothing point to next, time difference of when the two think recording started.",
        ylabel="Difference / s",
        xlabel=f"Samples at {grab.samples_per_second} samples/second.",        
    )
    return df_delta_to_next

df_delta_to_next = plot_delta_to_next(grab)

# Hunt the problem

df_delta_to_next[0.0001 < df_delta_to_next["from_one_to_next"]]

grab.datetime(26665984), grab.datetime(26683392)

def sanity_ns_to_index_and_back(grab):
    start_ns = grab.time_ns(0)
    end_ns = grab.time_ns(grab.num_of_samples -1)
    start_ns, end_ns, type(start_ns), type(end_ns)
    ns_and_back_delta = []
    t_nss = []
    # Cannot use numpy.linspace, as internally it calculates with floats
    # and then, rounding errors cause the initial value to be too low.
    delta_per_step = (end_ns - start_ns) // 2999
    t_ns = start_ns
    while t_ns < end_ns:
        indx = grab.sample_index_from_ns(t_ns)
        t_ns_back = grab.time_ns(indx)
        ns_and_back_delta.append((t_ns_back-t_ns)*1e-9)
        t_nss.append(t_ns)
        t_ns += delta_per_step
    
    df = pandas.DataFrame({
        "ns_and_back_delta": ns_and_back_delta
        },
        index = t_nss
    )
    df.plot(
        title="Going from ns to index and back to ns, resulting delta.",
        grid=True,
        figsize=(12,6)
    )
    return df
print(f"Compare with {1/grab.samples_per_second} s per sample")
sanity_ns_to_index_and_back(grab)

def sanity_index_to_ns_and_back(grab):
    indexs = numpy.linspace(0, grab.num_of_samples, 2999, endpoint=False, dtype=numpy.int64)
    for indx in indexs:
        t_ns = grab.time_ns(indx)
        indx_back = grab.sample_index_from_ns(t_ns)
        if indx != indx_back:
            raise RuntimeError(f"{indx} → {t_ns} → {indx_back}")
    for indx in indexs:
        timestamp = grab.datetime(indx)
        indx_back = grab.sample_index(timestamp)
        if indx != indx_back:
            raise RuntimeError(f"{indx} → {timestamp} → {indx_back}")
    print("All is well in nanoseconds.")

def sanity_index_to_datetime_and_back(grab):
    indexs = numpy.linspace(0, grab.num_of_samples, 499, endpoint=False, dtype=numpy.int64)
    for indx in indexs:
        time = grab.datetime(indx)
        indx_back = grab.sample_index(time)
        if indx != indx_back:
            raise RuntimeError("f{indx} → {time} → {indx_back}")
    print("All is well in timestamps")
    
sanity_index_to_ns_and_back(grab)
sanity_index_to_datetime_and_back(grab)

start_time = grab.start_time()
end_time = grab.end_time()
signal_start_time = start_time.replace(minute=58, second=0, microsecond=0)
signal_hour_start = signal_start_time.replace(hour=signal_start_time.hour + 1, minute=0, second=0, microsecond=0)
print(f"Recording from {start_time} to {end_time},\n"
    f"signal expected to start not before {signal_start_time}\n"
    f"or actually at about {signal_hour_start}"
)

import heapq

# The graph has more information if vmax is used to remove the few brightest values.
def find_reasonably_large_amplitude(fraction:float, fft_amplitudes):
    all_negative_amps = [-fft_amplitudes[f_index][t_index]
                for t_index in range(0,len(fft_amplitudes[0]))
                for f_index in range(0,len(fft_amplitudes))
    ]
    heapq.heapify(all_negative_amps)
    for disregard_count in range(0,int(len(all_negative_amps)*fraction)):
        heapq.heappop(all_negative_amps)
    return -all_negative_amps[0]

vmax = find_reasonably_large_amplitude(0.01, fft_amplitudes)
vmax

def plot(fft_amplitudes, fft_engine, from_t,
         vmax = vmax,
         zoom_t_start = None,
         zoom_t_end = None,
         zoom_f_start = 0,
         zoom_f_end = 3000
):
    delta_t_per_x = fft_engine.delta_t   
    if zoom_t_start is None:
        zoom_t_start = from_t
    if zoom_t_end is None:
        zoom_t_end = from_t + timedelta(seconds=(1+len(fft_amplitudes[0]))*delta_t_per_x)

    time_index_and_time = []
    for time_index in range(0, len(fft_amplitudes[0])):
        time = from_t + timedelta(seconds=time_index*delta_t_per_x)
        if zoom_t_start <= time < zoom_t_end:
            time_index_and_time.append((time_index, time))

    delta_f_per_y = fft_engine.delta_f
    freq_index_and_freq = []
    for freq_index in range(0, len(fft_amplitudes)):
        freq = freq_index * delta_f_per_y
        if zoom_f_start <= freq < zoom_f_end:
            freq_index_and_freq.append((freq_index, freq))

    fft_amplitudes_filtered = [
        [
            fft_amplitudes[freq_index][time_index]
            for time_index, time in time_index_and_time
        ]
        for freq_index, freq in freq_index_and_freq
    ]
    
    fig, axes = plt.subplots(figsize=(12, 8))
    cmap = plt.colormaps["plasma"]
    x_values = [(time - signal_hour_start).total_seconds()/60 for time_index, time in time_index_and_time]
    y_values = [freq for freq_index, freq in freq_index_and_freq]
    axes.pcolormesh(x_values, y_values, fft_amplitudes_filtered,
                    vmax=vmax, cmap=cmap, norm="linear")
    axes.set_xlabel("time into full hour / minutes")
    axes.set_ylabel("Frequency / Hz")
plot(fft_amplitudes, fft_engine, signal_start_time)

LOWER_AF = 2300
HIGHER_AF = 2325
ZOOM_VMAX = 4e7
# Finding first carrier
plot(fft_amplitudes, fft_engine, signal_start_time, ZOOM_VMAX,
     signal_hour_start.replace(minute=4, second=0, microsecond=0),
     signal_hour_start.replace(minute=5, second=10, microsecond=0),
     LOWER_AF, HIGHER_AF)
print(f"{fft_engine.delta_f:.5f} Hz per bin")

audio_frequencies = []
into_the_carrier_minute_timestamps = []

# Finding first audio frequency

def analyze_for_frequencies(grab: GrabSound, carrier_detect_from, carrier_detect_to, plot=False):
    frequencies, densities = scipy.signal.periodogram(
        grab.samples[grab.sample_index(carrier_detect_from):grab.sample_index(carrier_detect_to)],
        fs=grab.samples_per_second
    )    
    af_measured, density = functools.reduce(lambda b,n: b if n[1] < b[1] else n,
                                            zip(frequencies, densities))
    if plot:
        fds = [(f,d) for f,d in zip(frequencies, densities) if af_measured - 0.5 < f < af_measured + 0.5]
        df = pandas.DataFrame(
            {"relative Amplitude": [fd[1] for fd in fds]},
            index=[fd[0] for fd in fds]
        )
        df.plot(grid=True, xlabel="Frequency / Hz", figsize=(12,6))
    return frequencies, densities, af_measured

itcm = signal_hour_start.replace(minute=4, second=6)
into_the_carrier_minute_timestamps.append(itcm)

freq, densities, af_measured = \
    analyze_for_frequencies(grab, 
                            itcm,
                            itcm + timedelta(seconds=54), True)
print(f"Audio frequency measured in first minute: {af_measured:.2f}, "
      f"TX QRG: {rx_qrg + af_measured:.2f}")
audio_frequencies.append(af_measured)

# Finding second carrier
plot(fft_amplitudes, fft_engine, signal_start_time, ZOOM_VMAX,
     signal_hour_start.replace(minute=6, second=0, microsecond=0),
     signal_hour_start.replace(minute=7, second=30, microsecond=0),
     LOWER_AF, HIGHER_AF)

# Finding second audio frequency

into_the_carrier_minute_timestamps.append(signal_hour_start.replace(minute=6, second=21))

freq, densities, af_measured = \
    analyze_for_frequencies(grab,
                            into_the_carrier_minute_timestamps[-1],
                            into_the_carrier_minute_timestamps[-1] + timedelta(seconds=54),
                            True)
print(f"Audio frequency measured in minute after {into_the_carrier_minute_timestamps[-1]}: {af_measured:.2f}\n"
      f"TX QRG: {rx_qrg + af_measured:.2f}")
audio_frequencies.append(af_measured)

# Finding third carrier

plot(fft_amplitudes, fft_engine, signal_start_time, ZOOM_VMAX,
     signal_hour_start.replace(minute=8, second=30, microsecond=0),
     signal_hour_start.replace(minute=9, second=42, microsecond=0),
     LOWER_AF, HIGHER_AF)

# Finding third audio frequency

into_the_carrier_minute_timestamps.append(signal_hour_start.replace(minute=8, second=42))

freq, densities, af_measured = \
    analyze_for_frequencies(grab,
                            into_the_carrier_minute_timestamps[-1],
                            into_the_carrier_minute_timestamps[-1] + timedelta(seconds=54),
                            True)
print(f"Audio frequency measured in minute after {into_the_carrier_minute_timestamps[-1]}: {af_measured:.2f}\n"
      f"TX QRG: {rx_qrg + af_measured:.2f}")
audio_frequencies.append(af_measured)

# Finding forth carrier

plot(fft_amplitudes, fft_engine, signal_start_time, ZOOM_VMAX,
     signal_hour_start.replace(minute=10, second=48, microsecond=0),
     signal_hour_start.replace(minute=12, second=0, microsecond=0),
     LOWER_AF, HIGHER_AF)

# Finding forth audio frequency

into_the_carrier_minute_timestamps.append(signal_hour_start.replace(minute=10, second=57))

freq, densities, af_measured = \
    analyze_for_frequencies(grab,
                            into_the_carrier_minute_timestamps[-1],
                            into_the_carrier_minute_timestamps[-1] + timedelta(seconds=54),
                            True)
print(f"Audio frequency measured in minute after {into_the_carrier_minute_timestamps[-1]}: {af_measured:.2f}\n"
      f"TX QRG: {rx_qrg + af_measured:.2f}")
audio_frequencies.append(af_measured)

# Finding fifth carrier

plot(fft_amplitudes, fft_engine, signal_start_time, ZOOM_VMAX,
     signal_hour_start.replace(minute=13, second=12, microsecond=0),
     signal_hour_start.replace(minute=14, second=24, microsecond=0),
     LOWER_AF, HIGHER_AF)

# Finding fifth audio frequency
into_the_carrier_minute_timestamps.append(signal_hour_start.replace(minute=13, second=15))

freq, densities, af_measured = \
    analyze_for_frequencies(grab,
                            into_the_carrier_minute_timestamps[-1],
                            into_the_carrier_minute_timestamps[-1] + timedelta(seconds=54),
                            True)
print(f"Audio frequency measured in minute after {into_the_carrier_minute_timestamps[-1]}: {af_measured:.2f}\n"
      f"TX QRG: {rx_qrg + af_measured:.2f}")
audio_frequencies.append(af_measured)

# Finding sixth carrier

plot(fft_amplitudes, fft_engine, signal_start_time, ZOOM_VMAX,
     signal_hour_start.replace(minute=15, second=24, microsecond=0),
     signal_hour_start.replace(minute=16, second=42, microsecond=0),
     LOWER_AF, HIGHER_AF)

# Finding sixth audio frequency

into_the_carrier_minute_timestamps.append(signal_hour_start.replace(minute=15, second=33))

freq, densities, af_measured = \
    analyze_for_frequencies(grab,
                            into_the_carrier_minute_timestamps[-1],
                            into_the_carrier_minute_timestamps[-1] + timedelta(seconds=54),
                            True)
print(f"Audio frequency measured in minute after {into_the_carrier_minute_timestamps[-1]}: {af_measured:.2f}\n"
      f"TX QRG: {rx_qrg + af_measured:.2f}")
audio_frequencies.append(af_measured)

audio_frequencies, into_the_carrier_minute_timestamps

avg_audio_frequency = sum(audio_frequencies)/len(audio_frequencies)
avg_audio_frequency, avg_audio_frequency + rx_qrg

class FindSignalEdge:

    # num_of_samples 128 experimentally arrived at as yielding reasonable values
    def __init__(self, frequency, num_of_samples=128, samples_per_second=48000):
        self.num_of_samples = num_of_samples
        self.index_adjustment = round(num_of_samples * 1.66) # Experimentally arrived
        ts = numpy.array([n/samples_per_second for n in range(0, num_of_samples)])
        ωts = numpy.float64(2*math.pi*frequency) * ts
        self.cos_samples = numpy.cos(ωts)
        self.sin_samples = numpy.sin(ωts)
        self.edge_rise = numpy.array([1 if n < num_of_samples else -1 for n in range(0, 2*num_of_samples)]) 
        self.edge_fall = numpy.array([-1 if n < num_of_samples else 1 for n in range(0, 2*num_of_samples)]) 
        
    def convolve(self, samples):
        sin_cor = scipy.signal.convolve(self.sin_samples, samples, mode="full")
        cos_cor = scipy.signal.convolve(self.cos_samples, samples, mode="full")
        ampl2 = sin_cor**2 + cos_cor**2
        return ampl2 / self.num_of_samples

    def look_for_edge(self, samples, rise):
        if rise:
            return scipy.signal.convolve(self.edge_rise, self.convolve(samples)) / self.num_of_samples
        else:
            return scipy.signal.convolve(self.edge_fall, self.convolve(samples)) / self.num_of_samples

    def find_edge_raw(self, samples, rise):
        lfe = self.look_for_edge(samples, rise)
        # Homebrew peak finder (scipy has something better, should figure out how to use it).
        value, lfe_index = functools.reduce(lambda b,n: b if n[0] <= b[0] else n,
                                                    zip(lfe, range(0, len(lfe))))
        return lfe_index

    def find_edge(self, grab, from_t: datetime, to_t: datetime, rise=True) -> datetime:
        from_t_index = grab.sample_index(from_t)
        to_t_index = grab.sample_index(to_t)
        edge_index_raw = self.find_edge_raw(grab.samples[from_t_index:to_t_index], rise)
        edge_index_adjusted = edge_index_raw - self.index_adjustment
        return grab.datetime(from_t_index + edge_index_adjusted)

fse = FindSignalEdge(avg_audio_frequency)
real_starts = [
    fse.find_edge(grab, icmt+timedelta(seconds=50), icmt+timedelta(seconds=60), rise=True) \
       - timedelta(seconds=60)
    for icmt in into_the_carrier_minute_timestamps
]
real_starts

def oszi_data(from_t, to_t):
    from_i = grab.sample_index(from_t)
    to_i = grab.sample_index(to_t)
    frames = grab.samples[from_i:to_i]
    times = [grab.datetime(i) for i in range(from_i, to_i)]
    df = pandas.DataFrame(
        {
            "samples": frames,
        },
        index= [(t.second-60 if 30 < t.second else t.second) + t.microsecond/1e6 for t in times]
    )
    df.plot(grid=True, figsize=(12,6), xlabel="seconds into minute")
    return df

# Manually fished out of the data using the two lines below and experimentation.
impulse_minutes_starts = [
    datetime(2026, 3, 6, 21, 2, 58, 0, tzinfo=UTC),
    datetime(2026, 3, 6, 21, 5, 17, 0, tzinfo=UTC),
    datetime(2026, 3, 6, 21, 7, 35, 0, tzinfo=UTC),
    datetime(2026, 3, 6, 21, 9, 53, 0, tzinfo=UTC),
    datetime(2026, 3, 6, 21, 12, 11, 0, tzinfo=UTC),
    datetime(2026, 3, 6, 21, 14, 29, 0, tzinfo=UTC),
]
start = impulse_minutes_starts[0]
oszi_data(start+timedelta(seconds=0), start+timedelta(seconds=0.6))
pass

# There was something wrong with at least the second impulse batch. Document the problem:
start = impulse_minutes_starts[2] # + timedelta(seconds=1)
oszi_data(start-timedelta(seconds=1), start-timedelta(seconds=0.4))
oszi_data(start, start+timedelta(seconds=0.6))
oszi_data(start+timedelta(seconds=58), start + timedelta(seconds=58.6))
oszi_data(start+timedelta(seconds=59), start + timedelta(seconds=59.6))
# oszi_data(start - delta, start + timedelta(seconds=0.5) - delta)
pass

[(t1-t0).total_seconds() for t0, t1 in zip(impulse_minutes_starts[:-1],impulse_minutes_starts[1:])]

# Find the PPS delay for the minute i, for i in 0 to 5

def establish_pps_delay(i: int):
    # First step: Start investigation to find the impulse minute:
    start_investigation = impulse_minutes_starts[i]

    # They tend to goof on the last impulse train:
    def find_intervals(start_t, expected_num_of_impulsetrains=59):
        result = []
        for i in range(0, expected_num_of_impulsetrains):
            now_t = (start_t + timedelta(seconds=i))
            result.append(fse.find_edge(grab, now_t, now_t + timedelta(seconds=0.2)))
        return result

    impulse_starts = find_intervals(start_investigation)
    median_µs = sorted([i_start.microsecond for i_start in impulse_starts])[len(impulse_starts)//2]
    filtered_impulse_starts = [
        i_start for i_start in impulse_starts if median_µs-2000 < i_start.microsecond < median_µs+2000
    ]
    filtered_impulse_starts_df = pandas.DataFrame(
        {
            "ms": [i_start.microsecond / 1e3 for i_start in filtered_impulse_starts]
        },
        index = filtered_impulse_starts
    )
    avg_µs = sum((i_start.microsecond for i_start in filtered_impulse_starts)) / len(filtered_impulse_starts)
    filtered_impulse_starts_df.plot(kind="kde", grid=True, figsize=(12,6), xlabel="delay / ms")
    return median_µs, avg_µs

pps_median_avgs = [establish_pps_delay(i) for i in range(0,6)]
pps_median_avgs

impulse_delay_into_second = sum((ma[1] for ma in pps_median_avgs)) / len(pps_median_avgs)
impulse_delay_into_second * 1e-3 # ms

ar_af_ls = [
    find_rising_and_falling_edges(i, timedelta(seconds=0), timedelta(seconds=0.2))
    for i in range(0,6)
]
lengths = [(c[1]-c[0])*1e-3 for c in ar_af_ls]
sum(lengths)/len(lengths), lengths, ar_af_ls

ar_af_ls = [
    find_rising_and_falling_edges(i, timedelta(seconds=0.24), timedelta(seconds=0.35))
    for i in range(0,6)
]
lengths = [(c[1]-c[0])*1e-3 for c in ar_af_ls]
sum(lengths)/len(lengths), lengths, ar_af_ls

ar_af_ls = [
    find_rising_and_falling_edges(i, timedelta(seconds=0.35), timedelta(seconds=0.5))
    for i in range(0,6)
]
lengths = [(c[1]-c[0])*1e-3 for c in ar_af_ls]
sum(lengths)/len(lengths), lengths, ar_af_ls, sum((c[2] for c in ar_af_ls))

avg_audio_frequency + RF_FREQUENCY, 43.239+1.1, impulse_delay_into_second * 1e-3 - TOTAL_DELAY_FMTLAB * 1e3

# This is the only change after the hand-in of my results:

fmtlab_official_frequency = 3565712.72
print(f"Frequency error:    {3565712.536 - fmtlab_official_frequency:.3f} Hz")
fmtlab_official_pps_offset = 11.4
print(f"PPS offset error:    {12.159 - fmtlab_official_pps_offset:.3f} ms")
fmtlab_official_pulse_length = 44.3
print(f"Pulse length error: {44.239 - fmtlab_official_pulse_length:.3f} ms")