#!/usr/bin/env python3
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import sys
from datetime import datetime
import subprocess
import os

csvfile = "/home/x/Desktop/sweep.csv"
reference_file = "/home/x/Desktop/reference.csv"

#Do the sweep through the bash command
# Build the rtl_power command
rtl_cmd = [
    "rtl_power",
    "-f", "50M:1760M:125k", #startfreq, stopfreq, bandwidth eg 125k
    "-i", "0.5",   # integration interval 0.5s
    "-g", "30", # gain
    "-1", #fast sweep
    csvfile
]

# Run rtl_power
print("Starting sweep...")
subprocess.run(rtl_cmd)
print(f"Sweep finished. Results in {csvfile}")

# Read first line to get header info (some rtl_power versions have no header)
# Use pandas to parse; skip initial bad rows if present
df = pd.read_csv(csvfile, header=None)

# Extract all low/high ranges and flatten the data
freqs = []
powers = []

for i in range(len(df)):
    lowHz = df.iloc[i, 2]
    highHz = df.iloc[i, 3]
    stepHz = df.iloc[i, 4]
    pwr = df.iloc[i, 6:].astype(float).values.flatten()
    # build freq axis for this row
    f = lowHz + np.arange(len(pwr)) * stepHz
    freqs.extend(f)
    powers.extend(pwr)

freqs = np.array(freqs) / 1e6  # MHz
powers = np.array(powers)

# Sort by frequency (just in case)
idx = np.argsort(freqs)
freqs = freqs[idx]
powers = powers[idx]

if os.path.exists(reference_file):
    ref_df = pd.read_csv(reference_file, sep=",")
    
    # Compute center frequencies (in MHz)
    ref_freqs = ((ref_df["start"] + ref_df["end"]) / 2) / 1e6
    ref_power = ref_df["median"].values

    # Interpolate reference powers to match current sweep frequencies
    ref_interp = np.interp(freqs, ref_freqs, ref_power)

    # Compute difference (current - reference)
    diff = powers - ref_interp

    # Plot reference vs current
    plt.figure(figsize=(12,5))         
    plt.plot(freqs, powers, color="red", label="Current Sweep", linewidth=1.2, alpha=0.7)
    plt.plot(ref_freqs, ref_power, color="black", label="Reference Sweep", linewidth=2)
    plt.xlabel("Frequency (MHz)")
    plt.ylabel("Power (dB)")
    plt.title("Current vs Reference Spectrum")
    plt.legend()
    plt.grid(True)
    plt.tight_layout()
    plt.show()

    # Plot difference (ΔdB)
    plt.figure(figsize=(12,4))
    plt.plot(freqs, diff, color="red", linewidth=1)
    plt.axhline(0, color="black", linestyle="--", linewidth=0.8)
    plt.xlabel("Frequency (MHz)")
    plt.ylabel("Power Difference (dB)")
    plt.title("Difference: Current Sweep − Reference Sweep")
    plt.grid(True)
    plt.tight_layout()
    plt.show()

    # Optional: print significant deviations
    threshold = 5  # dB
    sig_idx = np.where(diff > threshold)[0]
    if len(sig_idx) > 0:
        print("\nSignificant deviations (> +5 dB):")
        for i in sig_idx:
            print(f"{freqs[i]:8.3f} MHz  Δ = {diff[i]:6.2f} dB")
    else:
        print("\nNo significant deviations from reference.")
else:
    print(f"\n Reference file not found: {reference_file}")