Demodulation of climate index using annual, etc carrier, optimized for impulse train

extract_satellite_periods_harmonics.py

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

def load_signal(csv_path):
    df = pd.read_csv(csv_path)
    time = df.iloc[:, 0].values
    signal = df.iloc[:, 1].values
    return time, signal

def impulse_train_multiplier(time, carrier_freq, num_harmonics):
    """
    Generates a real impulse train approximation: sum of cosines at harmonics of carrier_freq.
    num_harmonics is the maximum |k| (excluding k=0 if undesired).
    """
    mult = np.zeros_like(time)
    for k in range(-num_harmonics, num_harmonics + 1):
        if k == 0:
            continue  # skip DC if desired, remove this line to include DC
        mult += np.cos(2 * np.pi * k * carrier_freq * time)
    return mult

def periodogram(time, signal, freq_min, freq_max, num_freqs):
    T = time[-1] - time[0]
    time_shifted = time - time[0]
    freqs = np.linspace(freq_min, freq_max, num_freqs)
    energies = []
    for f in freqs:
        carrier = np.exp(-1j * 2 * np.pi * f * time_shifted)
        product = signal * carrier
        energy = np.abs(product.sum())**2
        energies.append(energy)
    return freqs, np.array(energies)

def annotate_top_peaks(ax, freqs, energies, num_peaks=50, min_distance=5):
    """
    Annotate the chart with callouts for the top N peaks.
    min_distance: minimum index separation between peaks.
    """
    from scipy.signal import find_peaks

    # Find all peaks
    peaks, _ = find_peaks(energies, distance=min_distance)
    # Get top N by energy
    if len(peaks) > num_peaks:
        top_peaks = peaks[np.argsort(energies[peaks])[-num_peaks:]]
        # Sort top peaks by frequency for annotation order
        top_peaks = top_peaks[np.argsort(freqs[top_peaks])]
    else:
        top_peaks = peaks

    for idx in top_peaks:
        freq = freqs[idx]
        energy = energies[idx]
        ax.annotate(
            f"{freq:.2f}",
            xy=(freq, energy),
            xytext=(freq, energy * 1.05),
            arrowprops=dict(arrowstyle="->", color="red", lw=1),
            fontsize=8,
            color="red",
            ha='center'
        )

def extract_satellite_periods_with_harmonics(
    csv_path, carrier_freq, freq_span, num_harmonics=10, num_freqs=1024, plot=True, num_peaks=50
):
    """
    Main function to extract and visualize satellite periods using a carrier impulse train.
    Annotates the top peaks in the periodogram.
    """
    time, signal = load_signal(csv_path)
    multiplier = impulse_train_multiplier(time, carrier_freq, num_harmonics)
    signal_mixed = signal * multiplier
    freqs, energies = periodogram(
        time, signal_mixed,
        freq_min=-freq_span/2, freq_max=freq_span/2, num_freqs=num_freqs
    )
    if plot:
        plt.figure(figsize=(12, 6))
        ax = plt.gca()
        plt.plot(freqs, energies)
        plt.title(f'Periodogram with Harmonic Impulse Train (±{num_harmonics} harmonics)')
        plt.xlabel('Frequency Offset (1/yr)')
        plt.ylabel('Energy')
        plt.grid()
        annotate_top_peaks(ax, freqs, energies, num_peaks=num_peaks)
        plt.tight_layout()
        plt.show()
    return freqs, energies

if __name__ == "__main__":
    import argparse
    parser = argparse.ArgumentParser(description="Extract satellite periods using a carrier impulse train (harmonics).")
    parser.add_argument('csv_path', help='Path to CSV file with time,value columns.')
    parser.add_argument('carrier_freq', type=float, help='Carrier frequency (1/yr).')
    parser.add_argument('freq_span', type=float, help='Frequency span around carrier (1/yr).')
    parser.add_argument('--num_harmonics', type=int, default=10, help='Number of harmonics in impulse train.')
    parser.add_argument('--num_freqs', type=int, default=1024, help='Number of frequency points.')
    parser.add_argument('--num_peaks', type=int, default=50, help='Number of top peaks to annotate.')
    parser.add_argument('--no-plot', action='store_true', help="Don't plot the result.")
    args = parser.parse_args()
    extract_satellite_periods_with_harmonics(
        args.csv_path,
        args.carrier_freq,
        args.freq_span,
        num_harmonics=args.num_harmonics,
        num_freqs=args.num_freqs,
        plot=not args.no_plot,
        num_peaks=args.num_peaks
    )

2 year train, 1/y wide spectrum
> python .\extract_satellite_periods_harmonics.py nino4.csv 0.5 1.0

comparison with Drac/Anom model
NINO34

NINO4