Time-domain phase fitting works almost

Except that the initial guess seems to massively impact the fitted phase.
If the initial_phase is submitted, it seems to fit quite fine

fourier/05_sine_fitting.py

@@ -9,20 +9,22 @@ Then fit in t-domain to resolve \\varphi_0
 """
 
 import matplotlib.pyplot as plt
-import matplotlib.gridspec as gridspec
 import numpy as np
 if not True:
     import numpy.fft as ft
 else:
     import scipy.fftpack as ft
 import scipy.optimize as opt
+from scipy.signal import hilbert
+
 
 from mylib import *
 
 rng = np.random.default_rng()
 
 def guess_sine_parameters(samples, fft=None, fft_freqs=None, guess=[None,None,None,None]):
-    """Use crude methods to guess the parameters to a sine wave
+    """
+    Use crude methods to guess the parameters to a sine wave
     from properties of both samples and their fourier transform.
 
     Parameters:
@@ -31,13 +33,13 @@ def guess_sine_parameters(samples, fft=None, fft_freqs=None, guess=[None,None,No
 
     guess - arraylike or float or None
         If float, this is interpreted as a frequency
-        Order of parameters: [amplitude, frequency, phase, offset]
+        Order of parameters: [amplitude, frequency, phase, baseline]
         If one parameter is None, it is filled with an approximate value if available.
 
     Returns:
     -----------
     guess - arraylike
-        An updated version of init_guess: [amplitude, frequency, phase, offset]
+        An updated version of init_guess: [amplitude, frequency, phase, baseline]
     """
 
     if not hasattr(guess, '__len__'):
@@ -68,30 +70,20 @@ def guess_sine_parameters(samples, fft=None, fft_freqs=None, guess=[None,None,No
             continue
 
         if i == 0: # amplitude
-            guess[i] = np.std(samples) * (2 ** 1/2)
+            if False:
+                guess[i] = np.std(samples) * (2 ** 1/2)
+            else:
+                guess[i] = max(samples-np.mean(samples))
         elif i == 1: # frequency
             guess[i] = nearest_f
         elif i == 2: # phase
             guess[i] = nearest_phase
-        elif i == 3: # offset   samples
+        elif i == 3: # baseline   samples
             guess[i] = np.mean(samples)
 
     return guess
 
-def curvefit_sine(time, samples, init_guess, fitfunc=sine_fitfunc, bounds=(-np.inf, np.inf), **curve_kwargs):
+def fit_sine_to_samples(time, samples, samplerate=1, bandpass=None, guess=[None,None,None,None], fitfunc=sine_fitfunc, fft=None, freqs=None, bounds=None, restrained_fit=False, **curve_kwargs):
-    """Fit a sine to samples with a supplied initial guess.
-
-    This function sets bounds for the phase.
-    """
-    if bounds is None or bounds == (-np.inf, np.inf):
-        high_bounds = np.array([np.inf, np.inf, +1*np.pi, np.inf])
-        low_bounds = -1*high_bounds
-
-        bounds = (low_bounds, high_bounds)
-
-    return opt.curve_fit(fitfunc, time, samples, p0=init_guess, bounds=bounds, **curve_kwargs)
-
-def fit_sine_to_samples(time, samples, samplerate=1, bandpass=None, guess=[None,None,None,None], fitfunc=sine_fitfunc, fft=None, freqs=None,**curve_kwargs):
     if bandpass is not None or guess[1] is None or guess[2] is None:
         if fft is None:
             fft = ft.rfft(samples)
@@ -102,14 +94,64 @@ def fit_sine_to_samples(time, samples, samplerate=1, bandpass=None, guess=[None,
             fft[(freqs < bandpass[0]) | (freqs > bandpass[1])] = 0
             samples = ft.irfft(fft, samples.size)
 
-    old_guess = guess.copy()
     guess = guess_sine_parameters(samples, fft=fft, fft_freqs=freqs, guess=guess)
 
+    guess = np.array(guess)
+
+    if restrained_fit:
+        # Restrained fit
+        # only allow phase to be fitted
+        # Take the amplitude from the hilbert envelope of the (bandpassed) samples
+
+        # References for lambda
+
+        frequency = guess[1]
+        baseline = guess[3]
+        envelope = np.abs(hilbert(samples))
+        base_fitfunc = fitfunc
+
+        samples = samples/envelope
+
+        fitfunc = lambda t, amplitude, phase: base_fitfunc(t, amp=amplitude, phase=phase, freq=frequency, baseline=baseline)
+
+        old_guess = guess.copy()
+
+        guess = guess[[0,2]]
+
+        if bounds is None:
+            sample_max = max(samples)
+
+            low_bounds = np.array([0.8,-np.pi])
+            high_bounds = np.array([1.2, np.pi])
+        else:
+            low_bounds =  bounds[0][[0,2]]
+            high_bounds = bounds[1][[0,2]]
+
+        bounds = (low_bounds, high_bounds)
+
+    elif bounds is None :
+        high_bounds = np.array([np.inf, np.inf, +1*np.pi, np.inf])
+        low_bounds = -1*high_bounds
+
+        bounds = (low_bounds, high_bounds)
+
+    print(bounds, guess)
+
     try:
-        fit = curvefit_sine(time, samples, guess, fitfunc=fitfunc, **curve_kwargs)
+        fit = opt.curve_fit(fitfunc, time, samples, p0=guess, bounds=bounds, **curve_kwargs)
     except RuntimeError:
         fit = None
 
+    if len(bounds[0]) == 1 or restrained_fit:
+        # Restrained fitting was used
+        # merge back into guess and fit
+
+        guess = old_guess
+        fit = [
+                np.array([fit[0][0], old_guess[1], fit[0][1], old_guess[3]]),
+                fit[1]
+            ]
+
     return fit, guess, (fft, freqs, samples)
 
 def chi_sq(observed, expected):
@@ -125,7 +167,8 @@ def simulate_noisy_sine_fitting_SNR_and_residuals(
         N=1, snr_band=passband(), noise_band=passband(),
         t_length=1e-6, f_sample=250e6,
         noise_sigma=1, init_params=[1, 50e6, None, 0],
-        show_original_signal_figure=True, show_bandpassed_signal_figure=True
+        show_original_signal_figure=False, show_bandpassed_signal_figure=True,
+        restrained_fit=True
         ):
     residuals = np.empty( (int(N), len(init_params)) )
     real_snrs = np.empty( (int(N)) )
@@ -137,12 +180,13 @@ def simulate_noisy_sine_fitting_SNR_and_residuals(
             print("Iteration {} running".format(j))
 
         # set random phase
-        init_params[2] = 2*np.pi*rng.random()
+        init_params[2] = phasemod(2*np.pi*rng.random())
 
         samples = sine_fitfunc(time, *init_params)
-        noise = None
         if noise_sigma: # noise
             noise = rng.normal(0,noise_sigma, size=(len(samples)))
+        else:
+            noise = np.zeros(len(samples))
 
         real_snrs[j] = signal_to_noise(samples, noise, signal_band=snr_band, samplerate=f_sample, noise_band=noise_band)
 
@@ -162,7 +206,13 @@ def simulate_noisy_sine_fitting_SNR_and_residuals(
 
             axs1[1].legend()
 
-        fit, guess, (fft, freqs, bandpassed) = fit_sine_to_samples(time, samples+noise, f_sample, guess=[None,f_sine,None,None], bandpass=snr_band)
+        if False:
+            # use initial_params as guess
+            guess = init_params
+        else:
+            guess = [None, f_sine, None, None]
+        fit, guess, (fft, freqs, bandpassed) = fit_sine_to_samples(time, samples+noise, f_sample, guess=guess, bandpass=snr_band, restrained_fit=restrained_fit)
+
 
         if fit is None:
             residuals[j] = np.nan
@@ -170,27 +220,75 @@ def simulate_noisy_sine_fitting_SNR_and_residuals(
 
         residuals[j] = normalise_sine_params(init_params - fit[0])
 
+        # figures
         if show_bandpassed_signal_figure and (j==0 or N == 1):
-            fitted_sine = sine_fitfunc(time, *fit[0])
+            analytic_signal = hilbert(bandpassed)
+            envelope = np.abs(analytic_signal)
+            instant_phase = np.angle(analytic_signal)
 
-            fig2, axs2 = plot_signal_and_spectrum(
+            fit_params = fit[0].tolist()
-                    bandpassed, f_sample, "Bandpassed samples\nS/N:{:.2e}".format(real_snrs[j]),
+            fit_params[0] = envelope
-                    freq_unit='MHz', freq_scaler=freq_scaler
+            fitted_sine = sine_fitfunc(time, *fit_params)
-                    )
-            for ax in axs2[[1,2]]:
-                ax.axvline(f_sine/freq_scaler, color='r', alpha=0.4) # f_beacon
-                ax.axvspan(snr_band[0]/freq_scaler,snr_band[1]/freq_scaler, color='purple', alpha=0.3, label='signalband') # snr
-                ax.axvspan(noise_band[0]/freq_scaler, noise_band[1]/freq_scaler, color='orange', alpha=0.3, label='noiseband') # noise_band
 
-            l = axs2[0].plot(time, fitted_sine, label='fit')
-            axs2[0].text(1, 1, '$\chi/d.o.f. = {:.2e}/{:.2e}$'.format(chi_sq(fitted_sine, samples), dof(samples,4)), transform=axs2[0].transAxes, ha='right', va='top')
 
-            # indicate initial phase
+            if False:
-            axs2[2].axhline(init_params[2], color='r', alpha=0.4)
+                fig4, axs4 = plt.subplots(2,1, sharex=True)
-            axs2[2].axhline(init_params[2], color=l[0].get_color(), alpha=0.4)
+                fig4.suptitle("Bandpassed Hilbert")
+                axs4[1].set_xlabel("Time")
 
-            axs2[0].legend(loc='upper left')
+                axs4[0].set_ylabel("Instant Phase")
-            axs2[1].legend()
+                axs4[0].plot(time, instant_phase, marker='.')
+                #axs4[0].axhline(init_params[2], color='r')
+
+
+                axs4[1].set_ylabel("Instant Freq")
+                axs4[1].plot(time[1:], np.diff(np.unwrap(instant_phase)) / (2*np.pi*f_sample), marker='.')
+                #axs4[1].axhline(init_params[1], color='r')
+
+
+            ## Next figure
+            if True:
+                fig2, axs2 = plot_signal_and_spectrum(
+                        bandpassed, f_sample, "Bandpassed samples\nS/N:{:.2e}".format(real_snrs[j]),
+                        freq_unit='MHz', freq_scaler=freq_scaler,
+                        signal_kwargs=dict(alpha=0.8, time_unit='us')
+                        )
+                for ax in axs2[[1,2]]:
+                    ax.axvline(f_sine/freq_scaler, color='r', alpha=0.4) # f_beacon
+                    ax.axvspan(snr_band[0]/freq_scaler,snr_band[1]/freq_scaler, color='purple', alpha=0.3, label='signalband') # snr
+                    ax.axvspan(noise_band[0]/freq_scaler, noise_band[1]/freq_scaler, color='orange', alpha=0.3, label='noiseband') # noise_band
+
+                l = axs2[0].plot(time, fitted_sine, label='fit', alpha=0.8)
+                #axs2[0].text(1, 1, '$\chi/d.o.f. = {:.2e}/{:.2e}$'.format(chi_sq(fitted_sine, samples), dof(samples,4)), transform=axs2[0].transAxes, ha='right', va='top')
+
+                axs2[0].plot(time, envelope, label='envelope')
+
+                # indicate initial phase
+                axs2[2].axhline(init_params[2], color='r', alpha=0.4)
+                axs2[2].axhline(fit[0][2], color=l[0].get_color(), alpha=0.4)
+
+                axs2[0].legend(loc='upper left')
+                axs2[1].legend()
+
+
+            if True:
+                fig5, axs5 = plt.subplots(2,1, sharex=True)
+                fig5.suptitle("Bandpassed Samples vs Model")
+                axs5[0].set_ylabel("Amplitude")
+                axs5[0].plot(bandpassed, label='samples', alpha=0.8)
+                axs5[0].plot(fitted_sine, label='fit', alpha=0.8)
+                axs5[0].plot(envelope, label='envelope')
+
+                axs5[0].plot(samples, label='orig sine', alpha=0.8)
+
+                axs5[0].legend()
+
+                axs5[1].set_ylabel("Residuals")
+                axs5[1].set_xlabel("Sample")
+                axs5[1].plot(samples - fitted_sine, label="Sine - Model", alpha=0.8)
+                axs5[1].plot(bandpassed - fitted_sine, label="Bandpassed - Model", alpha=0.8)
+
+                axs5[1].legend()
 
             print("init:", init_params)
             print("fit :", fit[0])
@@ -203,32 +301,34 @@ if __name__ == "__main__":
     from argparse import ArgumentParser
     from myscriptlib import save_all_figs_to_path_or_show
 
-    rng = np.random.default_rng(1)
-
     parser = ArgumentParser(description=__doc__)
     parser.add_argument("fname", metavar="path/to/figure[/]", nargs="?", help="Location for generated figure, will append __file__ if a directory. If not supplied, figure is shown.")
     parser.add_argument("-n", "--n-rand", dest='N', default=1, type=int, nargs='?', help='Number of random sines to fit')
+    parser.add_argument('--seed', default=1, type=int, help='RNG seed')
     args = parser.parse_args()
     default_extensions = ['.pdf', '.png']
 
     if args.fname == 'none':
         args.fname = None
 
+    rng = np.random.default_rng(args.seed)
+
     report_N_nan = True
+    restrained_fitting = True
 
     f_sine = 53.123456 # MHz
-    sine_amplitude = 0.2
+    sine_amplitude = 1
-    sine_offset = 0
+    sine_baseline = 0
-    init_params = np.array([sine_amplitude, f_sine, None, sine_offset])
+    init_params = np.array([sine_amplitude, f_sine, None, sine_baseline])
 
     N = int(args.N)
     f_sample = 250 # MHz
     t_length = 10 # us
-    noise_sigma = 1
+    noise_sigma = 0.01
 
     f_delta = 1/t_length
     noise_band = (30,80) # MHz
-    snr_band = (f_sine -2*f_delta, f_sine + 2*f_delta)
+    snr_band = (f_sine -50*f_delta, f_sine + 50*f_delta)
 
     time = sampled_time(f_sample, end=t_length)
 
@@ -236,7 +336,7 @@ if __name__ == "__main__":
 
     ###### End of inputs
 
-    residuals, real_snrs, _ = simulate_noisy_sine_fitting_SNR_and_residuals(N=N, snr_band=snr_band, noise_band=noise_band, t_length=t_length, f_sample=f_sample, noise_sigma=noise_sigma, init_params=init_params)
+    residuals, real_snrs, _ = simulate_noisy_sine_fitting_SNR_and_residuals(N=N, snr_band=snr_band, noise_band=noise_band, t_length=t_length, f_sample=f_sample, noise_sigma=noise_sigma, init_params=init_params, restrained_fit=restrained_fitting)
 
     # Filter NaNs from fit attempts that failed
     nan_mask = ~np.isnan(residuals).any(axis=1)
@@ -248,26 +348,43 @@ if __name__ == "__main__":
     real_snrs = real_snrs [ nan_mask ]
 
     ## Plot Signal-to-Noise vs Residuals of the fit paramters
-    fig, axs = plt.subplots(1,4, sharey=True)
+    fig, axs = plt.subplots(1,1 + 2*( not restrained_fitting), sharey=True)
-    fig.suptitle("S/N vs Residuals, S/N Band ({:.2e},{:.2e})MHz".format(snr_band[0]/freq_scaler, snr_band[-1]/freq_scaler))
+
+    if not hasattr(axs,'__len__'):
+        axs = [axs]
+
+    fig.suptitle("S/N vs Residuals\nS/N Band ({:.2e},{:.2e})MHz \namp/sigma: {}".format(snr_band[0]/freq_scaler, snr_band[-1]/freq_scaler, sine_amplitude/ noise_sigma))
     axs[0].set_ylabel("S/N")
+    j = 0 # plot counter
     for i in range(len(init_params)):
+        if restrained_fitting and i in [0,1,3]:
+            continue
+
         unit_scaler = [1, 1][i==1]
         unit_string = ['', '[MHz]'][i==1]
-        xlabel = ["Amplitude", "Frequency", "Phase", "Offset"][i]
+        xlabel = ["Amplitude", "Frequency", "Phase", "Baseline"][i]
 
-        axs[i].set_xlabel(xlabel + unit_string)
+        if i == 2:
-        axs[i].plot(residuals[:,i]/unit_scaler, real_snrs, ls='none', marker='o')
+            #axis_pi_ticker(axs[j].xaxis)
+            axs[j].set_xlim(-np.pi, np.pi)
+
+
+        real_snrs[np.isnan(real_snrs)] = 1 # Show nan values
+
+        axs[j].set_xlabel(xlabel + unit_string)
+        axs[j].plot(residuals[:,i]/unit_scaler, real_snrs, ls='none', marker='o', alpha=max(0.3, 1/len(real_snrs)))
+
+        j += 1
 
     ## Plot Histograms of the Residuals
     if True and N > 1:
         for j in range(len(init_params)):
-            if j == 0 or j == 3:
+            if j == 3 or restrained_fitting and j == 1 or j == 0:
                 continue
 
             unit_scaler = [1, freq_scaler][j==1]
             unit_string = ['', '[MHz]'][j==1]
-            xlabel = ["Amplitude", "Frequency", "Phase", "Offset"][j]
+            xlabel = ["Amplitude", "Frequency", "Phase", "Baseline"][j]
 
             title = xlabel + " residuals"
             title += "\n"
@@ -290,7 +407,7 @@ if __name__ == "__main__":
                 maj_div = max(1, 2**np.ceil(np.log2(np.pi/(xmax-xmin)) + 1 ))
                 min_div = maj_div*12
 
-                axis_pi_ticker(ax.xaxis, major_divider=maj_div, minor_divider=min_div)
+                #axis_pi_ticker(ax.xaxis, major_divider=maj_div, minor_divider=min_div)
 
         # Plot histogram between phase and frequency
         if True and N > 10:

fourier/mylib/util.py

@@ -16,12 +16,12 @@ def phasemod(phase, low=np.pi):
 # Alias phase_mod to phasemod
 phase_mod = phasemod
 
-def sine_fitfunc(t, amp=1, freq=1, phase=0, off=0, t_delay=0):
+def sine_fitfunc(t, amp=1, freq=1, phase=0, baseline=0, t_delay=0):
     """Simple sine wave for fitting purposes"""
-    return amp*np.cos( 2*np.pi*freq*(t-t_delay) + phase) + off
+    return amp*np.cos( 2*np.pi*freq*(t-t_delay) + phase) + baseline
 
 def sin_delay(f, t, phase=0):
-    return sine_fitfunc(t, amp=1, freq=f, phase=phase, off=1, t_delay=0)
+    return sine_fitfunc(t, amp=1, freq=f, phase=phase, baseline=1, t_delay=0)
 
 def sampled_time(sample_rate=1, start=0, end=1, offset=0):
     return offset + np.arange(start, end, 1/sample_rate)