From a32ae6b2efdf7c759dbfd3b3be96fbd8858287e7 Mon Sep 17 00:00:00 2001 From: Eric Teunis de Boone Date: Tue, 14 Nov 2023 15:22:45 +0100 Subject: [PATCH] PDFs: modified script from Harm for PhasorSum distributions --- simulations/12_noise_phase.py | 154 ++++++++++++++++++++++++++++++++++ 1 file changed, 154 insertions(+) create mode 100644 simulations/12_noise_phase.py diff --git a/simulations/12_noise_phase.py b/simulations/12_noise_phase.py new file mode 100644 index 0000000..6b403b7 --- /dev/null +++ b/simulations/12_noise_phase.py @@ -0,0 +1,154 @@ +#!/usr/bin/env python3 + +__doc__ = """ +Phase and Amplitude distributions for a phasor in the presence of noise. + +Author: Harm Schoorlemmer +""" + +import numpy as np +import matplotlib.pyplot as plt +import matplotlib +import scipy.stats as stat +from scipy import special + +from lib.util import MethodMappingProxy as MethodProxy + +def expectation(x,pdfx): + dx = x[1]-x[0] + return np.sum(x*pdfx*dx) + +def variance(x,pdfx): + mu = expectation(x,pdfx) + dx = x[1]-x[0] + return np.sum((x**2*pdfx*dx))-mu**2 + +def phase_distribution(theta,sigma,s): + theta = np.asarray(theta) + ct = np.cos(theta) + st = np.sin(theta) + k=s/sigma + pipi=2*np.pi + return (np.exp(-k**2/2)/pipi) + ( + (pipi**-0.5)*k*np.exp(-(k*st)**2/2)) * ( + (1.+special.erf(k*ct*2**-0.5))*ct/2) + +def phase_distribution_gauss(theta,sigma,s): + theta = np.asarray(theta) + k=s/sigma + return (2*np.pi)**-0.5*k*np.exp(-(k*theta)**2/2) + +def amplitude_distribution(a,sigma,s): + a = np.asarray(a) + k = s/sigma + return stat.rice.pdf(a,s/sigma,scale=sigma) + +def amplitude_distribution_gauss(a,sigma,s): + a = np.asarray(a) + k=s/sigma + return (2*np.pi)**-0.5*np.exp(-((a-s)/sigma)**2/2) + +signal_max= 4 +amp_max = signal_max*2 +thetas = np.linspace(-np.pi,np.pi,500) +amplitudes = np.linspace(0,amp_max,500) +signals = np.linspace(0,signal_max,5) +sigma = 1 + + +## figure 1 +if True: + fig, ax = plt.subplots(1,2,figsize=(2*8,1*8)) + _fig1, _ax0 =plt.subplots(1,1, figsize=(1*8, 1*8)) + _fig2, _ax1 =plt.subplots(1,1, figsize=(1*8, 1*8)) + + ax0 = MethodProxy(ax[0], _ax0) + ax1 = MethodProxy(ax[1], _ax1) + for s in signals: + pdfs_label='s/$\sigma$ ='+str(s) + phase_vals= phase_distribution(thetas,sigma,s) + amp_vals= amplitude_distribution(amplitudes,sigma,s) + phase_vals_g = phase_distribution_gauss(thetas,sigma,s) + + ax0.plot(amplitudes,amp_vals, label=pdfs_label) + ax1.plot(thetas,phase_vals, label=pdfs_label) + ax0.legend() + _ax1.legend()# only in the separate figure + ax0.set_xlabel(r'$a$') + ax0.set_xlabel(r'$\theta$') + ax1.set_ylabel(r'$p(a)$') + ax1.set_ylabel(r'$p(\theta)$') + # ax[0].grid() + # ax[1].grid() + MethodProxy(fig, _fig1, _fig2).tight_layout() + + fig.savefig('pdfs.pdf') + _fig1.savefig('pdfs-amplitudes.pdf') + _fig2.savefig('pdfs-phases.pdf') + + plt.close(_fig1) + plt.close(_fig2) + +## figure 2 +amplitudes = np.linspace(0,amp_max*5,500) +signals = np.linspace(0.1,signal_max*5,101) +if False: + fig2, ax2 = plt.subplots(2,2,figsize=(2*8,2*8)) + ax2 = fig2.get_axes() + V_theta = [variance(thetas,phase_distribution(thetas,sigma,s)) for s in signals ] + E_theta=[expectation(thetas,phase_distribution(thetas,sigma,s)) for s in signals ] + V_theta_g = [variance(thetas,phase_distribution_gauss(thetas,sigma,s)) for s in signals ] + E_theta_g=[expectation(thetas,phase_distribution_gauss(thetas,sigma,s)) for s in signals ] + V_a = [variance(amplitudes,amplitude_distribution(amplitudes,sigma,s)) for s in signals ] + E_a=[expectation(amplitudes,amplitude_distribution(amplitudes,sigma,s)) for s in signals ] + V_a_g = [variance(amplitudes,amplitude_distribution_gauss(amplitudes,sigma,s)) for s in signals ] + E_a_g=[expectation(amplitudes,amplitude_distribution_gauss(amplitudes,sigma,s)) for s in signals ] + + ax2[0].plot(signals,E_a,label='$p(a)$') + ax2[0].plot(signals,E_a_g,ls='dashed',label='Gaussian approx.') + ax2[0].set_xscale('log') + ax2[0].set_yscale('log') + ax2[0].set_ylabel('$\mu_a$') + + ax2[1].plot(signals,V_a,label='$p(a)$') + ax2[1].plot(signals,V_a_g,ls='dashed',label='Gaussian approx.') + ax2[1].set_xscale('log') + ax2[1].set_ylabel('$\sigma_a^2$') + + ax2[2].plot(signals,E_theta,label=r'$p(\theta)$') + ax2[2].plot(signals,E_theta_g,ls='dashed',label='Gaussian approx.') + ax2[2].set_xscale('log') + ax2[2].set_ylim(-1.1,1.1) + ax2[2].set_ylabel(r'$\mu_\theta$') + + ax2[3].plot(signals,V_theta,label=r'$p(\theta)$') + ax2[3].plot(signals,V_theta_g,ls='dashed',label='Gaussian approx.') + ax2[3].set_xscale('log') + ax2[3].set_yscale('log') + ax2[3].set_ylabel(r'$\sigma_\theta^2$') + for a in ax2: + a.grid(which='both') + a.set_xlabel(r'$s/\sigma$') + a.legend() + + fig2.tight_layout() + fig2.savefig('expectation_variance.pdf') + +## figure 3, beacon timing accuracy +if True: + fig3, ax3 = plt.subplots(1,1,figsize=(1*8,1*8)) + ax3 = fig3.get_axes() + sigma_t = [variance(thetas,phase_distribution(thetas,sigma,s)) for s in signals ] + lfs=np.linspace(np.log10(50.),4,1) + for lf in lfs: + freq = (10**lf)*1e6 + sigma_t = np.asarray(sigma_t)**0.5/(2*np.pi*freq) + ax3[0].plot(signals,sigma_t/1e-9,'o-') + ax3[0].set_ylim(0,2.5) + ax3[0].set_xlabel(r'$s/\sigma$') + ax3[0].set_ylabel(r'$\Delta t$(ns)') + + fig3.tight_layout() + fig3.savefig('timing_accuracy.pdf') + +plt.show()